Groundwater Temperature as an Indicator of the Vulnerability of Karst Coastal Aquifers

Abstract

Coastal karst aquifers show a three-dimensional vulnerability, which consists of the whole of the “intrinsic vulnerability” and the “groundwater vulnerability to seawater intrusion”. The results of a study carried out in the Salento karst coastal aquifer (southern Italy) show that temperature, as well as being a reliable tracer of groundwater flow, is also an effective indicator of vulnerability in anisotropic media. The trend of isotherms related to a cross-section of the aquifer thermal field, combined with geological, geomorphological, and hydrogeological information, allows the role of faults and dolines in the mass transport from ground surface to be inferred. Isotherm trends may also give information on the permeability distribution along faults. A specific temperature value evidence the saltwater top, thus indicating the groundwater vulnerability to salinization.

1. Introduction

The European Water Framework Directive (WFD, Directive 2000/60/EC of the European Parliament and of The Council of 23 October 2000) institutes a framework for Community action in the field of water policy. The implementation of WFD, besides targeting toward a general good qualitative and quantitative status of groundwater, also entails the delineation of protection zones for public drinking water supplies, which include the recognition of the intrinsic vulnerability of aquifers.

The notion of intrinsic vulnerability of aquifers was first introduced by Margat in 1968 [1]; since then several definitions have appeared. Here we will use the definition proposed by Dörfliger et al. (1999) [2] (p. 2): “Intrinsic vulnerability represents the inherent hydrogeological and geological characteristics, which determine the sensitivity of groundwater to contamination by human activities”. The intrinsic vulnerability considers the geological, hydrological, and hydrogeological properties of an area that influence the pollutants’ transport from the ground surface to the saturated zone of the uppermost aquifer. The intrinsic vulnerability does not consider the nature of the contaminants and the contamination scenario. Literature on intrinsic vulnerability mapping refers to many decades of research: a useful review of methods and applications to a variety of hydrogeological settings is in Witkowski et al. [3]. The karstic aquifers are a special case because of their unique features, the Report of COST (Scientific and Technological Cooperation of the European Commission) Action 620 [4] is a thorough treatise on methods for the intrinsic vulnerability assessment of karst aquifers. Notwithstanding such efforts, some studies show that a different set of indicators and/or different mapping protocols of a same karst aquifer may lead to questionable results [5].

To ascertain the vulnerability of water resources in coastal aquifers is much more difficult compared to continental aquifers. Coastal karst aquifers and coastal aquifers in general show, in fact, a three-dimensional vulnerability, which consists of the whole of the “intrinsic vulnerability” and the “groundwater vulnerability to seawater intrusion”. Seawater intrusion is induced by a hydraulic disequilibrium between groundwater and seawater and causes groundwater salinization. Useful general reviews about seawater intrusion (functioning, monitoring, modelling, and remediation practices) are in Custodio and Bruggeman (1987) [6], FAO (1997) [7], Bear et al. (1999) [8], Cheng and Ouazar (2003) [9] and Tulipano et al. (2005) [10]. Notwithstanding the wide knowledge and early expert warnings on the consequences of increasing anthropic pressures and climate change on the qualitative and quantitative status of groundwater in coastal aquifers [11], most coastal aquifers today show serious salinization problems.

To date there is not any specific method to define the vulnerability in karst coastal aquifers. Literature mostly considers methods for coastal porous aquifers (or porous aquifer systems) of limited thickness and extension [12,13,14,15,16,17,18,19,20,21,22,23,24,25] that mainly represent variations of the original GALDIT method [26,27]. GALDIT is the acronym for Groundwater occurrence, Aquifer hydraulic conductivity, water Level above mean sea level, Distance from the coastline, Impact status of existing seawater intrusion and aquifer Thickness, selected as vulnerability indicators. One only method deals with a coastal karst aquifer [28] considered as porous medium, but it does not prove suitable for karst aquifers. Discontinuities in hard rock coastal aquifers are, in fact, difficult to synthesize in a vulnerability indicator [29].

To date, the topic of the three-dimensional vulnerability of coastal aquifers, and even more so of coastal karst aquifers, lacks a systematization. With the aim of addressing this issue, we describe some basic indicators of the three-dimensional vulnerability of coastal aquifers in general and discuss the limits of major indicators for their application to karst and coastal karst aquifers. Then, we consider a domain (the coastal karst aquifer of Salento, southern Italy) as a demonstration site for applying a specific approach to vulnerability issues. The Salento aquifer is representative of permeability features, influenced by both karst processes and tectonics, widely found in other karst coastal aquifers. The outline of groundwater protection areas and vulnerability for the Salento aquifer, in fact, encounters problems in clarifying the role of both karst forms and main discontinuities in the vulnerability issues. The type of groundwater monitoring, based on multi-parametric logs carried out in a monitoring net projected for the control of seawater intrusion, allows setting a method for solving such questions. The method is based on interpolation of temperature logs to reconstruct the thermal field and on the interpretation of its trends. The object of the study is to prove the potential of temperature in tracing the processes as mass transport and seawater intrusion and in the recognition of the role of faults and karst forms in the context of the vulnerability of the anisotropic medium.

