Integrating Hydrogeological and Microbiological Data and Modelling to Characterize the Hydraulic Features and Behaviour of Coastal Carbonate Aquifers: A Case in Western Cuba

: Carbonate aquifers are the primary source of freshwater in Cuba. Unfortunately, coastal groundwater is often contaminated by seawater intrusion. The main aim of the present study was to test the e ﬃ cacy of an experimental modelling approach, ranging from hydrogeology / geomorphology to microbiology, to better characterise both the hydraulic features and behaviour of a coastal carbonate aquifer and acquire useful information to prevent groundwater salinization. The interdisciplinary approach was an e ﬀ ective tool in order to understand (i) the hydraulic role played by some fault zones; (ii) the inﬂuence of discontinuous heterogeneities on groundwater ﬂow and saltwater wedge shape; (iii) mixing processes between di ﬀ erent water bodies (groundwater, surface water, seawater); (iv) the role of karst conduits in inﬂuencing the step-like halocline within the mixing zone between fresh groundwater and seawater.


Introduction
The Cuban land territory encompasses 110,992 km 2 , 66% of which are carbonate aquifers [1]. These aquifers are the primary source of freshwater and therefore are essential to the future needs of Cuba [2]. Unfortunately, in the Guanahacabibes Peninsula (Western Cuba), the contamination of the fresh groundwater by saline seawater may be intensified by different natural and anthropogenic factors. The main factors are related to the climatic regime, sea tides, the conditions of recharge of the aquifer(s), excessive aquifer exploitation by pumping for water supply and agricultural activities [3][4][5]. Therefore, the implementation of suitable monitoring and protective actions is fundamental for the preservation of these aquifers and for assuring the future use of this resource.
Coastal carbonate aquifers often have hydraulic links to the sea resulting in dominant conduit flow conditions, submarine freshwater springs and/or natural seawater intrusion into the aquifer through karst conduits [6,7]. Groundwater flow and exchanges between aquifers and seawater are influenced by more or less complex karst networks, often characterised by different levels of karstified limestone, due to many base-level variations imposed by relative changes in the sea level. However, the evolution and the hydrogeological behaviour of these coastal systems are different according to whether the sea level was characterised by rising or fall phases. According to Fleury et al. 2007 [7], during marine regression, carbonate dissolution increases due to the resulting elevation difference from the new sea level. Conversely, during marine transgression two different cases can be observed: (i) the submerged karst conduits continue to discharge the groundwater, generating submarine springs, or (ii) they let in seawater. On the whole, the existence of complex karst networks and the possible hydraulic interconnections between the conduits and the surrounding fissured rocks cause coastal carbonate aquifers to be highly heterogeneous, and more complex than porous or fractured aquifers. In that scenario, the zones characterised by higher hydraulic conductivity influence the distribution of the mixing zone, and both the morphology and the geometry of karst conduits are factors of utmost importance in controlling groundwater temperature and electrical conductivity gradients in the saturated zone [8].
At the study site, Peninsula de Guanahacabibes, Western Cuba (see Figure 1), hydrochemical and isotopic investigations were used in a previous paper [9] to analyse in detail the geochemical processes characterising the local interaction and mixing between rain-, surface-, ground-and seawater. The main aim of the present study was to test the efficacy of an experimental modelling approach, ranging from hydrogeology/geomorphology to microbiology, to better characterise both the hydraulic features and behaviour of a coastal carbonate aquifer, therefore acquiring useful information to prevent groundwater salinization.

Hydraulic Head Measurements
The hydraulic head was measured monthly at site 1, from June to December 2011, through a water level meter, in wells shown in Figure 2, in order to analyse the groundwater flow pathway at the aquifer scale.
The piezometric head was also measured once at site 2, in June 2011, through a water level meter, into a network of 22 wells (5-15 meters deep and fully screened) drilled close to Cabo San Antonio ( Figure 3) and one karst cave (cr), in order to analyse in more details the influence of discontinuous heterogeneity on groundwater flow. For this purpose, the head measurements have been carried out in a dry period (with no influence from local recharge), and in a few hours. The short timespan was chosen to minimize the influence of tide fluctuations (up to about 0.5 m at Cabo San Antonio) on the head measurements, therefore allowing the interpolation of comparable groundwater levels and reconstructing reliable contour lines. Geological map of the study area. Legend: 1. carbonate formations (Vedado, Jaimanitas, Cayo Piedra); 2. silicoclastic formation (Guane); 3. alluvial deposit; 4. biogenic or undifferentiated deposit; 5. river; 6. fault; 7. lake; 8. monitoring well; 9. seawater sampling point; 10. groundwater sampling in karst cave; 11. water sampling in river, lagoon and cienaga.
Since bacteria are abundant, ubiquitous and thrive even in extreme habitats, microbiological investigations together with hydrogeological, geochemical, isotopic, and geophysical investigations can represent a valuable tool for studying aquifer system behaviour [10]. Interesting results have been obtained in carbonate environments of southern Italy, where the potential use of microorganisms as tracers in hydrogeological studies has been assessed regarding the analysis of recharge and flow processes, and in hypersaline and hydrothermal systems as well [10,11]. These achievements have been possible thanks to the development and application of molecular biological methods in the field of hydrogeology that have permitted increasing numbers of studies pertaining to microbial communities of aquifer systems. For instance, next-generation sequencing (NGS) technology has enabled researchers to investigate the composition and function of microbial populations with unprecedented resolution and throughput.

