Hydrostratigraphic Framework and Physicochemical Status of Groundwater in the Gioia Tauro Coastal Plain (Calabria—Southern Italy)

Abstract

In this study, we analysed the Gioia Tauro Plain (Tyrrhenian coast, southern Italy) in terms of hydrostratigraphy and the physicochemical status of groundwater. We investigated the hydrostratigraphic framework of the area identifying a deep aquifer (made by late Miocene succession), an aquitard (consisting of Pliocene clayey and silty deposits) and a shallow aquifer (including Late Pleistocene and Holocene marine and alluvial sediments) using subsoil data (boreholes and geophysics). Our reconstruction showed that the structural geology controls the spatial pattern of the aquitard top and the shallow aquifer thickness. Furthermore, we evaluated the hydraulic conductivity for the shallow aquifer using an empirical method, calibrated by slug tests, obtaining values ranging from 10⁻⁴ to 10⁻⁵ m/s with a maximum of 10⁻³ m/s located close to inland dune fields. The piezometric level of the shallow aquifer recorded a significant drop between the 1970s and 2021 (−35 m as the worst value). It is the effect of climate and soil use changes, the latter being the increased water demand for kiwi cultivation. Despite the overexploitation of the shallow aquifer, shallow groundwater is fresh (736 µS/cm as mean electrical conductivity) except for a narrow coastal area where the electrical conductivity is more than 1500 µS/cm, which can be due to the seawater intrusion. What was more complex was the physicochemical status of the deep aquifer characterised by high temperature (up to 25.8 °C) and electrical conductivity up to 10,520 µS/cm along the northern and southern plain boundaries marked by tectonic structures. This issue suggested the dominant role of the local fault system that is likely affecting the deep groundwater flow and its chemical evolution.

1. Introduction

The sustainable development and management of coastal plains have received significant attention worldwide during the last few decades (see, e.g., [1,2,3,4]). However, coastal areas and coastal plains suffer the adverse consequences of geological processes and climate change hazards, including rises in sea level [5]. Moreover, the increasing human-induced pressure worsens the impact of land-use and hydrological changes, including reclamation works, embankments, dams, and groundwater overexploitation [6]. The latter is a key issue in coastal governance, since it must face growing population demand. Thus, salinisation processes, subsidence, and/or flood damages are favoured [7,8,9,10,11,12], and in this respect, the Mediterranean coastal areas are prone to seawater intrusions [13,14,15,16,17].

The study of the underground relationships among fresh and salt water, groundwater flow, and geological and aquifer features is still a scientific methodological issue to be improved for a wide range of coastal plains and coastal aquifers. Therefore, this study aims to offer a contribution, based on hydrogeological and geochemical data, to investigate a geologically complex coastal plain located in Southern Italy with a multidisciplinary approach.

The Gioia Tauro Plain (GTP) is positioned along the Tyrrhenian coast of the Calabria Region. GTP is an important industrialised and coastal agricultural area of about 500 km². The harbour of Gioia Tauro is one of the most important container traffic hubs in the Mediterranean basin and the largest transhipment terminal in Italy. According to [18], about 160,000 inhabitants live in this area. However, the Plain still presents portions completely devoid of basic services, such as sewers or connections to aqueducts for domestic use. The lack of such infrastructures increases the overexploitation of groundwater. Due to well-developed anthropic activities and the absence, in some areas, of primary services, the GTP groundwater resources are constantly at risk of depletion and quality degradation.

Furthermore, GTP is characterised by a very articulated geological and structural framework that, mainly along with the bound areas, favours the development of water chemistry not fully explainable by simple water–rock interaction processes occurring in the sedimentary aquifers, and seawater intrusion. As a result, warm and saline waters were observed close to the northern and southern boundaries of GTP [19,20]. In addition to the ongoing anthropic pressure (mainly due to intensive farming), this issue makes the GTP groundwater system of high interest.

The main aim of the work is the definition of the hydrostratigraphic framework and the physicochemical status of GTP groundwater. We provide a detailed characterisation in terms of the hydrogeological framework, hydraulic properties, and geochemical processes occurring in the study area.