2. The 3D Vulnerability of Coastal Aquifers

Figure 1 outlines a homogeneous and isotropic unconfined coastal aquifer of an island and the general framework of its 3D vulnerability. To define the path of a pollutant from an external source towards groundwater, we consider the origin of pollution, the target, and the pathway to reach the target from the origin.

For intrinsic vulnerability as such (Figure 1 and Figure 2), the origin is the ground surface, the target is the water table, and the pathway is downward-oriented through the unsaturated zone. To fully define the three-dimensional vulnerability of coastal aquifers we consider additional origins, targets, and paths. Saline sources are alongside the sea bed (the ocean-land interface) and the freshwater–saltwater interface (origin), and their target is the lower surface of groundwater. The pathways have components from lateral- to upward-oriented along this lower surface. Thus, considering the island of Figure 1, the target is the whole (upper and lower) surface of the 3D groundwater lens floating on saltwater.

2.1. Sources of Pollution/Salts and Vulnerability Indicators

A whole discussion about indicators of coastal aquifer vulnerability is out of the scope of this paper. Hereafter, only some basic ideas are reported about all coastal aquifers. Then the indicators are reviewed in the light of the specific features of karst formations. We follow the scheme of the “Type of pollution sources/processes” of Figure 2 that shows also the linked indicators.

For case (a) (Figure 2), the indicators refer to those factors that regulate the contaminant transport from ground surface, namely the topographic surface, the unsaturated zone, and the precipitation regime. Low effective infiltration, diffused recharge, large thickness of covers, low permeability of covers and absence of discontinuities should indicate low intrinsic vulnerability, while high effective infiltration, low thickness of covers, and the presence of discontinuities should be used as indicators of high intrinsic vulnerability. In coastal aquifers under arid climates, salt spray and flood irrigation practices should be considered important pollution sources, able to cause the long-term increase of the salinity of groundwater recharge [30,31].

For case (b) (Figure 2) the source is seawater. Lateral salt pollution may also derive from salt lagoons [32], basins for aquatic farming, coastal streams and canals, and swamps fed by seawater because of tidal movements and sea level long-term variation. The indicators are the groundwater hydraulic head, the permeability of the aquifer, the presence of covers along the coast, and the depth of the geological limit. Thus, a low or high vulnerability is expected for high or low hydraulic heads and low or high horizontal hydraulic conductivity (K_h) along the coast, respectively. As for lagoons, saline rivers and channels, their interaction with groundwater depends on time variations of their hydraulic features: groundwater normally represents an inflow to surface waters, but tides, storm surges and groundwater exploitation, or specific local conditions, can trigger head reversal with inflow of surface waters into groundwater.

However, more salt sources and processes can produce groundwater salt increase (Figure 1 and Figure 2). The available data show that most of salinization in coastal aquifers comes from salt waters laying beneath groundwater and not from seawater, normally influencing groundwater only in the coastal zone. The target is the bottom of fresh groundwater and the source is the underlying saltwater; the pathway has a vertical component with an upward sense (case (c), Figure 2). Salt waters are slow-moving saline fluids corresponding to seawater that entered the aquifers during earlier transgressions and resided for a period long enough for the fluid-composition to be modified through water–rock interaction processes [30,33,34,35,36,37,38]. Their chemical features resemble those of saline waters found in large sedimentary basins [33] and in some cases, their salt content is higher than present seawater. Normally, these salt waters come to light when mobilized by exploitation [34]: the process is termed saltwater intrusion or saltwater upconing [35]. In the case (c) the indicators are the hydraulic head, the hydraulic conductivity, and the depth of the geological limit. Thus, high hydraulic heads, low hydraulic conductivity, and low depths of the geological limit would prevent from upconing, while the contrary would favour upconing and salt transport towards freshwater.

Saline paleo-water can upwell through aquitards and affect the quality of unconfined overlying coastal aquifers [39,40,41,42,43,44] (case (d), Figure 2). These saline waters are often found on low-lying coastal plains or islands where the wide range in permeability of the sediments may bring about the incomplete flushing of relict salinity. The upwelling of deep saline groundwater may be due to a hydraulic head higher than that of groundwater in the unconfined aquifer or may be triggered by head depression of fresh groundwater due to exploitation. In this case, the suitable indicators are the difference between the equivalent freshwater heads of the aquifers, the salt content of saline paleo-water, and the thickness and the hydraulic conductivity of the semi-confining unit. A comparison between the water levels of the two aquifers needs horizontal flow (static wells) in both aquifers and the knowledge of the density of both the fresh and the saline ground waters for the conversion of water levels in “equivalent freshwater heads” [45].