Study Area
In the Guanahacabibes Peninsula (westernmost portion of Pinar del Río), limestones, mainly of biogenic origin, represent the prevailing outcropping rocks (Figure 1). Pliocene and Quaternary rocks and deposits overlay deformed Jurassic calcareous rocks pertaining to the structure of the Cordillera de Guaniguanico [12,13].
The westernmost portion of the Guanahacabibes Peninsula is characterised by the Vedado Formation (Pliocene-lower Pleistocene), a limestone that can reach about 200 m of thickness. The Vedado Formation passes laterally to Guane Formation (upper Pliocene-lower Pleistocene), a deltaic deposit, up to 50 m thick, characterised by conglomerates, sandstones and clay, arranged in not well-defined lenses, outcropping in the northern portion of peninsula.
Upper Pleistocene and Holocene deposits unconformably overlay the previous rocks. The Jaimanita Formation (upper Pleistocene) is made up of massive biogenic limestones of a littoral environment. It is a backreef deposit containing calcarenites, corals, coral pieces and also debris from the continental environment. In turn, oolitic calcarenites belonging to the Cocodrilo Formation (upper Pleistocene) unconformably overlay the previous limestones. The calcarenites were beach deposits, formed with different seawater depth conditions, characterised by cross-bedding typical of sand bars. The thickness is about 10 m. Finally, in the northern portion of the Guanahacabibes Peninsula, the Siguanea Formation (upper Pleistocene) is mainly made by quartzitic sandstones, representing an alluvial fan.
The general geomorphic features of the Guanahacabibes Peninsula were widely driven by the tectonic movements affecting the area since the Oligocene. An extensional tectonic regime with NW-SE strike-slip faults determined the relative lowering of the Guanahacabibes block (at south) with respect to the Guaniguanico mountain range (at north). The fault system determined the morphology of the investigated area, such as the form of the coastline in the central part (Ensenada de Cortés and Bahía de Guadiana) of the peninsula [13], where lagoons and swamps are widespread. This part seems to have functioned as a "E-W oriented graben", even if the "collapse" is rather due to a northward tilting (4 • -6 • ) of the entire peninsula. Thus, since the lower Pleistocene, the southern part of the peninsula was uplifted, owing to the tilting, with respect to the northern one that has undergone drowning phases in relationship with the relative sea-level changes.
The large amount of calcareous rocks caused the spread and dominance of dissolution processes resulting in a classic tropical karst. The terrain corresponding to the southernmost flatland of the Guanahacabibes peninsula (Vedado Formation outcrops) is a naked karst [14]. This area is characterised both by exo-and endokarst, with landforms spanning from small-scale (karren and microkarren) to medium-scale landforms; for instance, large shafts called cenotes and hoios are a kind of flooded collapsed depressions that show entrances of the cave at different levels along their vertical walls.
Instead karst terrains covered by Holocene deposits characterise the central part of the peninsula, between Ensenada de Cortés (in the east) and Bahía de Guadiana (in the west). Water wells drilled in the area reconstructed the stratigraphic sequence of the shallow deposits characterised by clay and silty clay (continental environment), marls and biodetritic limestones (littoral environment). They testify to the existence in the area of periodic shift in the environment from deposition in shallow water basins to drowning events due to relative sea-level rise.
The general climate of Cuba is humid and tropical. The average annual precipitation during the last 30 years is about 1780 mm, while the annual average temperature is 25.3 • C [4].
The natural vegetation is the tropical rainforest and along the coast the rainforests are accompanied by mangrove forests. In the region, the main economic driver is the agriculture and the main cultivated crops are citrus fruits and tobacco. Tobacco is a water-intensive and chloride-sensitive plant (maximal chloride concentration: 30 mg/L). The water is used for irrigation and domestic purposes; therefore it is of utmost importance to avoid seawater intrusion.
The study area is characterised by the presence of more or less deep-water bodies (laguna) or swamps (cienaga). The east to west system of swamps made by Cienaga de Los Remates and Cienaga de Jerusalem represents the current widest wetland of the area.
Just northward with respect to Cienaga Los Remates, lagoons of different sizes exist. They are basins with different water depths and occupy a large portion of the central part of Guanahacabibes Peninsula. These lagoons often show a rounded shape both when they are simple basins and when they are composite forms, due to the coalescence of two or more simple forms.
The shape of these surface water basins could be conditioned by the karst morphology developed in the shallow carbonate substratum (beneath Holocene deposits), as testified by the different endokarst levels along the walls of major shafts. Thus, the lagoons could be hosted in collapsed or subsidence depressions due to the connection with buried forms.
Geomorphological mapping allows the analysis of landforms in coastal areas [15][16][17]. Along the southeast and southwest of Cuba the continental shelf is shaped by several terraces represented by a staircase due to the relative changes in sea level which occurred during the Quaternary. The analysis of geomorphological features in hydrographic maps has allowed up to seven orders of marine terraces at the range of depths 8-10, 14-15, 18-20, 25, 30-32, 40 and 60-65 m below the sea level to be recognized [18]. They testify to defined sea level stands during which the development of correspondent coastal endokarsts developed where carbonate bedrock outcrops.