2. Geological and Hydrogeochemical Background

From a geological point of view, GTP is an Upper Pliocene-Pleistocene NNE-SSW oriented half-graben [21]. Three main fault systems mark the boundaries: (i) the NW-SE striking Nicotera-Gioiosa fault zone (NGFZ) with right lateral kinematics to the north, (ii) the high-angle normal and active NNE-SSW Cittanova Fault (CF) to the east, and (iii) the NW-SE trending Palmi-Locri fault zone (PLFZ) with mainly strike-slip kinematics to the south [22]. Furthermore, toward the coastline, GTP is crossed by the Palmi-Gioia Tauro Fault (PGF), a high-angle NE-SW striking fault west-dipping, and another east-dipping parallel fault, whilst the NE-SW oriented Sant’Eufemia-Laureana Fault (SLF), west-dipping, cuts the middle sector of GTP (Figure 1) [22].

Along CF, the igneous-metamorphic bedrock, building the Aspromonte and Serre Massifs, is juxtaposed with the sedimentary infill of the GTP basin [22,23]. CF represents an active fault as testified by the fault rupture that caused catastrophic earthquakes in the 4th century AD and 1783 (Mw = 7.0) [24,25].

In the southernmost sector, GTP is bounded to the west by the Palmi high, made by igneous-metamorphic bedrock, uplifted and eastward-tilted [23]. The presence of a structural high in the GTP western sector is also testified by the mechanisms of antecedence visible in the terminal segment of the Petrace and Budello rivers. At the same time, the tilting toward the east of the GTP is highlighted by the east-dipping of the recent sediments (Figure 1) [26].

The igneous-metamorphic bedrock is unconformably covered by an upper Miocene siliciclastic (including both shallow- and deep-water marine sediments and fluvio-deltaic deposits) and carbonate (shallow marine) succession [27], which crops out mainly along the Palmi high. The GTP basin is filled by a 600 m thick Plio-Pleistocene marine sedimentary succession [23]. The bottom of this succession is made by Pliocene clayey and silty deposits (Trubi Formation in [27]), 100 m thick [23] with a pumice-rich horizon in the upper part [28]. Upward, cross-bedded sands and calcarenites (late Pliocene–Pleistocene Calcareniti di Vinco in [27]), representing the infill of the Siderno paleo-strait [29], cover, with a slight angular unconformity, the Trubi Formation [23,27]. In the eastern sector of GTP, the Calcareniti di Vinco are covered, along an erosive contact, by the Sintema di Taurianova (Late Pleistocene) consisting of gravelly and sandy alluvial fan deposits characterised by lateral and vertical heteropic facies [27]. Toward the top, the succession is closed by sandy to gravelly terraced marine deposits (Late Pleistocene–Holocene), up to 20 m thick aeolian inland dunes (Holocene) and recent alluvial (in the eastern and middle sector) and coastal (in the western sector) deposits [27,30].

From a geochemical point of view, previous surveys [19,20], with the aim of identifying possible implications for fluid-fault relationships, provided data about the geochemical characteristics of the GTP groundwater. As a result, different geochemical water types and important anomalies in terms of total dissolved solids (TDS) and temperatures in the northern portion of the GTP, at the intersection between the main fault systems, were recognised. As reported in Figure 2, the Refs. [19,20] identified five main hydrochemical facies in the GTP: Ca-HCO₃, Na-HCO₃, Ca-SO₄, Na-Cl (high salinity and high temperature) and Ca-Cl. The authors attributed the presence of several water types to different processes: water rock interaction occurring in the sedimentary aquifers, Messinian halite dissolution, marine intrusion, and a prolonged time of water-rock interaction with the crystalline basement.

3. Methods

The hydrostratigraphic analysis was based on stratigraphic data integrated by field surveys. Borehole data (85 as total) were acquired from an open-access database of boring and well data [31], technical sheets of the analysed wells and scientific papers [28,32]. The open-access database and technical sheets supply, in addition to the stratigraphic information, the piezometric level and, in some cases, groundwater temperature.