Another source of salts (case (e), Figure 2) is the process of over-wash related to storm surges that cause flooding with localized salinization [29]. More serious is the case of salinization caused by a tsunami, which can cause the inundation of large areas [46,47]. Both storm surge and tsunami produce the downward infiltration of seawater, having the water table as target. The indicators of low or high vulnerability are the coastal morphology (plunging cliffs or low topographic relief along the coasts), the topography (high topographic relief along the coasts or wedge-shaped bays), and the bathymetry of continental shelves (high or low seawater depths in front the shoreline) [48,49].

Coastal aquifers are affected by salt sources of other origins as well (case (f), Figure 2): evaporite salts, but also thermal saline fluids, sea water concentrated by permafrost salt exclusion, reverse osmosis, mobilization of reject of brines from desalination plants or hydrocarbon exploitation, and leakage of brines may cause water quality damage [30,31,38,50,51,52,53]. The vulnerability, as regards these sources, depends on their position compared to groundwater: the simplification of such contexts is feasible only at the scale of the local problem.

2.2. Limits of Indicators

The indicators have been outlined according to their most general meaning. Some of them, however, require specific considerations about their measure. The questions, relevant for all coastal aquifers, relate to salinity sources, “water level”, and quality (salt content) of groundwater. Other issues ensue from the anisotropy of hydraulic conductivity due to discontinuities and karst forms.

In the light of the variety of sources, to consider seawater as the only source of groundwater salinization is a limited view. The identification of salt sources different from seawater is often difficult [54], especially in case of source overlapping. However, their negative effects can sometimes outweigh those of seawater and their recognition becomes very important for the management since the remedies to be adopted can be very different.

Another problem concerns the “water level”, used instead of hydraulic head. Because of the frequent lack of specific monitoring nets, measures of “water level” are carried out in existing wells, whose completion and position within the flow system are often unknown. Often, wells are not in hydrostatic conditions, and water levels represent only an average of the hydraulic head along the well screens. The 3D variation of groundwater density in coastal aquifers is an additional factor that causes measured water levels to deviate from hydraulic heads. A thorough analysis of both the general concept of hydraulic head and its correct measurement in coastal aquifers is in Post et al. [45] and Post and von Asmuth [55].

Some methods for mapping groundwater vulnerability use the salt content of groundwater as indicator of its qualitative status. In a coastal aquifer, the three-dimensional distribution of salinity depends on aquifer heterogeneity [56] and on the historical evolution of exploitation and its current regime. If salt content refers to samples obtained by pumping wells, it will represent the mean value of the groundwater thickness influenced by pumping, thus inducing large errors in the evaluation of real qualitative status of groundwater. Correct information on salinity stratification derives only from the direct measure (electrical conductivity logs) in hydrostatic and totally screened wells located far from discharge zones. Wells close to coastline do not give correct information because of the upward flow of the discharge zone [57,58]. In case of vertical currents, the correct information at each elevation could be gained only by using packers to isolate all the well surface out the aquifer thickness to be sampled [59]. The average salinity referred to a defined groundwater thickness could be a reasonable indicator.

Karst aquifers show both anisotropy and heterogeneity. Their hydrogeological features depend on their structure [60,61]. Faulting, fracturing, and karst processes change the porosity (and permeability) of formations, which has a great impact on the patterns of migration and storage of underground fluids, where groundwater flow depends on arrangement, connectivity and opening of discontinuities and karst forms. In karst coastal aquifers, the mixing of freshwater and saltwater triggers water-rock interaction, involving calcite dissolution and precipitation, and dedolomitization, that cause a further change in porosity over the geological time [33,62,63]. Coupled models of dissolution and transport show that dissolution potential depends on the mixing rate of freshwater and saltwater [64]. Fully coupled reactive transport models [65] demonstrate that dissolution takes place along the entire mixing zone towards the freshwater side, while maximum dissolution occurs near the coastal outflow face. The widening of the transition zone, brought about by spatial heterogeneity and the temporary effects linked to density variations, may cause local intensification of mixing and karstification [66]. The coupling of tidal oscillation with heterogeneity can, instead, reduce mixing [67]. In coastal cave genesis, dissolution generated by fresh and sea water mixing is higher than dissolution by meteoric water infiltration [68]. The dissolutional potential of brackish water in the transition zone of freshwater lenses is at the base of the “flank margin model”, which explains the development of dissolution caves along the coast of the Bahama Islands [69]. Saline/freshwater-mixing processes during sea-level high-stands are known as one the factors in the spelogenesis of Puerto Princesa Underground River (Philippine) [70].