Hydraulic Head Measurements
The hydraulic head was measured monthly at site 1, from June to December 2011, through a water level meter, in wells shown in Figure 2, in order to analyse the groundwater flow pathway at the aquifer scale.   The piezometric head was also measured once at site 2, in June 2011, through a water level meter, into a network of 22 wells (5-15 meters deep and fully screened) drilled close to Cabo San Antonio ( Figure 3) and one karst cave (cr), in order to analyse in more details the influence of discontinuous heterogeneity on groundwater flow. For this purpose, the head measurements have been carried out in a dry period (with no influence from local recharge), and in a few hours. The short timespan was chosen to minimize the influence of tide fluctuations (up to about 0.5 m at Cabo San Antonio) on the head measurements, therefore allowing the interpolation of comparable groundwater levels and reconstructing reliable contour lines.    In well 18 ( Figure 2) a pressure transducer with a data-logger was used to measure hourly the hydraulic head fluctuations during some important rainfall events. Well 18 is about 40 meters deep, drilled within fissured and karstified limestone from 6-40 meters below ground, and fully screened in the limestone units. Four seawater samples were taken several meters from the coastline, at Los Cayuelos (slc), Bahia Guadiana (sbg), Maria La Gorda (sml) and Bahia Cortes (sbc).
Four surface water from swamps were sampled at Ciénaga Los Remates (LR1 and LR2), from a lagoon at Laguna El Pesquero (LP) and at Cuyaguateje River (rc).
Groundwater samples were collected from 19 water wells, all distributed in the Sandino municipality.
In wells 18 and 33 as well as Faro Roncali (fr) and Pacheco (p), samples were collected at different depths below the hydraulic head.
In wells 18 and 33, located close to the coast, both the temperature (T) and the electrical conductivity (EC) were measured on a monthly basis, from June to December 2011 through a borehole probe, in order to reconstruct both the thermocline and the halocline within the saturated zone, and analyse their evolution over time.