The lithological information from boreholes was converted into hydrogeological units (HU), mainly in terms of grain size. The latter is defined as a geologically homogeneous unit (but not necessarily isotropic) with characteristic hydrodynamic properties [33].

Considering the good efficiency of the electrical geophysical prospecting in the reconstruction of the Quaternary sediment geometries (e.g., [34,35]), borehole data were integrated by 133 vertical electrical soundings (VES) inferred by SP26 [36]. Boreholes and VES data were used to define the top and bottom elevation of the Hus and referred to as meters above sea level (m asl).

According to static and dynamic piezometric levels available for 201 wells [31] (Table S1) and measured in 24 wells (in this work, Table S1), we obtained values of hydraulic conductivity (k) using an empirical method, which is a simplification of an analytical solution [37]. It uses the following equations:

$$\mathrm{T}=1.22\times \mathrm{Sc}$$

(1)

$$\mathrm{Sc}=\mathrm{Q}/∆\mathrm{h}$$

(2)

$$\mathrm{T}=\mathrm{k}\times \mathrm{b},$$

(3)

where T is the transmissivity (m²/s); Q is the well discharge (m³/s); Δh is the drawdown (m); Sc is the specific capacity (m²/s), equal to Q/Δh; and b is the saturated aquifer thickness, drilled by each well (m). To validate the k-values obtained by the empirical method, a slug test survey was realised using 12 selected wells (diameter between 0.25 and 0.40 m) with lithostratigraphic boring data. The slug test consists of a rapid increase (rising-head test) or decrease (falling-head test) of the static piezometric level into the low-diameter borehole, usually a piezometer; this effect is usually obtained by quickly adding/removing a volume of water or quickly throwing/removing a slug, respectively [38]. We realised that falling-head tests using a high-power submersible pump were able to produce a very quick drawdown (together with a sufficient withdrawal volume) beneath the static piezometric level. This modification of the common-used procedure was adopted to solve the practical problems due to the large diameter of the available boreholes; the pump usage was calculated to be very short and quick (less than 30 s). The following return of the water level to the initial position was monitored. We applied the Hvorslev method [38], measuring the maximum drop of the water level (H0) and the subsequent relative recovery (HU1, …, Hn).

The load ratio (H/H0) and the relative times, obtained by the tests, were plotted in a semi-log (ln) plot, and the tl (time lag) was taken as the time which corresponds to H/H0 = 0.37.

Considering that the volume (V) needed to balance the active load is:

$$\mathrm{V}=\mathrm{A}\times \mathrm{h}$$

(4)

$$\mathrm{tl}=\frac{\mathrm{V}}{\mathrm{Q}}=\frac{\mathrm{A}\times \mathrm{h}}{\mathrm{F}\times \mathrm{K}\times \mathrm{h}}=\frac{\mathrm{A}}{\mathrm{F}\times \mathrm{K}}$$

(5)

$$\mathrm{K}=\frac{\mathrm{A}}{\mathrm{F}\times \mathrm{tl}}$$

(6)

F, representing a shape factor, was calculated by the following equation, which complies with the condition of the tested wells characterised by filters distant from the aquifer bottom:

$$\mathrm{F}=\frac{2\times \mathsf{\pi }\times \mathrm{L}}{\ln \left(\frac{\mathrm{L}}{\mathrm{R}}+\sqrt{1+{(\frac{\mathrm{L}}{\mathrm{R}})}^2}\right)}$$

(7)

where k is the hydraulic conductivity (m/s); h is the dynamic level (m); A is the basal area of the tested sector (m²); L is the length of the filters (m); R is the radius of the filters (m); tl is the time lag. The k-values obtained by the 12 slug tests and the relative values calculated by the empirical method are reported in

Table 1. K-values obtained using slug test and pumping test data in correspondence to the same wells.