Fault zones in carbonate rocks are characterized by wide zones of fractured rock (the so-called “damage zone”) that accompany narrow fault zones (localized deformation zones, “fault-core”). In carbonate rocks, secondary porosity for fractures and faults can be further modified by diagenetic and karstic processes: fault and fracture areas are places where dissolution and precipitation, both capable of altering the permeability structure, occur. If carbonate rocks are massive, fracturing and dissolution can cause the increase of secondary permeability such that storage capacity is determined by the only volume stored in the fault areas. The hydrogeological behaviour of fractured rocks in general is difficult to determine, as tectonics results in anisotropic permeability structures in correspondence to narrow strain bands. Faults can behave as barriers to cross-fault flow or as effective conduits for fault-parallel flow: however, small changes in the fault zone architecture and hydraulic parameters can largely modify the hydrogeological behaviour of faults [71,72,73,74,75,76,77,78,79,80,81].

In karst aquifers, fissures, fractures, faults, surface karst forms, and endorheic basins influence recharge processes and pollutant transport from the ground surface. For karst aquifers (case (a), Figure 2), the density of dolines [82] (also named sinkholes [83]) and the presence of endorheic basins may be indicators of focused recharge. The density of fissures, fractures, and faults may indicate vertical preferential mass transport, provided that the hydraulic behaviour of related permeability structures is known. Because of reactivation of basement faults or recent tectonics, faults and sinkholes can interest the covers lying on karst basements, with the chance for infiltration water to by-pass the covers. Another specific indicator of vulnerability is the presence of the epikarst [84]. It has an important role in the transfer of pollutants from the ground surface to the saturated zone [61]: on one hand, it can store and concentrate recharge flows, thus increasing the natural protection of groundwater against reactive pollutants (i.e., subject to processes of degradation and abatement) and, on the other hand, it can transmit stored and concentrated non-reactive pollutants under the effects of extreme events or high recharge flows, thus increasing the intrinsic vulnerability of the aquifer [85].

General considerations about “water level” apply to karst aquifers as well, but the hydraulic head measures in karst aquifers are even more difficult than in other aquifers because of anisotropy of hydraulic conductivity. In fact, hydraulic head distributions in karst and karst coastal aquifers depend significantly on the position of high- or low-conductivity elements and on the permeability structure [86]. Heterogeneous and anisotropic permeability distribution and non-constant fluid density limit the use of water levels as in homogeneous and isotropic aquifers [87]. Thus, field hydraulic data and potentiometric maps must be considered with great caution, because their interpretation needs specific approaches [88].