Biomolecular Investigations
Several methods have been employed to reveal microbial community composition and responses to environmental changes in many environments and in different contexts [10,11,[19][20][21][22]. An understanding of the temporal and spatial structures, functions, interactions and population dynamics of microbial communities is critical for many aspects of life, including sustainable agriculture, environmental protection and human health [23].
In the present paper a 16S rDNA polymerase chain reaction-denaturing gradient gel electrophoresis (PCR-DGGE)-based approach and next-generation sequencing (NGS) technology were used to investigate microbial communities and to assess their potential as natural (bio)tracers of mixing processes between different water bodies (groundwater, surface water, seawater), as well as of hydrochemical processes.

DNA Extraction
Within 8 h post-collection, water samples (1 L) were filtered through mixed esters of cellulose filters (S-Pak TM Membrane Filters, 47 mm diameter, 0.22 µm pore size, Millipore Corporation, Billerica, MA, USA). Immediately after filtration, the filters were stored at −80 • C until nucleic acid extraction. DNA was extracted from filters using a PowerWater ® DNA Isolation Kit (MO BIO Laboratories, Inc., Carlsbad, CA, USA).
PCR amplification was performed using the following program: 95 • C for 2 min, 20 cycles of denaturation at 92 • C for 45 s, annealing at 50 • C for 2 min, and extension at 72 • C for 1 min and 45 s, and a single final extension at 72 • C for 5 min.
Touchdown PCR was performed with an initial denaturation step of 94 • C for 5 min, 10 touchdown cycles of 94 • C for 1 min, 65 • C (−1 • C per cycle) for 1 min and 72 • C for 3 min, 5 cycles of 94 • C for 1 min, 55 • C for 1 min and 72 • C for 2 min, and a final elongation step of 72 • C for 4 min.
The presence of PCR products was confirmed by analysing 5 µL of product on 1.5% agarose gels and staining with ethidium bromide.
A cluster analysis of the DGGE patterns was performed using FPQuest Software Version 5.1 (BIO-RAD Laboratories, Inc., Hercules, CA, USA). The similarity in the profiles was calculated based on the Dice coefficient with the UPGMA (Unweighted Pair Group Method using Arithmetic averages) clustering algorithm.

Dual Index 16S rRNA Gene Amplicon Library Preparation and Bioinformatics Analysis
Next-generation sequencing (NGS) protocol was performed at BMR Genomics srl (Padova, Italy). For each sample collected at different depths from wells 18 and 33, the V3-V4 regions of the 16S rRNA gene were amplified using the primers 331F and 797R [25].
The classification step used ClassifyReads, a high-performance implementation of the Ribosomal Database Project (RDP) Classifier described in [26]. This process involved matching short subsequences of the reads (called words) to a set of 16S rRNA gene reference sequences (the taxonomic database used was an Illumina-curated version of the May 2013 release of the Greengenes Consortium Database). The accumulated word matches for each read were used to assign reads to a particular taxonomic classification.
To compare the microbiota among the samples, a Principal Component Analysis (PCA) was carried out by R software package. For this purpose, the function PCA from the FactoMineR package was exploited to examine the relative abundance values of operational taxonomic units (OTUs) at the genus level. The first two principal components were used to generate a scatter plot.

Groundwater Numerical Model
A numerical model of the aquifer, based on a simple conceptualization of the study case, has been developed in order to explain the electrical conductivity profile observed at the monitoring well 18. In particular, a numerical simulation of the seawater intrusion in a dual permeability coastal karst aquifer with two conduit networks has been performed using a three-dimensional variable-density groundwater flow and multi-species transport SEAWAT model. SEAWAT is a coupled version of MODFLOW [27] and MT3DMS [28][29][30].
The objective of the numerical model is to describe the variation of the salinity concentration observed in well 18 as a function of the depth and it is based on the parameter values of a porous medium aquifer with two karst conduits with large values of both the hydraulic conductivity and the effective porosity.
We have chosen the numerical grid to describe two conduits with the minimum possible number of cells. A refinement in the y and z directions was performed around both conduits to delineate them and the size of the cells around the conduits are 1 m × 1 m. We did not refine cell sizes along the x direction since we are interested to the region near the sea where the salt intrusion wedge takes place. See Figure 4 for details. The hydraulic conductivity of the fractured medium is assigned to be 1045 m/day, based on a mean velocity of 10 m/day for a fractured medium [31]. For the conduit system a very high value of the hydraulic conductivity of 2.4 × 10 6 m/day was assigned based on the value of 2371 m/day as a mean velocity for a typical karst conduit (usually in the range of km/day). As a matter of fact, an average velocity of 1728 m/day was statistically determined by regression analyses of more than 3000 tracing tests worldwide [32].
On the whole, the hydraulic conductivities used in the present work are not so far from those utilized in recent papers focused on the simulation of groundwater flow and saltwater wedge shape in karst environments (e.g., up to 4572 m/day for fractured media and up to 9.1 × 10 5 m/day for conduits in [33]). The effective porosity is set to 0.1 in the fractured medium and 1.0 along the conduit. The longitudinal dispersivity is assigned 10 m in the porous medium and very small (0.3 m) in the conduit since advection is dominating compared to dispersion (in agreement with previous simulations in karst aquifers [33]).
The finite difference grid of the three-dimensional model covering a total area of 10 km 2 and 45 m depth is set to a rectangular parallelepiped of 10 km long (with 400 column), 1 km width (with 33 rows) and a thickness of 26 layers in the cross section, for a total number of 343.200 cells (Figure 4). The cell sizes were chosen as much as possible while keeping a relatively low computational cost and time.
The horizontal discretization for each cell is set uniformly to 2 m except for columns 8, 9, 10, 11 and 16, 17 where a fine resolution of 1 m each is implemented in correspondence with both conduits. An upper conduit is situated at layer 10 (−14.5 meters above sea level (m.a.s.l.)) and row 16 where refinement is also implemented. A second conduit is located below to the first one and parallel to it at layer 16 (−24.5 m.a.s.l.) and row 16. Both conduits span the entire horizontal area of 10 km long. For simplicity, the size of the horizontal conduits is assumed to be constant all over the parallelepiped. The outlet of both karst conduits, of an area of 1 m 2 , are in contact with the sea boundary on the left side where the saline wedge would exist.
Assigned hydraulic heads of 0 m.a.s.l. and 10 m.a.s.l. are set as boundary conditions, respectively at the sea boundary (x = 0 m) and at the inland boundary (x = 10,000 m). A constant salt concentration is fixed to 37 kg/m 3 at the sea boundary. The concentration at the inland boundary condition is 0.0 kg/m 3 as freshwater. The values of the density fluid in the SEAWAT VDF packages goes from 1000 kg/m 3 corresponding to freshwater to 1026.0 kg/m 3 which corresponds to a typical seawater density.

Groundwater Flow Pathway and Regime
The groundwater flow pathway was reconstructed only in the eastern part of the study area (site 1 in Figure 1), where the monitoring points were sufficient to draw a piezometric map. It depicts radial flow towards the surrounding coastlines ( Figure 2). The hydraulic head intercepts the ground surface close to some topographic depressions, therefore causing the groundwater to appear at the surface. This phenomenon is responsible for the so-called cienagas, not so far from the coast, whose existence and thickness strictly depends on groundwater surface fluctuations. The higher head and radial flow, observed close to the lagoons, suggest that they can feed the groundwater, despite the low-permeability media, which crop out within the same area ( Figure 1). The lagoons are fed by the Cuyaguateje River. The ingression of waters from El Pesquero lagoon to the surrounding groundwater is already described in [9] using isotopic and chemical investigations, and further supported by biomolecular ones (see hereafter).
At local scale (site 2 in Figures 1 and 3), the groundwater flow pathway is characterised by several "drainage axes", in agreement with the high development of karst conduits in the area. However, despite the heterogeneity, the hydraulic heads measured in the observation wells are all compatible with one groundwater body, pervasively flowing within the network of fissures and karst openings. Moreover, taking into consideration the faults' patterns (Figure 3), there is no relationship between drainage axes and fault zones, therefore suggesting that these faults do not act as conduits that enhance flow, differently to findings in other karst environments (e.g., [35]). On the contrary, taking into account the increase of the hydraulic gradient from 0.02% (far from the fault zones) to 1.9-4.2% (close to the fault zones) along the south-western edge of the test site, it seems that a lower-permeability fault core is associated to the local fault zones, which partially impedes fluid flow, causing the gradient to locally rise. This relationship between hydraulic gradient increase and lower-permeability fault core is in agreement with the results of experimental studies carried out in other carbonate aquifers, fissured and karstified, where the groundwater flow lines are approximately orthogonal to the direction of the fault core [36][37][38][39][40][41][42][43].
Concerning the groundwater regime, the precipitations result in rapid effective infiltration, and fast hydraulic head increase ( Figure 5). At the same time, rapid depletion is observed just after the end of rainfall, therefore showing an irregular regime. The rapid head rise due to huge effective infiltration also causes the mixing zone between fresh-and seawater to move downwards, as expected (see abrupt variation in EC at the well bottom during recharge events close to well 18, in Figure 5).

16S rDNA PCR-DGGE Analysis of Sea-, Surface-, and Groundwater Microbial Communities
The main findings obtained from the cluster analysis of 16S rDNA PCR-DGGE profiles can be summarized as follows ( Figure 6). Firstly, significant differences were found between bacterial communities in sea-, surface-and groundwater samples collected at the study site. Moreover, important differences can be found also between samples of the same water type (from the hydrochemical point of view), suggesting the existence of different factors, site-or zone-specific, influencing bacterial communities. For example ( Figure 6, box A), the upper cluster groups almost all the wells located in the central part of the study area as well as the water sample collected in the El Pesquero lagoon (LP) and seawater (sbc); this latter point is in agreement with the groundwater flow net, as shown in Figure 2. This clusterisation further supports the hypothesis of a hydraulic interconnection between the lagoon system and the surrounding groundwater, as suggested in function of the hydraulic head distribution and the previous isotopic and hydrochemical investigations [9].
Concerning the vertical zoning along the saturated zone, close to the coast, well 18 is characterised by three sub-clusters within the same main cluster ( Figure 6, box B). The sub-clusterisation is strictly linked to the sampling depth (4-31 meters below the hydraulic head (mbhh)). This result is in agreement with a progressive influence of mixing between ground-and seawater on microbial community composition. Conversely, the communities detected at each sampling depth in well 33 (10-45 mbhh) do not cluster all together. This observation is in agreement with the more complex variation in chemical features observed within the saturated zone investigated in well 33. A detailed geochemical characterisation of these waters suggested [9] that the samples from well 33 showed Ca enrichment and a Mg, Na and K depletion compared with the curves representing binary mixing between seawater and freshwater or rainwater, on which the samples from well 18 are clustered. Moreover, the calcium excess and sulfate deficit observed in well 33 samples have been explained by additional calcium carbonate dissolution, probably occurring in the presence of bacterial reduction in organic-rich sediments. In fact, the sulfate decline observed in well 33 occurred at a chloride content that cannot be explained by dilution.

16S rDNA PCR-DGGE Analysis of Sea-, Surface-, and Groundwater Microbial Communities
The main findings obtained from the cluster analysis of 16S rDNA PCR-DGGE profiles can be summarized as follows ( Figure 6). Firstly, significant differences were found between bacterial communities in sea-, surface-and groundwater samples collected at the study site. Moreover, important differences can be found also between samples of the same water type (from the hydrochemical point of view), suggesting the existence of different factors, site-or zone-specific, influencing bacterial communities. For example (Figure 6, box A), the upper cluster groups almost all the wells located in the central part of the study area as well as the water sample collected in the  In order to support the hypotheses made for well 33 on the basis of geochemical analyses, the 16S rRNA gene sequencing was performed using Illumina technology. The microbial community was then characterised at different depths to analyse in more detail the reasons why the same progressive variation detected in well 18 was not observed.

Next-Generation Sequencing Results
MiSeq runs produced a total of 2,496,007 raw data (reads), including V3 and V4 regions of the 16S rRNA gene. After quality filtering of raw data, 2,392,542 sequences were used for subsequent analysis. Classified OTUs belonged to 30 phyla, 64 classes, 127 orders, 277 families, 811 genera and 2708 species. The top five dominant phyla Proteobacteria, Firmicutes, Actinobacteria, Bacteroidetes and Nitrospirae accounted, on average, for 97.86% and 93.27% of sequences retrieved from wells 18 and 33, respectively ( Figure 7). Overall, Proteobacteria represented the main phylum in all the samples with relative abundance values ranging from 73.13-81.66% in well 18 and from 40.11-71.61% in well 33. Within this phylum, Alphaproteobacteria are numerically dominant in marine ecosystems [44], but they are also, as well as Betaproteobacteria, copious and ubiquitous in freshwater [45][46][47][48]. On the other hand, Gammaproteobacteria are more abundant in oceans [49] or saline lakes [50].

Next-Generation Sequencing Results
MiSeq runs produced a total of 2,496,007 raw data (reads), including V3 and V4 regions of the 16S rRNA gene. After quality filtering of raw data, 2,392,542 sequences were used for subsequent analysis. Classified OTUs belonged to 30 phyla, 64 classes, 127 orders, 277 families, 811 genera and 2708 species. The top five dominant phyla Proteobacteria, Firmicutes, Actinobacteria, Bacteroidetes and Nitrospirae accounted, on average, for 97.86% and 93.27% of sequences retrieved from wells 18 and 33, respectively (Figure 7). Overall, Proteobacteria represented the main phylum in all the samples with relative abundance values ranging from 73.13%-81.66% in well 18 and from 40.11%-71.61% in well 33. Within this phylum, Alphaproteobacteria are numerically dominant in marine ecosystems [44], but they are also, as well as Betaproteobacteria, copious and ubiquitous in freshwater [45][46][47][48]. On the other hand, Gammaproteobacteria are more abundant in oceans [49] or saline lakes [50].

Microbial Communities in Well 18
In well 18, Alphaproteobacteria ranged from 5.23%-33.77%. An increasing trend of relative abundance with depth was observed for Gammaproteobacteria (from 8.87% at 4 m below hydraulic head (bhh) to 67.49% at 28 mbhh) as opposed to Betaproteobacteria which reached their lowest values at 28 and 31 mbhh (Figure 8). Actinobacteria, a phylum comprising of microorganisms widely distributed in both terrestrial and aquatic (including marine) ecosystems [51], were retrieved at percentages ranging from 3.17%-6.99% whereas Bacteroidetes and Firmicutes represented 7.52%-16.75% and 3.25%-13% of sequence reads, respectively. PCA of the NGS data at genus level clearly showed that in well 18 a great diversification among microbial communities occurred with depth ( Figure 9A). At 28 and 31 mbhh, microbial communities seemed to be more shaped by the presence of the genera Marinobacter, Idiomarina and Pseudoidiomarina which include the species Marinobacter sediminum, Marinobacter santoriniensis, Idiomarina zobellii and Idiomarina fontislapidosi, etc. isolated from saline habitats with a wide range of salinities, such as coastal and oceanic waters, coastal sediments, and inland hypersaline wetlands. On the other hand, communities retrieved at 4 mbhh were more

Microbial Communities in Well 18
In well 18, Alphaproteobacteria ranged from 5.23-33.77%. An increasing trend of relative abundance with depth was observed for Gammaproteobacteria (from 8.87% at 4 m below hydraulic head (bhh) to 67.49% at 28 mbhh) as opposed to Betaproteobacteria which reached their lowest values at 28 and 31 mbhh (Figure 8). Actinobacteria, a phylum comprising of microorganisms widely distributed in both terrestrial and aquatic (including marine) ecosystems [51], were retrieved at percentages ranging from 3.17-6.99% whereas Bacteroidetes and Firmicutes represented 7.52-16.75% and 3.25-13% of sequence reads, respectively. PCA of the NGS data at genus level clearly showed that in well 18 a great diversification among microbial communities occurred with depth ( Figure 9A). At 28 and 31 mbhh, microbial communities seemed to be more shaped by the presence of the genera Marinobacter, Idiomarina and Pseudoidiomarina which include the species Marinobacter sediminum, Marinobacter santoriniensis, Idiomarina zobellii and Idiomarina fontislapidosi, etc. isolated from saline habitats with a wide range of salinities, such as coastal and oceanic waters, coastal sediments, and inland hypersaline wetlands. On the other hand, communities retrieved at 4 mbhh were more influenced by the genera Rhodobacter, Hydrogenophaga, Niastella, Novosphingobium, Giesbergeria, Methylotenera, Ramlibacter and Sphingomonas, comprising species isolated from freshwaters, soils and also marine surface waters. Overall, the distribution of samples along the water column reflected the increasing trend of EC in the well, with microbial communities in the upper layer more linked to freshwaters and soil ecosystems than those of the deeper layers which were more influenced by genera typical of marine and saline environments. A transition zone could be represented by communities retrieved from 16-27 mbhh that grouped together and were placed between water samples representing the extremities of the salt gradient. influenced by the genera Rhodobacter, Hydrogenophaga, Niastella, Novosphingobium, Giesbergeria, Methylotenera, Ramlibacter and Sphingomonas, comprising species isolated from freshwaters, soils and also marine surface waters. Overall, the distribution of samples along the water column reflected the increasing trend of EC in the well, with microbial communities in the upper layer more linked to freshwaters and soil ecosystems than those of the deeper layers which were more influenced by genera typical of marine and saline environments. A transition zone could be represented by communities retrieved from 16-27 mbhh that grouped together and were placed between water samples representing the extremities of the salt gradient.

Microbial Communities in Well 33
In well 33, Alphaproteobacteria represented 3.64%-9.77% of reads. Betaproteobacteria were the dominant class of microbial communities up to 25 mbhh (>50%) whereas at the greatest depths their abundance was reduced by more than half. Gammaproteobacteria were retrieved with values ranging from 8.66%-21.88% (Figure 8). Deltaproteobacteria, which were detected also in well 18, reached the highest values of relative abundance at 10 , 25 and 45 mbhh (5.93%, 4.39% and 5.97%, respectively; Figure 8). Within this class, several genera such as Desulfovibrio, Desulfonatronum, Desulfomonile, Desulfosarcina, Desulfococcus, Desulfotignum etc., known to include sulfate reducer bacteria, were found in water samples.
In addition, with reference to sulfate reduction, in well 33 the phylum Nitrospirae was retrieved at relatively high abundances at 10 and 25 mbhh (7.55% and 5.52%), being largely represented by the genus Thermodesulfovibrio (species Thermodesulfovibrio aggregans). This genus was established in 1994 following the isolation of a novel thermophile, Thermodesulfovibrio yellowstonii, from a hydrothermal vent in Yellowstone Lake [52]. Other species, such as Thermodesulfovibrio islandicus [53], Thermodesulfovibrio hydrogeniphilus [54], Thermodesulfovibrio aggregans and Thermodesulfovibrio thiophilus [55], were subsequently isolated from terrestrial hot springs in Iceland and Tunisia and thermophilic methanogenic sludge samples. However, molecular signatures of this group have been detected in deep aquifer systems of the Great Artesian Basin, Australia [56], in formation waters and oil from several reservoirs in China [57][58][59] and in deep granitic groundwater from the Grimsel Test Site in Switzerland [60]. All known Thermodesulfovibrio isolates are strict anaerobes and use a limited number of organic electron donors, which are incompletely oxidized to acetate. They all are capable of dissimilatory reduction of sulfate and thiosulfate, and some species can also use sulfite, Fe(III)-NTA, and nitrate as electron acceptors [61].
Finally, Actinobacteria were retrieved at percentages ranging from 7% to 10.32%, Bacteroidetes ranged between 0.90% and 4.27% and Firmicutes between 7.56% and 38.76%, with the highest relative abundance values recorded at depths over 25 mbhh.
PCA indicated two different clusters ( Figure 9B). The first cluster comprised microbial communities detected at 10 and 25 mbhh, mainly influenced by members of the genera Methyloversatilis, Thauera, Sulfuritalea, Thermodesulfovibrio, Propionivibrio and Methylocaldum isolated from different types of habitats such as soils, freshwaters and marine waters. The impact on microbial community structure of the two genera involved in the sulfur cycle Sulfuritalea (chemolithoautotrophic growth under anoxic conditions by the oxidation of reduced sulfur compounds and hydrogen) and Thermodesulfovibrio (thermophilic bacteria that reduce sulfate and other sulfur compounds [52,62]) explains the sulfate deficit observed in this well. Microbial communities retrieved beyond 30 mbhh were grouped in the second cluster. In this case, the genera

Microbial Communities in Well 33
In well 33, Alphaproteobacteria represented 3.64-9.77% of reads. Betaproteobacteria were the dominant class of microbial communities up to 25 mbhh (>50%) whereas at the greatest depths their abundance was reduced by more than half. Gammaproteobacteria were retrieved with values ranging from 8.66-21.88% (Figure 8). Deltaproteobacteria, which were detected also in well 18, reached the highest values of relative abundance at 10, 25 and 45 mbhh (5.93%, 4.39% and 5.97%, respectively; Figure 8). Within this class, several genera such as Desulfovibrio, Desulfonatronum, Desulfomonile, Desulfosarcina, Desulfococcus, Desulfotignum etc., known to include sulfate reducer bacteria, were found in water samples.
In addition, with reference to sulfate reduction, in well 33 the phylum Nitrospirae was retrieved at relatively high abundances at 10 and 25 mbhh (7.55% and 5.52%), being largely represented by the genus Thermodesulfovibrio (species Thermodesulfovibrio aggregans). This genus was established in 1994 following the isolation of a novel thermophile, Thermodesulfovibrio yellowstonii, from a hydrothermal vent in Yellowstone Lake [52]. Other species, such as Thermodesulfovibrio islandicus [53], Thermodesulfovibrio hydrogeniphilus [54], Thermodesulfovibrio aggregans and Thermodesulfovibrio thiophilus [55], were subsequently isolated from terrestrial hot springs in Iceland and Tunisia and thermophilic methanogenic sludge samples. However, molecular signatures of this group have been detected in deep aquifer systems of the Great Artesian Basin, Australia [56], in formation waters and oil from several reservoirs in China [57][58][59] and in deep granitic groundwater from the Grimsel Test Site in Switzerland [60]. All known Thermodesulfovibrio isolates are strict anaerobes and use a limited number of organic electron donors, which are incompletely oxidized to acetate. They all are capable of dissimilatory reduction of sulfate and thiosulfate, and some species can also use sulfite, Fe(III)-NTA, and nitrate as electron acceptors [61].
Finally, Actinobacteria were retrieved at percentages ranging from 7% to 10.32%, Bacteroidetes ranged between 0.90% and 4.27% and Firmicutes between 7.56% and 38.76%, with the highest relative abundance values recorded at depths over 25 mbhh.
PCA indicated two different clusters ( Figure 9B). The first cluster comprised microbial communities detected at 10 and 25 mbhh, mainly influenced by members of the genera Methyloversatilis, Thauera, Sulfuritalea, Thermodesulfovibrio, Propionivibrio and Methylocaldum isolated from different types of habitats such as soils, freshwaters and marine waters. The impact on microbial community structure of the two genera involved in the sulfur cycle Sulfuritalea (chemolithoautotrophic growth under anoxic conditions by the oxidation of reduced sulfur compounds and hydrogen) and Thermodesulfovibrio (thermophilic bacteria that reduce sulfate and other sulfur compounds [52,62]) explains the sulfate deficit observed in this well. Microbial communities retrieved beyond 30 mbhh were grouped in the second cluster. In this case, the genera Acinetobacter and Lactobacillus seemed to have a major role in the sample "spatial distribution". Acinetobacters are mostly free-living saprophytes while Lactobacillus is the largest genus within the group of lactic acid bacteria (LAB). A very large number of species belonging to these genera can be found ubiquitously in nature, for example in water and soil, on vegetables and human skin, and in wastewater and human faeces. Discovery at these depths could indicate contaminant ingress into the aquifer system. Although further investigations are necessary to define exactly the sources and transport modalities of these contaminants, we hypothesize that their presence could be determined by the intrusion of sewage-polluted seawater rather than infiltration through the soil since well 33 is located within an area with a very low anthropogenic impact.

Saltwater Wedge Shape
In wells 18 and 33, the EC increases with increasing depth, showing mixing processes between ground-and seawater. In both the observation points, it is possible to detect a transition layer (TL). This layer is wider into the well 33 (from around −20 to −40 m.a.s.l.) and thinner into the well 18 (from around −25 to −31 m.a.s.l.; Figure 10A,B. Moreover, during the entire monitoring period, in well 18 it is possible to detect a step-like shape, with the steps always located at the same altitude (around −14 and −23 m.a.s.l.). The groundwater temperature also varies with depth, with shallower values strictly related to the seasonal temperature ( Figure 10C,D).

Numerical Results
A simple numerical model has been implemented in order to analyse in detail the role of karst conduits on the halocline and then on the saltwater wedge shape.
A series of preliminary transient simulations were performed in order to identify the best set of hydraulic parameters that reproduce the EC profile of Figure 10 at a local scale. Figure 11 shows the salinity concentration as a function of the depth (black line) simulated through the numerical model; the parallel grey lines show the position of the two karst conduits.
As can be seen from Figure 11 and Figure 12, a step-like shape may be identified in the salt According to previous studies in fissured media (e.g., [63][64][65][66]), the location of the sharp variations in EC with depth identifies the location of major and more transmissive openings in a lower permeability medium.
Merging this information with the analysis of the local bathymetry, good agreement can be observed between two terraces (−14, −15, and −25 m.a.s.l.) identified by [31] and the steps in EC detected in well 18. This observation suggests that during one of the marine transgressions, a karst system (previously developed as a function of a different base-level) was submerged within the current saturated zone. Therefore, a karstified layer is actually interbedded within fissured limestone, and then karst conduits are the higher permeability openings causing the observed step-like shape of EC profiles. Therefore, different interactions between ground-and seawater are expected within the saturated zone.

Numerical Results
A simple numerical model has been implemented in order to analyse in detail the role of karst conduits on the halocline and then on the saltwater wedge shape.
A series of preliminary transient simulations were performed in order to identify the best set of hydraulic parameters that reproduce the EC profile of Figure 10 at a local scale. Figure 11 shows the salinity concentration as a function of the depth (black line) simulated through the numerical model; the parallel grey lines show the position of the two karst conduits.  As can be seen from Figures 11 and 12, a step-like shape may be identified in the salt concentration. One at the sea boundary (x = 0) and the other two in concomitance with both conduits similar to those appearing in Figure 10.   Apparently, the presence of the two karst conduits prevent the formation of the typical salt wedge intrusion and, instead, show a step-like shape. In this case, the freshwater coming from the aquifer conduit pushes back the seawater even if it may enter into the karst conduits depending on the choice of the different parameters (e.g., the hydraulic conductivity or the initial value of the heads).
Electrical conductivity (as an indirect measure of the salinity) observed at well 18 is superimposed with a different x-axis scale to the modelled concentration in Figure 11. It is clear that the model is able to reproduce the two-step progress of the observed electrical conductivity. These results support the hypothesis (two conduits at about −15 and −25 m.a.s.l.) taken during the development of the conceptual model.

Discussion and Conclusions
The efficacy of an experimental modelling approach was tested in order to better characterise both the hydraulic features and behaviour of a coastal carbonate aquifer and acquire useful information to prevent groundwater salinization ( Figure 13). The prevention of seawater intrusion at the study area is of utmost importance because the main economic driver is the agriculture and the main cultivated crops are chloride-sensitive plants.
The hydrogeological behaviour at the basin scale is influenced mainly by the location of the lagoons (fed by the Cuyaguateje River), that govern the hydraulic head in the central part of the study site, causing the groundwater to flow radially towards the surrounding coast. This phenomenon has been further supported by the results of biomolecular investigations. The phreatic surface is locally intercepted by morphological depressions, and several wetlands (cienagas) are generated parallel to the southern and the western coastlines.
At the local scale, both layered and discontinuous heterogeneities play an important role in governing the groundwater flow and significantly influence the shape of the saltwater wedge.
Concerning the discontinuous heterogeneity, the detailed reconstruction of the groundwater flow net close to Cabo de San Antonio has shown that lower permeability fault cores probably influence the hydraulic gradient, causing a significant increase (two orders of magnitude) at the local scale. From an experimental point of view, in the present study more attention was given to (bio)tracers investigated through biomolecular analyses. As a matter of fact, microbial communities were able to refine knowledge and/or further support hypotheses based upon hydrogeological and hydrochemical investigations, with emphasis on (i) interconnections between surface-(lagoons and cienagas) and groundwater; (ii) mixing processes between ground-and seawater, including the differences in the halocline features in different areas (see examples of wells 18 and 33); and (iii) hydrochemical processes within the saturated zone close to sea (see example of well 33).
The conceptual model based upon the interdisciplinary investigations was then used to implement a schematic three-dimensional numerical model, with emphasis on the simulation of both the extent and the shape of seawater intrusion in the dual permeability (fissured limestone and karst conduits) of the coastal carbonate aquifer. The numerical model further demonstrated that karst conduits play an important role in influencing the distribution of water salinity within the saturated zone, according to the results of other recent studies coupling density-dependent flow and transport models [33]. In the case study, the simple aquifer numerical model can reproduce the shape of the EC profile observed in monitoring well 18 and consequently strengthens the hypotheses for the conceptual model.  Regarding the layered heterogeneity, the existence of karstified layers in fissured limestone causes the halocline to be step-like within the mixing zone between fresh groundwater and seawater. According to the literature and the site-specific results of the numerical model implemented within the present study (see hereafter), the location of the sharp variations in EC with depth identifies the location of major and relatively higher transmissive openings in a lower permeability medium. At the study site, the location of major karst conduits was also supported by the analysis of the local bathymetry and then the relationships between karst evolution and marine transgressions.
From the methodological point of view, this study shows that the knowledge about mixing processes in this type of heterogeneous media can be further improved using multilevel sampling systems, to avoid possible misunderstandings due to intra-well mixing in fully screened wells. At the same time, based upon the results of this work and a previous study [9], the interdisciplinary characterization of the Guanahacabibes Peninsula carbonate aquifer further demonstrates the great effectiveness of approaches that merge hydrogeological, hydrochemical, isotopic and biomolecular investigations to understand the behaviour of complex systems.
From an experimental point of view, in the present study more attention was given to (bio)tracers investigated through biomolecular analyses. As a matter of fact, microbial communities were able to refine knowledge and/or further support hypotheses based upon hydrogeological and hydrochemical investigations, with emphasis on (i) interconnections between surface-(lagoons and cienagas) and groundwater; (ii) mixing processes between ground-and seawater, including the differences in the halocline features in different areas (see examples of wells 18 and 33); and (iii) hydrochemical processes within the saturated zone close to sea (see example of well 33).
The conceptual model based upon the interdisciplinary investigations was then used to implement a schematic three-dimensional numerical model, with emphasis on the simulation of both the extent and the shape of seawater intrusion in the dual permeability (fissured limestone and karst conduits) of the coastal carbonate aquifer. The numerical model further demonstrated that karst conduits play an important role in influencing the distribution of water salinity within the saturated zone, according to the results of other recent studies coupling density-dependent flow and transport models [33]. In the case study, the simple aquifer numerical model can reproduce the shape of the EC profile observed in monitoring well 18 and consequently strengthens the hypotheses for the conceptual model.