Slug Test	k (m/s)	k (m/s) Pumping Test
ST1	9.85 × 10−4	3.98 × 10−3
ST2	3.73 × 10−4	8.27 × 10−4
ST3	7.07 × 10−4	1.85 × 10−3
ST4	5.70 × 10−5	3.82 × 10−4
ST5	2.20 × 10−4	8.56 × 10−4
ST6	5.04 × 10−4	1.33 × 10−3
ST7	1.58 × 10−3	5.43 × 10−3
ST8	3.94 × 10−4	1.17 × 10−3
ST9	2.41 × 10−4	3.46 × 10−4
ST10	6.77 × 10−4	2.76 × 10−4
ST11	1.14 × 10−3	4.04 × 10−3
ST12	3.94 × 10−4	9.43 × 10−4

.

The methodology of groundwater characterisation was based on groundwater level and physicochemical parameters, as described by the refs. [39,40]. The water level and temperature, pH, Eh, and electrical conductivity (at 25 °C) were measured during March–April 2021 sampling campaigns by a phreatimeter and multiparametric probe (Hanna Instruments—HI9828), respectively. In the phreatic wells, the physicochemical parameters were measured directly in situ. Two pH buffers, with nominal pH values of 4.01 and 7.01 at 25 °C, were used for pH calibration. In addition, ZoBell’s solution [41] was used to calibrate the mV-meter for Eh measurement. The survey consisted of 101 wells with depths ranging from 4 to 200 m: 83 wells from the shallow phreatic aquifer, and 18 wells from the deep confined aquifer. The main statistical values referring to groundwater physicochemical features of surveyed wells are listed in

Table 2. The main statistical features of the physicochemical parameters of groundwater belonging to the main aquifers (83 and 18 wells for HU3 and HU1, respectively) characterising the GTP. EC: electrical conductivity; TDS: Total Dissolved Solids; p.s.u.: practical salinity unit (as ‰).

		T	pH	EC (25 °C)	Eh	TDS	Salinity
		(°C)		μS/cm	mV	ppm	p.s.u.
HU1	Min	14.8	7.5	230	−253	115	0.1
	Max	25.8	8.9	10,520	238	5264	5.9
	Median	18.2	8.16	509	−68	254	0.25
	Mean	19.1	8.2	1375	−49.3	686	0.7
	St. Dev.	2.9	0.4	2436	152.1	1219	1.4
HU3	Min	11.7	5.9	228	−106	114	0.1
	Max	19.7	8.9	4049	303	2025	2.2
	Median	17.55	7.04	509.5	153.6	254.5	0.25
	Mean	17.3	7.1	736	126	371	0.4
	St. Dev.	1.4	0.5	595	112	298	0.3

.

The water level measurements, carried out on the 88 selected wells representative of HU3, were elaborated together with a digital elevation model (DEM) and well geometrical features to obtain the piezometric head (m a.s.l.) of each sampling point.

In summary, the new investigations included: (1) hydrostratigraphic data acquisition; (2) water level measurements in 88 selected wells representative of the main phreatic aquifer; (3) 12 slug tests; (4) and physicochemical measurements carried out on 101 water samples.

Piezometric head, hydrostratigraphic, hydraulic and physicochemical data were processed using interpolation techniques to obtain a continuous spatial reconstruction. The method of spatial interpolation (Inverse Distance Weighting: IDW) was applied to produce a Digital Surface Model, and thickness, iso-concentration and iso-level maps. The IDW method estimates values of cells by weighting values (points) of geometric data in the proximity of each processed area (cell). The points closer to the cell centre weigh more in the process of weighting.

4. Results and Discussion

4.1. Lithology and Geometry of the Hydro-Stratigraphic Units

According to the geological data, the sedimentary sequences were divided into three hydrogeological units (HU1 to HU3, from bottom to top), derived from boreholes and in situ surveys. HU1 represents the deep aquifer hosted in Late Miocene siliciclastic deposits and limestones, which unconformably cover the igneous-metamorphic bedrock (mainly consisting of granitoids). The Late Miocene succession, mainly cropping out along the northern (Monte Poro massif) and southern (Palmi high) GTP boundaries, includes lithified arenites and conglomerates, and limestones (microbialites), all characterised by fracture permeability. The HU1 is lowered by the faults surrounding the plain, its top reaching a maximum depth of about 0.8 s of TWT (two-way travel time), as shown by exploration seismic profiles [22].