The permeability along the coasts depends on the density of fissures and fractures, and from karst development. The uneven distribution of permeability allows either focused (coastal and submarine springs) or diffuse discharge that, depending on hydraulic condition, ease or impede the ingression of seawater inland [89,90,91,92]. The paleo-karst features (as the conduits of submarine springs and the horizontal karstification levels [93,94,95] may enhance local communicability with seawater along the subterranean estuaries [96,97]. For case (b) (Figure 2), the distribution of fractures, faults, and karst structures along the coast together with their hydraulic behaviour become specific indicators. For case (c) (Figure 2), the radial horizontal (K_r) and the vertical hydraulic conductivity (K_z), and the anisotropy ratio (K_z/K_r) will favour/inhibit the upconing.

3. Material and Methods

Salento is a peninsula located in the Apulia region (southern Italy). It extends from NW to SE between the Ionian Sea and the Adriatic Sea (Figure 3a).

In the 1960s a pioneering net of 14 observation wells (OW) was set in the Salento aquifer with the aim of studying seawater intrusion processes [98]. The OWs are special wells that include about 50 m of saltwater under the transition zone; well screens extend over the whole saturated zone thickness. From 1970 to 1986 monitoring of variable frequency involved temperature (T) and electrical conductivity (EC) logs, while a few sampling surveys allowed the chemical features of fresh-, brackish- and salt waters at different depths to be ascertained. In 1994, the Apulia Regional government included the OW net of Salento in a larger net consisting of a total of 120 wells; from 1995 to 1996 groundwater sampling, multi-parametric logs and water level measurements were carried out with quarterly frequency [99,100]. From 2007 to 2011, monitoring followed a half-yearly frequency on the same net [101]. From 2016, the Salento net includes 35 wells with a half-yearly or four-month frequency monitoring depending on parameters.

T data from the multi-parametric logs related to the regional monitoring net of winter 1995 allow the reconstruction of horizontal and vertical thermal sections of Salento aquifer at regional scale [89]. Such sections gave a regional outline of the thermal field. T probes have a range between −5 °C to 50 °C, accuracy of ±0.20 °C and resolution of 0.01 °C; temperatures were recorded at 0.5 m intervals. The reconstruction along vertical sections requires more data than horizontal sections because of the increase of scale (medium). In karst aquifers the isotherm trend may give information on location of recharge areas, permeability variations, preferential flow paths, role of discontinuities, location of conduits, connection with overlying aquifers, and groundwater–surface water relationships [85,102,103,104,105,106]. The rationale is that the hydraulic field disturbs the conductive thermal field by forced advection, leading to a temperature distribution, which is the result of a time-dependent exchange of heat between water and surrounding rocks.

Consequently, aquifer thermal fields can be used as analogues of flow systems. The most significant disturbances are observed in fractured/karst aquifers, characterized by high K_v/K_h anisotropy ratio of the hydraulic conductivity. The significance of the hydrogeological disturbances depends on topography, climate, 3D water table configuration, basin geometry, magnitude and spatial distribution of permeability, hydraulic anisotropy and depth of active flow, and thickness [107,108,109]. In coastal aquifers the trend of isotherms may give additional information on the spatial extent of seawater intrusion and saltwater upconing, provided that the effects of solar radiation in case of very shallow groundwater, geothermal anomalies, and human water injections are considered. In this way, the thermal field can help clarifying some vulnerability issues.

The OW LR (location in Figure 3a; the code LR stands for Red Lake, where the toponym probably derives from local old bauxite quarries) is selected for the following elaborations because of the longest series of monitoring data and since it is in a zone with horizontal groundwater flow. This last feature allows the well to reproduce the density stratification of groundwater [57,105]. EC logs allow the measurement of freshwater thickness, position and thickness of transition zone, saltwater top, and average EC of fresh and transition thicknesses. The lab measures of EC, TDS (Total Dissolved Solids), and density on “standard” solutions obtained by mixing freshwater and saltwater sampled in the OWs allow the EC-TDS-density calibration. Density profiles allow calculating the elevation of the theoretical sharp interface (TSI) and the thickness of the complementary fresh water column (FWC) [110]. TSI is used as synthetic indicator of the TDS vertical distribution. It is calculated on the basis of TDS profiles, by splitting the whole water column overlying the top of saltwater in two components: a freshwater component having the lowest salinity measured in time in the observation well, and a saltwater component having the salinity of the saltwater found at the bottom of the well. The addition of the thickness of saltwater component to the measured elevation of the top of saltwater (identified through the analysis of the TDS profile) gives the elevation of the TSI. This is equivalent to considering freshwater and saltwater as non-miscible fluids. The decrease or rise of the TSI elevation, respectively, indicate the increase or decrease of the thickness of the freshwater component. The time series of the TSI and the FWC suffer of a few uncertainties in their calculation due to variation over time of EC probes, calibrations, and operators.

Data on yearly effective infiltration, yearly irrigation, and groundwater stress are from Portoghese et al. [111]. Chloride concentrations in groundwater for the period 1970–2015 come from the monitoring records of AQP (Acquedotto Pugliese) that distributes drinking water in the Apulia region. Data on water levels (Figure 3a) refer to measures in the wells of the regional monitoring net during winter 1995. The geological and geomorphological data are from the online databases of Apulia Regional Government, and the River Basin District of southern Italy, with some integrations of published scientific data.

Geologic and Hydrogeological Setting of the Study Area

The coastal karst aquifer of Salento belongs to the Apulia carbonate platform, a succession of Jurassic–Cretaceous carbonate rocks from about 3 to 5 km thick, which is the upper part of the Apulia foreland emerged at the end of Cretaceous [112,113]. Cretaceous limestone and dolomitic limestone form the geological basement, which outcrops in large areas with patchy covers of clay, sand, and calcarenite of Miocene to Pleistocene age [114] (Figure 3a). To the NW, a large tectonic scarp W–E trending (North Salento Fault Zone, NSFZ) separates Murgia from Salento [115].

The basement shows structural highs and lows separated by major normal faults (sub-vertical normal, strike-slip and oblique-slip) of Plio-Pleistocene age, striking NNW–SSE and subordinately NW–SE with maximum extensional offset of 200 m (Figure 3b). Brittle deformation during Middle and Late Pleistocene caused rare faults with small displacements, and joints [116]. The fault system determines in the southernmost part of the peninsular rhomboidal basins and ridges. The SW Salento shows major folds with NNW–SSE to NW–SE trending axes, while minor folds are ENE and NE trending. The karst surface shows karst plains, fracture zones, dolines (1452 from a recent census [117]), and hundreds of endorheic basins of fractal dimension (Figure 3b). Many surface karst landforms are significantly controlled by tectonic discontinuities [118,119]; many dolines are found on covers and align along Quaternary faults [116]. Figure 3b shows only part of dolines because of their high density (from 0 to 7 per km²). Subsurface karst forms develop both along the vertical discontinuities and in near-horizontal planes according to lithological-stratigraphic heterogeneity and the paleo-geographic history.

Precipitation is 638 mm/year, while the effective infiltration is 132 mm/year (yearly average 1960–2010, [111]). The flow system is regional. Fresh groundwater floats on saltwater as a lens with a salt content varying between 0.2 and 0.5 g/L; groundwater discharges to the sea through brackish coastal springs whose mean salt content varies in the range 3.5–20 g/L. The interconnection of preferential horizontal levels, vertical karst forms, and major and minor discontinuities cause high anisotropy of the hydraulic conductivity; the permeability is high at regional scale. Groundwater is locally confined by the presence of low permeability carbonate levels or because of the dislocation of the carbonate basement under mean sea level. The water levels (mean values) at regional scale (Figure 3a) vary between 5 m AMSL (above mean sea level) in the NW and SE sectors of the Peninsula and less than 0.5 m AMSL close to the coastline. The mean hydraulic gradient is 0.2‰. Figure 3a shows at NW borders water levels higher than 5 m AMSL; these levels highlight the groundwater flow coming from the adjacent karst coastal aquifer of Murgia, which recharges the Salento aquifer with about 9 m³/s [120].

Because groundwater is the only available water resource for drinking purposes and irrigation, in the absence of freshwater springs the exploitation takes place by more than 100,000 (official and unofficial) wells, with unavoidable consequences of groundwater pollution and salinization [121]. A decrease of groundwater discharge and an increase of groundwater salinization are expected in relation to sea-level rise [122].

4. Results and Discussion

Figure 4 refers to an area of the Salento Peninsula (borders of the area in Figure 3) that is characterized by a high density of dolines, fractures and faults, with important basement disruption. The area is of fundamental importance for the drinking water exploitation, which is why the information on vulnerability and protection areas is required. The track AA’ (Figure 3) with a horizontal distance of 48.8 km, is selected for the reconstruction of a cross-section of the thermal field.

The lithological map of the selected area (Figure 4a) shows the location of wells, while Figure 4b shows the endorheic basins (outlined as macro-areas), dolines, and tectonic lineaments. Figure 4c outlines the elevation of the top of the carbonate basement and the (also inferred) normal faults. Figure 4d shows the water level profile from the slicing of the contour map of Figure 3a. Figure 4d also shows the horizontal distance of the wells used for the thermal section of Figure 4e, along with the value of the average TDS of ground waters from EC logs. These TDS values refer, indeed, to different saturated thicknesses explored by each well.

Figure 4e shows the vertical cross-section A-A’ of the Salento aquifer thermal field. It is based on interpolation of T logs carried out during winter 1995 in 4 deep wells (well codes P1UGE, P1TAU, 6-FEOGA, and 14RFOT) of the regional monitoring net and other 5 logs (well codes 7TAU, 3TAU, 1SAN, 2PALM. 3BAGN) from a database consisting of 380 T logs related to wells of the Irrigation Consortium “Ugento Li Foggi” of Salento. Dotted lines show the penetration of the wells in the saturated zone. The data interpolation considered also the consistency with other thermal cross-sections (not shown), whose tracks cross the AA’ track.

In the central part of the section, cold waters (with a temperature between 15 °C and 17 °C) manifest their influence up to −180 m AMSL. The convexity of isotherms is downward oriented, pointing at a relation of cold waters with the ground surface. The cold area matches a basement depression (a horst with grabens at both sides). The horst-bounding faults intersect the Quaternary and Miocene-Plio-Pleistocene covers and many deep dolines align with fault tracks. Faults and dolines locate in endorheic basins. To the left of the cross-section, the downward inflection of the 17.2 and 17.4 °C isotherms suggests a recharge flow coming from the outcropping dolomite and limestone. At about 13 km of horizontal distance, an upward inflection of the 17.4 isotherm separates this recharge flow from the middle one, suggesting a rise of saltwater. A few main faults cross the basement and the ground surface shows a lower density of dolines compared to the middle area. To the right of the section, the downward inflection of the 17 °C isotherm outlines another recharge area, whose influence reaches the −50 m AMSL. Here, next to an outcropping limestone horst, the basement top is between −100 m and −120 m AMSL. The corresponding covers are not crossed by faults and show only fracture zones in endorheic basins. An upward inflection of the 17.4 °C isotherm to the left of the horst suggests another rise of saltwater. The groundwater TDS along the track AA’ indicates, in fact, that the downward inflection of isotherms corresponds to low salinity ground waters, while the upward inflection matches warmer and higher salinity waters.

The inflection or the steps of hydraulic gradients from water level profiles may provide a first estimate of the hydraulic properties of faults [77,123]. The available water level profile for the section AA’ derives from data on a small scale, which do not focus specifically on fault zones. Even with these limitations and considering what was said about the limits of indicators, we would use this water level profile to infer the hydrogeological role of faults in the area. At about 14 km of horizontal distance, close to the P1TAU well, the water level starts decreasing towards right (Figure 4d). This might indicate the presence of a low conductivity fault. The hydraulic gradient changes again at about 32 km of horizontal distance, suggesting the presence of another low conductivity fault. The upward convexity of the isotherm 17.4 °C between the 3TAU and P1TAU and between the 2PALM and 3BAGN wells suggests an anisotropic permeability distribution within the low-conductivity faults, which at the same time may favour the upconing of saltwater and prevent infiltration from the surface. The decreased water levels in the central part of the profile suggest the presence of a draining zone between two local flow systems (Figure 3a and Figure 4b). Thus, the faults in the central zone might behave as barrier to horizontal flow along the track AA’ and as effective conduits for fault-parallel flow; along this direction the faults favour the mass transport from surface.

In the central part of the section the covers do not prevent the direct infiltration to the underlying groundwater. The faults in the outcropping limestone to the left of the section show from medium (recharge area) to low conductivity (before the highest water level), while higher-conductivity faults correspond to water level decrease. The successive outcropping limestone horst is bounded to the left by low conductivity faults, while faults show a higher conductivity to the right (recharge area). The next covers are devoid of dolines, and faults are in the buried basement. The temperatures of recharge areas at both sides of the section are higher than in the central area. The temperature of the coldest fresh waters is between 15 °C and 16 °C, close to the mean atmosphere temperature during the recharge period (12 °C during the autumn-winter season). These temperatures indicate high permeability of formations, while warmer fresh waters (T > 17 °C) correspond to less permeable zones.

Figure 4e shows the elevation of the predicted Ghyben–Herberg (GH) interface compared to the isotherm trend. The temperature of saltwater in the Salento aquifer varies with latitude. In correspondence with the section, the available observation wells show that saltwater temperature is permanently between 17.8 °C and 17.9 °C. The 17.8 °C isotherm ends near the coastlines, thus outlining the top of the saltwater. The predicted GH interface gives a misleading indication, which opposes to the TDS actual stratification found in the wells of the section and to the indication of the 17.8 °C isotherm. Differences, as a first approximation, can be attributed just to the hydrogeological role of faults. The gap between the real salinity distribution and GH elevation occurs even in the absence of significant anisotropy and vertical hydraulic gradients.

Figure 5a shows the TDS profiles related to the LR observation well (in Figure 3a), located inland, far from the discharge zones. The profiles show a transition zone of variable thickness, spanning from −70 to max −160 m AMSL (max thickness of 90 m), which moves vertically over time. Figure 5e shows the elevation of the GH interface, predicted using the actual densities of both the freshest water and saltwater measured in the well using the water level. This water level (Figure 5e) is actually an “environmental head” (EH) [124] and not an hydraulic head because it refers to a column of variable density: in this case, being the water column of LR in stagnant conditions (without vertical currents), the water column inside the well represents the condition of the aquifer outside the well. Actually, about a half of the water column above the GH interface is brackish, and the thickness of freshwater is smaller than predicted.

Figure 5e compares the elevation of this GH interface with the elevation of the TSI. TSI is a more realistic synthetic index of density distribution than the predicted GH interface because it considers the real groundwater density stratification. However, the predicted fresh water column does not correspond to the reality as well. The use of the “equivalent freshwater head” (EFWH) instead of the EH (EFWH is higher than EH) in predicting GH and TSI causes the GH interface elevation to decrease significantly, while TSI elevation changes little. The gap between the GH and TSI increases when adding the correct formulation of the hydraulic head.

TDS profiles and calculated parameters obviously vary with time because of the evolution of FW-SW equilibrium under climate and human pressure variations. Figure 5c shows the groundwater stress (GW) index, obtained as the ratio of the yearly irrigation and the effective infiltration (Figure 5d, data from Portoghese et al. [96]). It shows high values during the droughts around 1990 and 2003 [125,126]. Between 1987 and 1991, the GW stress increases up to 0.7 (severe stress) due to the concurrent decrease of the effective infiltration and the increase of irrigation. This triggers upward salt transport (LR well is located inland) and the salinization of fresh groundwater.

Another peak of GW stress (moderate stress) occurs between 2000 and 2004. The TSI of 1995 and 1996, following the first high GW stress, moves upward about 13 m compared to 1974. Despite the following increase of recharge, TSI shows only a modest downward shift. After the second period of GW stress, the TSI does not vary significantly.

Unluckily, the monitoring of the LR well has two wide time gaps, which occur just during the drought periods. More detailed monitoring relates to the drinking water well net. Chloride concentrations (Figure 5b) measured in groundwater samples drawn under pumping from four representative wells (location in Figure 3a) of this net show a significant increase during the two drought periods with a certain delay compared to groundwater levels decrease. All data show that the succession of droughts and the continuous exploitation do not allow recovering the pristine water quality of the 1970s.

5. Conclusions

Temperature of groundwater demonstrates its potential as natural tracer of groundwater flow, recharge areas, anisotropy of hydraulic conductivity, and seawater/saltwater intrusion. The downward and upward convexity of isotherms outline the zones of recharge and upwelling of saltwater, respectively. The groundwater salinity distribution validates this association. The former trend provides information on the intrinsic vulnerability (likelihood of mass transport from the ground surface), while the latter relates to the vulnerability of groundwater to salinization, both inland and close to the coastline.

The trend of isotherms may be associated with the permeability features of faults as well. In this case the validation of the hypotheses suffers from the lack of detail of water level distribution and uncertainty about the meaning of water levels in a karst medium. Anyway, the changes of the hydraulic gradient matching the upward isotherm convexity suggest a low conductivity of corresponding faults, while the decrease of water level in connection with downward convexity zones may represent faults of higher conductivity. The coupling of hydraulic and thermal information allows hypothesising that the faults crossing the covers, corresponding to the buried karst, behave at the same time as barrier for horizontal flow along the section track, and as conduits for fault-parallel flow. On the basis of these findings, we can infer that groundwater under the exposed karst (characterized by some low conductivity faults) is less vulnerable than groundwater in the faulted and displaced basement, overlaid by a significant thickness of covers crossed by high conductivity faults. The difference in vulnerability between the exposed and the covered karst is also marked by a different density of dolines. A high density of dolines occurs in the area of the covered karst: they are often aligned with faults [127]. These alignments have been observed in other regions of the world in correspondence with active faults linked to extensional tectonics [128]. However, the role of faults in the vulnerability issues of karst aquifers depends on the permeability structure of each fault and this prevents any generalization on their use as vulnerability indicators. A low or high density of dolines might be a good indicator of low and high intrinsic vulnerability, respectively. However, a correct assessment of doline typology should be undertaken before using their density as vulnerability indicator [82].

The elevation of the 17.8 °C isotherm can be used in the study area as a reliable homolog of the elevation of the bottom of the transition zone. In Salento, the saltwater temperature varies little with latitude, but it does not fluctuate over time in a same zone; reference values can be assessed with suitable observation wells. At other latitudes, saltwater may show different temperatures.

The reason why temperature logs are more effective and reliable than EC/TDS logs in recognizing the salinity stratification is that, always providing that wells are in static condition, EC logs represent groundwater outside the well only coinciding with well-screens, while temperature logs give reliable information on the “ambient temperature” even along the casing. Groundwater, indeed, flows around the casing, but water within the borehole easily equilibrates its temperature with that of the aquifer. Based on this fact, we could measure temperature in wells with vertical currents by using a liner system to seal temporary the boreholes: in stagnant wells, the “ambient temperature” correctly represents the fracture structure of the aquifer [129,130].

The information on the salinity stratification can be gained only in zones of horizontal flow through observation wells open to the aquifer and crossing its thickness up to saltwater. Measures should be then repeated over time at these elevations to obtain information on the evolution of salinity. The transport–density-dependent flow models might help to reconstruct the 3D distribution of salinity, but we must solve the issues linked to the anisotropy of the karst medium. Even neglecting the anisotropy, such models require a thorough hydrogeological knowledge and numerous field data to set the assumptions and validate their results [131].

The use of water levels as indicators of vulnerability in karst coastal aquifers appear questionable in the light of the variable permeability of rock mass and the hydrogeological influence of faults. We used the water level profile to validate hypotheses about the role of faults, but the uncertainty on water levels may be considered low because of the scale effect on permeability. The position of the GH interface compared to temperature trends and salinity data highlights the questionable information of water levels when used to predict the position of saltwater and, consequently, the vulnerability to salinization.

The results of the study, even if limited to a part of the Salento aquifer, suggest that the definition of the three-dimensional vulnerability of a coastal karst aquifer is not immediate, and the object cannot be achieved on large scales on the basis of an oversimplification of the issues at stake. Results shows how discontinuities and karst landforms play a predominant role in determining vulnerability. This role is not unique and must be understood zone by zone before deciding about vulnerability issues. Temperature, combined with a coherent physics of flow and mass transport and site-specific “hydrogeological knowledge”, can be of great help in this recognition. The future development of the study will consider information about stable isotopes and water chemistry in the search for additional indicators of vulnerability. Finally, it should be noted that the groundwater vulnerability to salinization is not invariant over time as the intrinsic vulnerability, because it depends on the qualitative state of groundwater, which, in turn, depends on previous states because of the non-linear behaviour of coastal aquifers.. This requires fit-to-purpose monitoring to evaluate the time-evolution of all necessary elements.

Unfortunately, the proposed methodology requires measures along wells, and this limits its direct application to karst aquifers where only springs are available for monitoring.