HU2 (aquiclude) comprises the Pliocene clayey and silty deposits of the Trubi Formation. The thickness of HU2 roughly ranges from 100 to 200 m. This unit crops out along the GTP boundaries and is lowered by the faults in the middle sector. The Digital Surface Model of the HU2 top (Figure 3a) highlights a complex pattern. In the northern sector, the Digital Surface Model (Figure 3a), and the AA’ and BB’ transversal profiles (Figure 4) show a depth of HU2 top between −60 and −40 m a.s.l. in the first 2 km eastward from the coastline, while it uplifts up to −20 m a.s.l. in the following 2 km. Proceeding eastward, the HU2 top reaches a maximum depth between −60 and −80 m a.s.l. Beyond this, it progressively uplifts up to 80 m a.s.l. in correspondence with the eastern GTP boundary. In the southern sector, the HU2 top is characterised by a deepening from the coastal to the middle sector (from −20 to −80 m a.s.l. in CC’ profile) and a subsequent uplift toward est. In the DD’, EE’ and FF’ longitudinal profiles (Figure 4), the deepening of the HU2 top (up to −80 m a.s.l.) can be recognised between the northern and southern GTP boundaries. The EE’ profile clearly shows a slope break of the HU2 top profile between 6 and 7 km from the north. At the whole plain scale, the Digital Surface Model of the HU2 top shows a main N-S striking depozone in the middle sector of the GTP.

The HU2 top deepening, from east to west, can be related to the CF and SLF, NE-SW striking, and WNW-dipping normal fault [22,23]. Furthermore, in our reconstruction, it is possible to observe an uplift of the HU2 top located between the NE-SW striking PGF and the parallel fault, respectively, with regard to the WNW- and ESE-dipping between the coastline and middle sector of the GTP. The presence of this structural high in the GTP, in continuity with the Palmi high, is in agreement with bibliographic sources [22,23,26,30]. The HU2 top deepening increases eastward again from PGF. In addition, in the southern GTP sector, the transversal profile CC’ (Figure 4) displays eastward tilting, which matches the geological cross-section proposed by the Reference [23]. Along the south of the GTP boundary, the Digital Surface Model shows a deepening of the HU2 top across PLFZ, a NW-SE striking and NE-dipping fault zone [22], while in the southwestern sector, the HU2 crops out on the Palmi High. According to our reconstruction, along the northern GTP boundary, the HU2 top deepening occurs in correspondence with the NGFZ NW-SE striking and SW-dipping main fault. Moreover, along a secondary fault (NW-SE striking and SW-dipping) of NGFZ, crossed by the Mesima River mouth [22], the Digital Surface Model and the longitudinal profiles (in particular, EE’) (Figure 4) record an additional deepening of the HU2 top toward the south.

The shallow and unconfined aquifer (HU3) groups the deposits from the Late Pliocene to Holocene. The HU3 bottom consists of the Calcareniti di Vinco (Late Pliocene-Pleistocene), including cross-bedded sands and calcarenites up to 70 m thick and locally cropping out along the GTP boundary. Upward, in the inner GTP sector, the HU3 includes the alluvial fan deposits (up to 30–35 m thick) of the Sintema di Taurianova (Late Pleistocene) made by clast-supported conglomerates, sands, and gravelly, silty and clayey sands. Instead, in the central portion of GTP, on top of the Calcareniti di Vinco, the HU3 continues with terraced marine deposits mainly composed of sand and silty sand with locally clayey silt intercalations (Middle-Upper Pleistocene) and terraced alluvial fan gravels (Late Pleistocene–Holocene). Finally, the HU3 top includes different Holocene deposits: coarse-grained alluvial sediments along the main rivers and streams, coarse aeolian sands of the two inland dune fields (located southward from the Nicotera Marina and Rosarno villages, respectively), modern coastal sediments, and detrital and gravitative deposits. Locally, in the sector northward of the Mesima mouth, the Holocene deposits are characterised by the presence of peat (recorded by boreholes, Figure 1).

The isopach map shows that HU3 (Figure 3b) is as thick as >100 m in a N-S striking area in the central part of GTP, while close to Taurianova, it reaches about 260 m. The HU3 depth is ca. 80 m around Polistena and close to Gioia Tauro’s harbour area, which is in the western-middle and eastern-middle sectors of the GTP, respectively.

We also investigated the spatial pattern of the HU3 thickness (Figure 5b), which was >100 m in the SE and middle sectors of GTP in correspondence to the depozone related to the so-called Gioia Tauro Extensional Fault System [22]. The maximum thickness, up to 260 m, was recorded on the hanging wall of CF between Taurianova and Polistena (Figure 1) where a negative Bouguer anomaly (−10 mGal) was shown by the Ref. [42]. Toward the coastline, high values (up to 100 m) were determined close to the Gioia Tauro harbour. In contrast, low values were observed in the surroundings of Gioia Tauro (20–40 m) and in the NW GTP sector (20–60 m) on the footwall of the NW-SE striking fault crossed by the Mesima River mouth.

In the area between the mouth of the Petrace River and Scinà (inset an in Figure 1), the hydrostratigraphic pattern differs from the rest of GTP, being characterised by HU3 and lying directly on the igneous-metamorphic bedrock. In this area, HU3 is hosted in fluvio-deltaic (mainly coarse-grained bedload) sediments at the mouth of the Petrace River and coastal deposits to the south. The Digital Surface Model of the bedrock top (Figure 5a) shows progressive deepening, up to about −65 m a.s.l., toward NW. In the isopach map (Figure 5b), the HU3 thickens up to 70 m along a prevalent SE-NW direction. In this area, the deepening toward NW of contact between the HU3 bottom and the igneous-metamorphic bedrock is clearly conditioned by NE-SW striking and WNW-dipping PGF. To summarise, our spatial reconstruction of the HU2 top depth shows a good correlation with the GTP structural setting.

4.2. Hydrography and Groundwater Surface

The GTP hydrograph pattern shows a relatively simple structure characterised by three main sub-basins: (i) The Mesima basin located in the northern sector with a SW, W, and subordinately NW drainage direction; (ii) the Petrace basin, in the southern sector, showing a N and NW drainage direction; and (iii) the Budello basin positioned in the central portion showing W-NW drainage [26]. The three sub-basins show significant differences in terms of flow rates and morphological features. The central basin is represented almost exclusively by the Budello River, a torrential watercourse only fed during winter (November to March; [26]). The small size of the Budello River is due to the limited extent of its catchment. In fact, the watercourses draining the Aspromonte Massif, which should feed the Budello River, deviate towards NW and SW to the Mesima and Petrace basins, respectively. This evidence explains the absence of drainage in the central area whose extension is >50 km². The Petrace sub-basin shows enrichment in branched waterways that become sub-parallels downstream. The waterways engrave the sediments converging then towards the Petrace River that cuts the northern portion of the Palmi High in the coast’s proximity. Finally, the Mesima Basin is characterised by the interaction with the medium-course drainage of the Mesima River, which flows for about 20 km towards SSW before reaching the coast. The Mesima River receives, close to Rosarno, the contribution of the Metramo River that flows behind (east side) the Budello basin (Figure 1).

The groundwater flow direction coincided with the hydrography, as shown in the iso-piezometric map constructed with the new data (Figure 6a). The well locations, final piezometric contour line map (blue lines), and relative flow directions are plotted in Figure 6a. The map shows a NW-SE-oriented prevalent flow direction along the eastern sector of GTP in the proximity of the recharge area (Aspromonte massif). During its migration towards the coast (NW), the flow separates into two main sub-flows that are W-E- and NW-SE-oriented, respectively, due to the presence of a wide watershed in the southern part of Rosarno. The W-E-oriented sub-flow feeds both the Budello and Petrace Basins, whereas the NW-SE-oriented one migrates northwards, feeding the Mesima Basin. By comparing the iso-piezometric maps reported in Figure 6a (new 2021 vs. 1970s [36] piezometric heads), the presence of comparable main flow directions can be highlighted.

As reported in Figure 6b, a substantial decrease in the overall piezometric level (up to −35 m from the 1970s to 2021) was observed in the middle sector of the plain among the urban centres of Rosarno, Rizziconi, Taurianova, Cittavova and Polistena. The main rise (up to 25 m from the 1970s to 2021) was recorded in correspondence with the Petrace River. A substantial unchanged level apparently characterises the coastal sector between Gioia Tauro and Nicotera Marina. The issue could be linked to, in addition to a general increase of the annual average temperatures, the continuous transformations in the last decades of GTP in terms of land use (e.g., change in crops: from olive and citrus groves to kiwi one) and urban growth, which produced an increase of groundwater overexploitation. Furthermore, the central portion of the plain, as reported above, is characterised by the presence of a wide watershed where the main flow undergoes the separation in two main sub-flows, respectively oriented as W-E and NW-SE. This asset could represent a further explanation of recorded decreases.

4.3. Hydraulic Properties of the Hydrogeological Aquifer Units

The hydraulic conductivity (k) of HU3 was calculated using 225 pumping test data. These data were validated by means of the values obtained through the slug test survey (Table 1). The k inferred by slug tests were mainly (9 out of 12) characterised by values of the order of magnitude of 10⁻⁴ m/s, two (ST7 and 11) of 10⁻³ m/s, and one (ST4) of 10⁻⁵. This variability is likely due to the different deposits hosting the upper part of HU3. The k-values of 10⁻⁴ m/s (ST1, 2, 3, 5, 6, 8, 9, 10, 12) were computed toward the coastline, where the top of HU3 consists of coastal (marine and aeolian) medium sands and alluvial (mainly sand and silty sand) deposits. On the other hand, the highest value of 10⁻³ m/s pertained to the inland dune fields, consisting of coarse-grained sands.

The lowest k-value of 10⁻⁵ characterises the central portion of GTP, close to Rizziconi (Figure 1), where the top of HU3 includes terraced marine deposits. The reliability of the k-values obtained in the present work is supported when compared with those reported by the Ref. [43] for the HU3 lithologies (Figure 7).

Considering this good correlation, we used the 225 empirical values to analyse the HU3 k spatial pattern at the GTP scale (Figure 8). The k-values are about 1 × 10⁻⁴ in the middle and coastal sectors, with local areas characterised by higher values (10⁻³), mainly found close to the inland dune fields. The lowest k-values (≤10⁻⁵) are observed in the eastern and northern GTP sectors.

4.4. Groundwater Characterisation

In the study area, the water temperature varied from 14 to 25.8 °C. The pH and Electrical Conductivity (EC) values showed great variability. The former was comprised between 5.9 and 8.9, while EC ranged, according to the classification of [44], from very fresh (238 μS/cm) to very saline (10,520 μS/cm). The Eh values spanned from −253 mV to 303 mV (Table 2 and Figure 9).

From a general point of view, the main anomalies were recognised along the northern boundary of the Plain, in the proximity of one of the main fault zones (NGFZ). In this sector, the distribution maps of Figure 9 highlight a good correlation among temperature, pH, and EC (Figure 9). Furthermore, the northern sector recorded the most negative Eh values.

The southern portion of the Plain showed comparable values, with the central sector characterised by low EC (<1000 μS/cm) and positive Eh values and relatively low temperatures (<18 °C). Only one restricted area, along the coastline and confined by the Palmi high (at Scinà village, Figure 1), revealed anomalous EC values (>1500 μS/cm not associated with anomalous temperatures) for the HU3. In this sector, characterised by a peculiar hydrostratigraphic framework, the presence of a restricted recharge area and numerous water wells could favour seawater intrusion.

Analysing the two main aquifers separately (HU3 and HU1), they had slightly comparable spatial distribution in terms of physicochemical parameters but were probably attributable to different processes. In detail, the deep aquifer (HU1), whose features were provided by artesian wells (18 wells; Figure 9), highlighted the greater spatial heterogeneity. Waters hosted in HU1 showed the highest temperature (up to 25.8 °C), EC (mean value: 1375 μS/cm, up to 10,520 μS/cm), pH (>8), and negative Eh close to the north-eastern portion of the Plain that allows for the hypothesis of a dominant role of faults to the fluid’s circulation and evolution. These anomalies were not registered in the remaining artesian wells located in distal areas from the main faults. This evidence agrees with data provided by the ref. [19] who observed that the warmest and most saline waters (25 < T < 38 °C; TDS > 0.5 g/L, sometimes reaching 10–20 g/L) discharged preferentially along the WNW-ESE belt running from Nicotera to Galatro (northern sector of the Plain). The sector located near the Palmi High also showed anomalous, though lower, temperature and salinity values (maximum value of 1900 μS/cm and 22 °C, respectively). In this sector, the tectonic structure, represented by the PLFZ, is likely playing a decisive role in the circulation, with a mixing of deep and surface rates. Although the faults’ dominant role in circulation and chemical evolution of the deep systems is clear, the main processes that favour the formation of waters characterised by high temperatures and salinities still remain uncertain. The interaction with evaporitic successions, not present in the studied area, can be excluded, as well as a possible seawater intrusion considering the great distance of the wells from the coastline. Possible explanations are found in the presence of both “connate” waters and deep mature waters deriving from a long circulation path in the basal crystalline units, which likely originates high-salinity and NaCl-rich waters [45,46]. The well-developed fault systems could play an important role in the final circulation and rise to the surface of this type of water (e.g., [47]). Further studies should be carried out to provide more information about the ongoing geochemical processes.

As observed for HU1, in the northern portion of the plain (northward of the mouth of Mesima River), waters representative of the shallow HU3 showed relatively high EC (>1000 μS/cm) and slightly negative Eh values. The northern portion of the Plain is characterised by the presence of peat in the Holocene deposits that can favour anoxic conditions in groundwater. Furthermore, the complex structural setting located between the main faults of the NGFZ (the regional tectonic system along which is located the Galatro thermal area [48], green star in Figure 1) can favour a possible exchange between the two main aquifers [49,50]. Finally, the middle and inner portions of GTP (around Rizziconi, Taurianova and Polistena: Figure 1), as also observed for the deep HU1 aquifer, showed the lowest temperature and EC values.

5. Concluding Remarks

This study provided a detailed characterisation of the GTP groundwater’s hydrostratigraphic framework, hydraulic properties, and physicochemical status. The multidisciplinary approach allowed to provide a detailed reconstruction of the main hydrogeological asset of GTP consisting of shallow (HU3) and deep (HU1) aquifers separated by an aquitard (HU2) (Figure 10).

Furthermore, the direct influence of the structural setting on the HU3 thickness spatial pattern was observed. For the HU3, the hydraulic conductivity was calculated by an empirical method ranging from 10⁻⁴ and 10⁻⁵ m/s with peaks of 10⁻³ m/s in correspondence of inland Holocene aeolian deposits. The water level measurements allowed for identification of the main flow direction and estimation of the piezometric level variation from the 1970s to 2021, showing a general downwards trend with a maximum of up to −35 m.

The groundwater physicochemical parameters showed anomalous high values of temperature and electrical conductivity (up to 25.8 °C and 10,520 µS/cm, respectively) in the deep aquifer (HU1) along the main tectonic lineaments, which likely played a key role in fluid circulation and evolution (Figure 10). Conversely, shallow groundwater showed lower electrical conductivity (mean value of 736 µS/cm), except at Scinà village (also characterised by a peculiar hydrostratigraphic framework) where values > 1500 µS/cm were recorded and linked to the seawater intrusion.

The results achieved in the present study can be regarded as essential guidelines for optimum management of groundwater protection in the study area. Furthermore, the obtained results (representing a starting point) make the GTP groundwater system an example of interaction between shallow and deep waters, with an important role played by regional faults and the interaction of fault systems.

To improve the hydrostratigraphic reconstruction, passive seismic surveys and VES surveys will be realised where few borehole data and aquitard top depth data are available. In addition, new stratigraphic data are expected to help analyse the lateral variation of the shallow aquifer due to the heteropic facies typical of the alluvial environment. Major, minor and trace elements and water isotopes can also provide additional information for defining the main geochemical processes acting in the circulating waters and understanding the deep relationship between tectonic structures and geochemical features.