Evidence for Groundwater Salinity Origin Based on Hydrogeochemical and Isotopic ( 2 H, 18 O, 37 Cl, 3 H, 13 C, 14 C) Approaches: Sousse, Eastern Tunisia

: The key processes responsible for the rise in groundwater salinization in the Mio–Pliocene aquifer system of Sousse (Tunisia, eastern coastline) were identiﬁed through a multidisciplinary approach based on the use of geochemical, stable ( 2 H, 13 C, 18 O and 37 Cl) and radioactive ( 3 H and 14 C) isotope methods. In the study region, the mineralization of groundwaters is related to water–rock interaction ascribed to the dissolution of minerals in evaporite rocks, as well as to saltwater intrusion. Both processes explain the development of groundwaters in which Cl and Na dominantly determine the groundwater quality deterioration state. The isotopic and geochemical signatures of the studied groundwaters are clearly explained by the (i) occurrence of saline basins (sebkhas adjacent to the study region), (ii) type of rocks found below the ground surface, and (iii) cation exchange between clays and groundwaters.


Introduction
Salinization of water resources is one of the most prevalent processes that reduces water quality and threatens further water abstraction. This issue is worsened in coastal aquifers, particularly in dry and semi-arid areas where freshwater is threatened by seawater intrusion. The cause of rising groundwater salinity may be natural or anthropogenic, such as water-rock interaction (WRI), which involves the direct dissolving of evaporites by fresh groundwater, concentration of dissolved salt by evaporation, and infiltration of saline surface water.
Additionally, human activities can accelerate water quality deterioration because of pollution and sources of saline solutions (industrial liquid or solid waste, salt mining activities, or agricultural backflow) [1].
Once different salinization mechanisms are identified, it is possible to retrace the origin of groundwater, follow its evolution, and predict its future change [2].
The 1300 km of Tunisian coastline has a coastal width of 20 to 60 km and a total surface area of around 40,000 km 2 [3].
Along the northeast coastal zone of the Tunisian Sahel, the preservation of water quantity and quality are major concerns.
Environmental isotope tracers and age indicators have been used in arid regions, such as North Africa, for assessing groundwater systems, particularly in terms of residence time

Main Features of the Study Area
The Oued Laya aquifer is situated in the Tunisian Sahel about 140 km south of the city of Tunis, and is close to the city of Sousse, one of Tunisia's most important seaside resort towns ( Figure 1). It is near the coast at the end of the littoral zone of the Gulf of Hammamet in the southwestern part of the Mediterranean Sea. The landscape is a coastal plain that gradually slopes (2%) in the direction of the coast. It extends 15 km from north to south between Zaouia Sousse and Kalaa Kebira, and 25 km in width from the sea to the south west. At its highest point the altitude is about 100 m [23]. The topography of the Sousse region is composed of plains and hills, including the anticline of Kalaa Kebira and Mouredine. The plain is low and flat and has a sandy coastline along the Mediterranean Sea [23].
The region of Sousse is characterized by a semi-arid climate. The annual average rainfall is 350 mm with significant variations; between 1901 and 2020, the extremes recorded were 110 and 700 mm in 2004 and 1969, respectively. This climate prevails in the Mediterranean regions with a mean temperature of 19 °C. September through February are the wettest months, while July and August are the driest [24,25].
The evapo-transpiration (ETP) is estimated to be higher than 1750 mm/year by the Riou approach [28]. The region is economically characterized by important agriculture, commerce, and tourist activity. The number of inhabitants increases considerably during the summer period, with a population approaching 300,000.
The groundwater in the Oued Laya aquifer is divided into two distinct aquifers. The first is a shallow aquifer that extends to about 60 m in depth and has a reservoir primarily composed by Mio-Pliocene sediments (with a few gypsum lenses scattered throughout the geological formations), covering a surface of about 217 km 2 and varying in thickness from 10 to 40 m. The second is a deep aquifer that lies between 90 and 250 m in the sand-

Main Features of the Study Area
The Oued Laya aquifer is situated in the Tunisian Sahel about 140 km south of the city of Tunis, and is close to the city of Sousse, one of Tunisia's most important seaside resort towns (Figure 1). It is near the coast at the end of the littoral zone of the Gulf of Hammamet in the southwestern part of the Mediterranean Sea. The landscape is a coastal plain that gradually slopes (2%) in the direction of the coast. It extends 15 km from north to south between Zaouia Sousse and Kalaa Kebira, and 25 km in width from the sea to the south west. At its highest point the altitude is about 100 m [23]. The topography of the Sousse region is composed of plains and hills, including the anticline of Kalaa Kebira and Mouredine. The plain is low and flat and has a sandy coastline along the Mediterranean Sea [23].
The region of Sousse is characterized by a semi-arid climate. The annual average rainfall is 350 mm with significant variations; between 1901 and 2020, the extremes recorded were 110 and 700 mm in 2004 and 1969, respectively. This climate prevails in the Mediterranean regions with a mean temperature of 19 • C. September through February are the wettest months, while July and August are the driest [24,25].
The evapo-transpiration (ETP) is estimated to be higher than 1750 mm/year by the Riou approach [28]. The region is economically characterized by important agriculture, commerce, and tourist activity. The number of inhabitants increases considerably during the summer period, with a population approaching 300,000.
The groundwater in the Oued Laya aquifer is divided into two distinct aquifers. The first is a shallow aquifer that extends to about 60 m in depth and has a reservoir primarily composed by Mio-Pliocene sediments (with a few gypsum lenses scattered throughout the geological formations), covering a surface of about 217 km 2 and varying in thickness from 10 to 40 m. The second is a deep aquifer that lies between 90 and 250 m in the sandstone deposits of the Miocene. In this area, there are no perennial streams due to the climatic specificity of the region (semi-arid) and geological features. During severe flooding and intense storms, a large number of minor streams (Oueds) may transport water. The main Oued that reaches the Tunisian city of Sousse's shore is Oued Laya [23][24][25].
According to official yearly data from 2020, the region's renewable water resources are predicted to be 3.3 × 10 6 m 3 /year with an annual abstraction of about 3.133 × 10 6 m 3 . The number of wells has almost doubled from 400 to 795 wells between 1980 and 2010. [29].
The agricultural, touristic, and industrial sectors place the most significant demands on the Oued Laya aquifer. Groundwater abstraction from the aquifer is mainly through a relatively large number of shallow wells and boreholes with a depth ranging from 10 to 60 m [29].

Geological and Hydrogeological Setting
In eastern Tunisia, the north-south axis separates the Atlasic area from the Sahel domain. [30,31]. According to the geological map, and due to the architectural characteristic of this domain sealed by Mio-Plio-Quaternary deposits, the research area looks like a weak tectonic area and is less affected on the surface ( Figure 2). matic specificity of the region (semi-arid) and geological features. During severe flooding and intense storms, a large number of minor streams (Oueds) may transport water. The main Oued that reaches the Tunisian city of Sousse's shore is Oued Laya [23][24][25].
According to official yearly data from 2020, the region's renewable water resources are predicted to be 3.3 × 10 6 m 3 /year with an annual abstraction of about 3.133 × 10 6 m 3 . The number of wells has almost doubled from 400 to 795 wells between 1980 and 2010. [29].
The agricultural, touristic, and industrial sectors place the most significant demands on the Oued Laya aquifer. Groundwater abstraction from the aquifer is mainly through a relatively large number of shallow wells and boreholes with a depth ranging from 10 to 60 m [29].

Geological and Hydrogeological Setting
In eastern Tunisia, the north-south axis separates the Atlasic area from the Sahel domain. [30,31]. According to the geological map, and due to the architectural characteristic of this domain sealed by Mio-Plio-Quaternary deposits, the research area looks like a weak tectonic area and is less affected on the surface (Figure 2).
The outcropping formations in the Oued Laya watershed include: the Quaternary, composed by aQ, dQ, cQ, and Qv; the Pliocene (Pl), located in the region of Akouda composed by marl and sandstone; and the Mio-Pliocene, in the Kalaa Kebira anticline formed by conglomerates, sand, and clay. However, there is significant geologic complexity in the subsurface linked to Triassic halokinesis along a major fault system and severe volcanic activity from about the Aptian to Late Cretaceous [30].
Along with all the significant geodynamic events of the Atlas Mountains, almost 8000 m of post-Permian deposits have been recorded in the basin, indicating a weak tectonic The outcropping formations in the Oued Laya watershed include: the Quaternary, composed by aQ, dQ, cQ, and Qv; the Pliocene (Pl), located in the region of Akouda composed by marl and sandstone; and the Mio-Pliocene, in the Kalaa Kebira anticline formed by conglomerates, sand, and clay.
However, there is significant geologic complexity in the subsurface linked to Triassic halokinesis along a major fault system and severe volcanic activity from about the Aptian to Late Cretaceous [30].
Along with all the significant geodynamic events of the Atlas Mountains, almost 8000 m of post-Permian deposits have been recorded in the basin, indicating a weak tectonic region sinking since the Jurassic. The El Hdadja fault and the Kairouan-Sousse fault serve as basin boundaries To the north and south, respectively [30].
The aquifer system of Oued Laya is composed by a shallow aquifer with Mio-Pliocene sandstone deposits in the Segui formation, which has a thickness varying between 10 and 40 m, and a sandy confined aquifer located in the Oum Douil formation, which has a thickness varying from 90 to 250 m [5].
A layer of clay about 30 to 40 m thick that is sufficiently continuous and impermeable to provide a hydraulic head for the deep aquifer separates the Miocene aquifer from the Mio-Pliocene aquifer ( Figure 3). This aquifer is located between two anticlines of Kalaa Kebira to the north and that of Moureddine to the south. The density of the wells is variable. It is very strong in areas near Oued Laya, Oued el Kharroub, and Oued El Kebir where the water is of fairly good quality. region sinking since the Jurassic. The El Hdadja fault and the Kairouan-Sousse fault serve as basin boundaries To the north and south, respectively [30]. The research area's geological formations range from the Cretaceous to the Quaternary [32][33][34].
The aquifer system of Oued Laya is composed by a shallow aquifer with Mio-Pliocene sandstone deposits in the Segui formation, which has a thickness varying between 10 and 40 m, and a sandy confined aquifer located in the Oum Douil formation, which has a thickness varying from 90 to 250 m [5].
A layer of clay about 30 to 40 m thick that is sufficiently continuous and impermeable to provide a hydraulic head for the deep aquifer separates the Miocene aquifer from the Mio-Pliocene aquifer ( Figure 3). This aquifer is located between two anticlines of Kalaa Kebira to the north and that of Moureddine to the south. The density of the wells is variable. It is very strong in areas near Oued Laya, Oued el Kharroub, and Oued El Kebir where the water is of fairly good quality.
According the results of pumping tests performed in boreholes, the transmissivity of the aquifer varies from 3 × 10 −6 m 2 /s to 1.8 × 10 −4 m 2 /s [32]. The shallow aquifer is recharged by rain and runoff. The Mediterranean Sea constitutes its natural outlets as shown in the piezometric map ( Figure 4). In the Miocene deep aquifer, water also flows from SW to NE. The shallow groundwater near the coastline has a higher salinity. Previous hydrogeological studies [5,[32][33][34][35] indicate that there is no hydraulic connection between the shallow unconfined and the confined aquifer. According the results of pumping tests performed in boreholes, the transmissivity of the aquifer varies from 3 × 10 −6 m 2 /s to 1.8 × 10 −4 m 2 /s [32].

SW-NE cross section
The shallow aquifer is recharged by rain and runoff. The Mediterranean Sea constitutes its natural outlets as shown in the piezometric map ( Figure 4). In the Miocene deep aquifer, water also flows from SW to NE. The shallow groundwater near the coastline has a higher salinity. Previous hydrogeological studies [5,[32][33][34][35] indicate that there is no hydraulic connection between the shallow unconfined and the confined aquifer.

Material and Methods
The geographical position and the groundwater table below ground level were recorded at 37 sampling points in the study area ( Figure 1). The pH and electric conductivity were measured at each sample location with a WTW multi-Parameter field meter (Accuracy ±1% of reading or ±1 µS/cm whichever is greater) which was calibrated in the laboratory with 3 standard solutions (84 µS/cm, 1413 µS/cm and 12.88 mS/cm).
Field campaigns were carried out between 2009 and 2013 (Table 1) and the groundwater sampling procedures followed those recommended by the Isotope Hydrology Section of the International Atomic Energy Agency (IAEA) [36].
The chemical composition of the groundwater does not fluctuate remarkably during the different seasons (dry/wet). Samples were collected in order to measure major elements and stable and radioactive isotopes (δ 2 H, δ 18 O, δ 13 C, δ 37 Cl, 3 H and 14 C) ( Table 1).
The water samples taken from the shallow aquifer (boreholes and dug wells) ( Chemical analyses of the water samples were performed at the Isotope Hydrology and Geochemistry Unit at the Centre National des Sciences et Technologies Nucléaires, Tunis, Tunisia (CNSTN). Ca, Mg, Na, and K (the major cations) were determined from filtered samples by atomic absorption spectrometry, and Cl, SO4, and NO3 (the major anions) were determined in filtered samples by ion chromatography on a Dionex DX 120 equipped with AG14 and AS14 Ion Pac columns and an AS-40 auto-sampler. In most samples the electric charge balance between major anions and major cations is better than ±5%.

Material and Methods
The geographical position and the groundwater table below ground level were recorded at 37 sampling points in the study area ( Figure 1). The pH and electric conductivity were measured at each sample location with a WTW multi-Parameter field meter (Accuracy ±1% of reading or ±1 µS/cm whichever is greater) which was calibrated in the laboratory with 3 standard solutions (84 µS/cm, 1413 µS/cm and 12.88 mS/cm).
Field campaigns were carried out between 2009 and 2013 (Table 1) and the groundwater sampling procedures followed those recommended by the Isotope Hydrology Section of the International Atomic Energy Agency (IAEA) [36].
The chemical composition of the groundwater does not fluctuate remarkably during the different seasons (dry/wet). Samples were collected in order to measure major elements and stable and radioactive isotopes (δ 2 H, δ 18 O, δ 13 C, δ 37 Cl, 3 H and 14 C) ( Table 1).
The water samples taken from the shallow aquifer (boreholes and dug wells) ( Chemical analyses of the water samples were performed at the Isotope Hydrology and Geochemistry Unit at the Centre National des Sciences et Technologies Nucléaires, Tunis, Tunisia (CNSTN). Ca, Mg, Na, and K (the major cations) were determined from filtered samples by atomic absorption spectrometry, and Cl, SO 4 , and NO 3 (the major anions) were determined in filtered samples by ion chromatography on a Dionex DX 120 equipped with AG14 and AS14 Ion Pac columns and an AS-40 auto-sampler. In most samples the electric charge balance between major anions and major cations is better than ±5%. Spectrometry (CF-IRMS) coupled to a gas bench. Oxygen isotope ratios were determined on CO 2 gas that was formed after equilibration with CO 2 with a known oxygen isotope composition for 24 h at 25 • C [37]. Hydrogen isotope ratios were determined on H 2 gas that was formed after reduction of the water sample with Cr at 850 • C [38]. δ 18 O and δ 2 H were measured relative to internal water standards that were calibrated using IAEA Vienna Standard Mean Ocean Water (VSMOW). Data were normalized following [39] and are reported relative to V-SMOW. Samples were measured at least in duplicate, and the analytical precision is ±0.1‰ for δ 18 O and ±1‰ for δ 2 H. The data are reported as δ 18 Tritium ratios ( 3 H) were determined at the C 2 TN-former ITN, Lisbon, Portugal. Tritium content was measured after electrolytic enrichment with liquid scintillation spectrometry using a Liquid Scintillation Counter (LSC, PACKARD Tri-Carb 2000 CA/LL) [40,41]. Tritium ratios are reported in Tritium Units (TU). One TU is defined by the isotope ratio 3 H/ 1 H = 10 −18 . δ 13 C and 14 C were measured for 23 samples. The analyses were carried out on the Total Dissolved Inorganic Carbon (TDIC) of the groundwater. This was collected in the field through precipitation of BaCO 3 at a pH above 8.3. 14 C was determined at the Isotope Hydrology and Geochemistry Unit of the CNSTN, Tunis, Tunisia. It was determined after synthesis of benzene from the BaCO 3 and subsequent liquid scintillation spectrometry [42] using a PACKARD Tri-Carb 2750 TR/LL. 14 C concentrations are expressed as percentage of modern carbon (pmC). Corrected ages of the water samples were defined following the IAEA models and are expressed in years. The δ 13 C values are reported in ‰ vs. V-PDB, with accuracy of ±0.1%, and were analyzed by isotope ratio mass spectrometry at the C2TN-former ITN, Lisbon, Portugal. δ 37 Cl was measured to Standard Mean Oceanic Chloride (SMOC) at Utrecht University (Department of Geochemistry) according to the method that was developed by Eggenkamp [20]. δ 37 Cl is determined on methyl chloride (CH 3 Cl) gas by Isotope Ratio Mass Spectrometry (IRMS). In this case, the ratio between the masses 50 and 52 is determined. CH 3 Cl was produced using the following method: KNO 3 was added to a water sample, as ionic strength adjustment buffer, together with citric acid and Na 2 HPO 4 to obtain buffered low and stable pH. This high ionic strength is important to obtain small AgCl crystals and a low pH is important to avoid the formation of (black) Ag 2 O after the addition of silver ions (in the form of AgNO 3 ) that will prohibit the formation of CH 3 Cl. The buffered solution was heated to a temperature of about 80 • C. Then, a solution of 0.2 M AgNO 3 was added so that AgCl precipitated from the mixture. This AgCl suspension was filtered over a Whatman GF/F glass microfiber filter. This filter was dried at 80 • C for one night and the next day the filter was put into a glass reaction tube. Above the filter a restriction was made. Then, methyl iodide (CH 3 I) was added to the reaction tube. Thereafter, the tube was attached to a vacuum line. The CH 3 I was frozen on the filter with liquid N 2 and the reaction tube was evacuated. The reaction tube was then closed at the restriction resulting in an ampoule that only contains CH 3 I and AgCl. The ampoule was placed, for two days, in an oven at 80 • C. After this time, all AgCl will have been reacted to CH 3 Cl. Then, the ampoule was put into a break seal connected to a gas chromatograph. With the use of a helium flow, the CH 3 Cl was separated from the remaining CH 3 I over a Porapak-Q 80-100 mesh-packed column. Pure CH 3 Cl was then collected in a gas tube and its Cl isotope ratio was measured on a VG ISOGAS SIRA 24 stable isotope ratio mass spectrometer. The chlorine isotope ratio is reported versus a sample of ocean water that was taken near Madeira (Portugal). The chlorine isotope composition of this ocean water sample was proven to be equal to the chlorine isotope composition of average ocean water and for that reason it could be taken as representative for Standard Mean Oceanic Chloride (SMOC) [43]. The precision (2σ) of chlorine isotope measurements was ±0.1‰.

Results and Discussion
Between 2009 and 2013, fieldwork campaigns were undertaken at 37 sampling sites (including boreholes and dug wells), in order to improve knowledge on the (i) source of groundwater salinity, (ii) hydrodynamics, and (iii) flow patterns. The shallow groundwater (surficial aquifer) is characterized by a dry residue mean value of 4420 mg/L. In this shallow unconfined aquifer system, the electrical conductivity (EC) data were measured during the sampling campaign. They show high values increasing from (1.38 mS/cm) to (12.7 mS/cm). The highest EC values are ascribed to sampling sites located near the central part of the study region, whereas the lower EC values are mostly found away from the central axis that follows the Oued ( Figure 5). On the other hand, in the deeper (confined) aquifer systems groundwater mineralization rises, achieving salinity values that make such groundwaters unsuitable for drinking, agriculture, or industrial uses (average EC values varying from 7.7 to 9.2 mS/cm).
The presentation of geochemical data in the form of graphical charts such as a Piper diagram [44] is helpful for recognizing and differentiating various water facies. The major ions of the water analyses are depicted in the Piper diagram ( Figure 6). The chemical composition of groundwater plotted in a Piper trilinear equivalence diagram ( Figure 6) shows a trend of water classified under a Na Cl SO 4 Ca water type. The presentation of geochemical data in the form of graphical charts such as a Piper diagram [44] is helpful for recognizing and differentiating various water facies. The major ions of the water analyses are depicted in the Piper diagram ( Figure 6). The chemical composition of groundwater plotted in a Piper trilinear equivalence diagram ( Figure 6) shows a trend of water classified under a Na Cl SO4 Ca water type.   The presentation of geochemical data in the form of graphical charts such as a diagram [44] is helpful for recognizing and differentiating various water facies. The ions of the water analyses are depicted in the Piper diagram ( Figure 6). The chemical position of groundwater plotted in a Piper trilinear equivalence diagram ( Figure 6) s a trend of water classified under a Na Cl SO4 Ca water type.

Na + /Cl − Correlation
In coastal areas, a strong correlation between Na-Cl content in groundwater systems is often observed, which can be due to first air masses condensation (water vapor), sometimes associated with the dissolution of marine aerosols deposited in the soil in coastal regions during recharge [45], or when seawater intrusion mechanisms occur.
In the Mio-Pliocene groundwater samples, a similar positive relation was found, i.e., the correlation coefficient Na/Cl was equal to 0.88 (n = 37). However, when the molar ratio between these two ions was calculated for the sampled water, the value obtained was almost 1, pointing to halite dissolution as the main source of Mio-Pliocene groundwater mineralization (Figure 7a). Assuming, by hypothesis, that seawater intrusion mechanism is the main process responsible for the salt content in the Mio-Pliocene fresh waters, the Na-Cl molar ratio should be smaller or equal to the 0.86 value of Mediterranean water composition [11].
Water 2023, 15, x FOR PEER REVIEW 10 of 23

Na + /Cl − Correlation
In coastal areas, a strong correlation between Na-Cl content in groundwater systems is often observed, which can be due to first air masses condensation (water vapor), sometimes associated with the dissolution of marine aerosols deposited in the soil in coastal regions during recharge [45], or when seawater intrusion mechanisms occur.
In the Mio-Pliocene groundwater samples, a similar positive relation was found, i.e., the correlation coefficient Na/Cl was equal to 0.88 (n = 37). However, when the molar ratio between these two ions was calculated for the sampled water, the value obtained was almost 1, pointing to halite dissolution as the main source of Mio-Pliocene groundwater mineralization (Figure 7a). Assuming, by hypothesis, that seawater intrusion mechanism is the main process responsible for the salt content in the Mio-Pliocene fresh waters, the Na-Cl molar ratio should be smaller or equal to the 0.86 value of Mediterranean water composition [11]. Considering that both hypotheses can occur in the research region, i.e., seawater intrusion and salts dissolution, the excess of sodium concentration must also be explained. In fact, the excess of Na is observed in some water samples, either with high or low mineralization, which can be ascribed to cation exchange mechanisms [26,46,47], often found in sedimentary basins, during seawater mixing processes, and also in deltaic regions with high amount of organic matter and clay minerals.

Br − /Cl − Relationship
Likewise, the Br/Cl ratio allows us to obtain unique information to identify possible seawater mixing processes within coastal aquifers, where the main origin of salts in the groundwater (seawater intrusion vs. evaporate dissolution) is unknown or doubtful. Due to the long residence time of the ocean water masses, combined with a relative stable Br-Cl ratio around 1.5 × 10 −3 [10], such a ratio can also be applied as a criterion in the characterization of the main origin of salt in groundwater systems. In our research area, the Br/Cl ratio of the sampled groundwater was plotted versus the Cl content (Figure 7b). Apart of a small group of samples, the graphic distribution does not point to a clear origin of the salts. This dispersion is most probably ascribed to the mixture of fresh water and seawater and/or to the dissolution of marine aerosols dissolved by precipitation during aquifer recharge.
Nonetheless, it is important to mention that most of the points plotted far from the Br/Cl seawater ratio line are away from the coastline and off of the most probable seawater intrusion area. Considering that both hypotheses can occur in the research region, i.e., seawater intrusion and salts dissolution, the excess of sodium concentration must also be explained. In fact, the excess of Na is observed in some water samples, either with high or low mineralization, which can be ascribed to cation exchange mechanisms [26,46,47], often found in sedimentary basins, during seawater mixing processes, and also in deltaic regions with high amount of organic matter and clay minerals.

Br − /Cl − Relationship
Likewise, the Br/Cl ratio allows us to obtain unique information to identify possible seawater mixing processes within coastal aquifers, where the main origin of salts in the groundwater (seawater intrusion vs. evaporate dissolution) is unknown or doubtful. Due to the long residence time of the ocean water masses, combined with a relative stable Br-Cl ratio around 1.5 × 10 −3 [10], such a ratio can also be applied as a criterion in the characterization of the main origin of salt in groundwater systems. In our research area, the Br/Cl ratio of the sampled groundwater was plotted versus the Cl content (Figure 7b). Apart of a small group of samples, the graphic distribution does not point to a clear origin of the salts. This dispersion is most probably ascribed to the mixture of fresh water and seawater and/or to the dissolution of marine aerosols dissolved by precipitation during aquifer recharge.
Nonetheless, it is important to mention that most of the points plotted far from the Br/Cl seawater ratio line are away from the coastline and off of the most probable seawater intrusion area.

Ca 2+ /Mg 2+ Relationship
The concentration of Ca 2+ and Mg 2+ in the shallow groundwaters varies widely, although in most situations, Ca 2+ predominates over Mg 2+ (Figure 8a). Groundwater samples with a Ca/Mg ratio greater than 1 may not have been mixed with saltwater (marine ratio = 0.2) [6], and mineralization from another source of Ca should be considered. Dispersed layers of gypsum exist in the subsurface geology of the reference area and these can be regarded as a possible "external" calcium supply for the aquifer. The Ca-Mg correlation observed in the groundwater samples reflects gypsum dissolution and indicates that the gypsum layers are a likely source of calcium for the system. Nevertheless, ion exchange within the aquifer system is also a possible mechanism in this system since it can explain the increase in sodium (Na/Ca = 2) in addition to the decrease in calcium in the groundwater content, as calcium may be exchanged on clay surfaces with sodium.

Ca 2+ /Mg 2+ Relationship
The concentration of Ca 2+ and Mg 2+ in the shallow groundwaters varies widely, although in most situations, Ca 2+ predominates over Mg 2+ (Figure 8a). Groundwater samples with a Ca/Mg ratio greater than 1 may not have been mixed with saltwater (marine ratio = 0.2) [6], and mineralization from another source of Ca should be considered. Dispersed layers of gypsum exist in the subsurface geology of the reference area and these can be regarded as a possible "external" calcium supply for the aquifer. The Ca-Mg correlation observed in the groundwater samples reflects gypsum dissolution and indicates that the gypsum layers are a likely source of calcium for the system. Nevertheless, ion exchange within the aquifer system is also a possible mechanism in this system since it can explain the increase in sodium (Na/Ca = 2) in addition to the decrease in calcium in the groundwater content, as calcium may be exchanged on clay surfaces with sodium.

Ca 2+ /SO4 2− Relationship
In line with the previous hypothesis, the reduction in Ca 2+ concentrations in comparison to SO4 2− in the more sulphate-rich sampled shallow groundwaters, as shown in (Figure 8b), is possibly the result of cation exchange reactions attributed to the presence of a fraction of clay and organic matter within the aquifer matrix as Na + is released, here resulting in calcium concentration as is present here. If we consider the gypsum lenses as an additional potential Ca 2+ source for groundwater mineralization, Ca 2+ might be removed by Na + exchange, but dissolved SO4 2− content would remain constant. Consequently, ion exchange and the dissolution of gypsum support the hypothesis of their contribution to the high salinity of groundwater.

Saturation Indices (SI)
Within the studied aquifer, the processes of ion exchange, precipitation, and dissolution are all actively occurring. According to the findings, the majority of the shallow aquifer samples are unsaturated regarding gypsum and anhydrite. It is assumed that unsaturated mineral phases (SI ≤ 0.1) tend to dissolve, whereas oversaturated mineral phases (SI ≥ 0.1) have a tendency to precipitate.
The program used to calculate the saturation indices (SI) was "Diagramme" [48]. The results reveal that SO4 2− concentrations rise when gypsum gradually becomes saturated ( Figure 9). This strong relationship appears to be connected to the geology of the research area, where we can witness dispersed gypsum deposits in the clay formations (see Figures  3 and 10). The presence of organic matter suggests that it may play a crucial function in

Ca 2+ /SO 4 2− Relationship
In line with the previous hypothesis, the reduction in Ca 2+ concentrations in comparison to SO 4 2− in the more sulphate-rich sampled shallow groundwaters, as shown in (Figure 8b), is possibly the result of cation exchange reactions attributed to the presence of a fraction of clay and organic matter within the aquifer matrix as Na + is released, here resulting in calcium concentration as is present here. If we consider the gypsum lenses as an additional potential Ca 2+ source for groundwater mineralization, Ca 2+ might be removed by Na + exchange, but dissolved SO 4 2− content would remain constant. Consequently, ion exchange and the dissolution of gypsum support the hypothesis of their contribution to the high salinity of groundwater.

Saturation Indices (SI)
Within the studied aquifer, the processes of ion exchange, precipitation, and dissolution are all actively occurring. According to the findings, the majority of the shallow aquifer samples are unsaturated regarding gypsum and anhydrite. It is assumed that unsaturated mineral phases (SI ≤ 0.1) tend to dissolve, whereas oversaturated mineral phases (SI ≥ 0.1) have a tendency to precipitate.
The program used to calculate the saturation indices (SI) was "Diagramme" [48]. The results reveal that SO 4 2− concentrations rise when gypsum gradually becomes saturated ( Figure 9). This strong relationship appears to be connected to the geology of the research area, where we can witness dispersed gypsum deposits in the clay formations (see Figures 3 and 10). The presence of organic matter suggests that it may play a crucial function in the ion exchange processes, particularly in deltaic regions [49]. Additionally, neither the groundwater flow direction nor the piezometric map of the unconfined aquifer indicates any probability of seawater intrusion (Figure 4). In some geological formations, the lithostratigraphic investigation reveals the dissolution of existent evaporates [30]. This supposition was primarily tested by geochemical data, which recommended water-rock interactions as a significant source of salinity in shallow groundwaters. This was further supported by the interpretation of the salinity map ( Figure 5) and chemical diagrams that show the interactions between the major elements.
Water 2023, 15, x FOR PEER REVIEW 12 of 23 the ion exchange processes, particularly in deltaic regions [49]. Additionally, neither the groundwater flow direction nor the piezometric map of the unconfined aquifer indicates any probability of seawater intrusion (Figure 4). In some geological formations, the lithostratigraphic investigation reveals the dissolution of existent evaporates [30]. This supposition was primarily tested by geochemical data, which recommended water-rock interactions as a significant source of salinity in shallow groundwaters. This was further supported by the interpretation of the salinity map ( Figure 5) and chemical diagrams that show the interactions between the major elements.  3 and Figure 10). The presence of organic matter suggests that it may play a crucial function in the ion exchange processes, particularly in deltaic regions [49]. Additionally, neither the groundwater flow direction nor the piezometric map of the unconfined aquifer indicates any probability of seawater intrusion (Figure 4). In some geological formations, the lithostratigraphic investigation reveals the dissolution of existent evaporates [30]. This supposition was primarily tested by geochemical data, which recommended water-rock interactions as a significant source of salinity in shallow groundwaters. This was further supported by the interpretation of the salinity map ( Figure 5) and chemical diagrams that show the interactions between the major elements.

Stable Isotopes (δ 2 H, δ 18 O)
In the Sousse region, the Mio-Pliocene aquifer recharge results from precipitation infiltration or can be associated with the infiltration of water used in agriculture activities with an isotopic enrichment due to evaporation (recycled water) or by the called irrigation return flow [50]. Assuming that no isotopic fractionation occurred (e.g., evaporation) the sampled groundwaters should have an isotopic composition similar to the weighted mean isotopic composition of the regional precipitation.
The isotopic composition of groundwater (δ 2 H, δ 18 O) should remain conservative during the groundwater flow [51], and for this reason isotopes ( 18 O and 2 H) are often used as tracers in hydrogeological investigations. In coastal areas dealing with salinization issues [51][52][53], the previous assumption will be extremely important to separate seawater intrusion processes from salt dissolution mechanisms. Several hydrogeological studies performed in coastal areas have demonstrated how isotopes can give unique information regarding groundwater mixing. Several authors have mentioned and used the isotopic composition together with physico-chemical parameters [2,5,26,46,51,54,55].
In our case study, the groundwater isotopic composition was plotted in an orthogonal diagram (Figure 11a), and varied between −5.40 to −3.90‰ in oxygen-18 and from −32.7 to −25.1‰ in deuterium. In Figure 10a, three lines are also plotted, representing the Global Meteoric Water Line (G-MWL) defined by Craig (1961), the Local Meteoric Water Line (Tunis-Carthage data [56]) and the fresh water-seawater mixing line. Analyzing the groundwater sample distribution, two different isotopic evolution trends can be noticed, i.e., a group composed of the most enriched waters, towards the seawater composition, suggesting mixing between fresh and sea water. The second group is composed by the samples plotted between the G-MWL and the regional precipitation (Tunis-Carthage Meteoric Water Line), with a more depleted composition, and with concentrations rather apart from the long term weighted means [56] of Tunis-Carthage precipitation (δ 18 O = −4.38 ± 0.92‰ and δ 2 H = −25.4 ± 6.5‰). An isotopic difference of around 1‰ in oxygen-18 values (mean value) is observed between the Mio-Pliocene waters and modern regional precipitation. The isotopic gap (depletion found in the second group) is not understood and is still in need of an explanation.
Within the environmental behavior of isotopes in precipitation, the isotopic depletion can be associated with different preferential recharge altitudes or with different environmental climatic conditions as in colder periods [57][58][59][60]. Topographic and geographic changes of the studied area do not seem to have a feasible explanation for the observed depletion. Similar isotopic differences have been reported in studies carried out in Western and Central Europe and the United Kingdom, in old groundwater systems recharged during the Last Glacial Maximum (±20,000 years ago). In these areas, a depletion between the groundwater systems and modern precipitation was also noticed [57,58,[61][62][63][64]. Considering these case studies, one possible hypothesis can be formulated to explain the isotopic depletion found in part of the water samples from Sousse Mio-Pliocene aquifer: these waters were infiltrated under different a climatic environment; a colder period, probably during the Last Glacial Maximum (LGM). A difference of 5 • C between LGM and nowadays corroborates the hypothesis that the isotopic depletion found in the groundwater, when compared with modern regional precipitation, was induced by different climatic condition during colder period of recharge [65]. To support the characterization/origin of the salt content in the groundwater, the δ 18 O values were plotted vs. the Cl − content (Figure 11b). A clear trend could not be identified; the samples are dispersed, creating a group. Nevertheless, although dispersed, the hypothesis of mixing with seawater should be considered to explain the isotopic enrichment and salinization. The samples that are low in mineralization and isotopically enriched can be due to evaporation during the sampling or even water evaporation within the borehole.

Chlorine Isotope
δ 37 Cl has been used regularly to understand the sources of the salinity in natural water samples. δ 37 Cl values of seawater and of water affected by dissolution of halite varies between −0.9‰ and +0.5‰ [16]. More recently, a study on δ 37 Cl values from a large set of primary halite samples showed that the δ 37 Cl of halite can be defined more specifically between −0.24‰ and +0.33‰. It also showed that the δ 37 Cl of ocean water has always been close to 0‰ (for the last 2 Ga) [66].

Chlorine Isotope
δ 37 Cl has been used regularly to understand the sources of the salinity in natural water samples. δ 37 Cl values of seawater and of water affected by dissolution of halite varies between −0.9‰ and +0.5‰ [16]. More recently, a study on δ 37 Cl values from a large set of primary halite samples showed that the δ 37 Cl of halite can be defined more specifically between −0.24‰ and +0.33‰. It also showed that the δ 37 Cl of ocean water has always been close to 0‰ (for the last 2 Ga) [66]. δ 37 Cl variations in the aquifer of Sousse ( Figure 12 and Table 1) are large, with δ 37 Cl values ranging from −0.6 to +0.3‰ vs. SMOC. However, most samples show δ 37 Cl values that are within the range of primary halite. δ 37 Cl variations in the study area depend largely on the well location. Only a limited δ 37 Cl range is found in the eastern section of the field area with values between −0.3 and 0‰. More inland δ 37 Cl values have a much larger range. Here, the δ 37 Cl variability is about equal to the full range that has been measured in this study. The data suggests that the δ 37 Cl of groundwater in the vicinity of the coastline is in part influenced by the presence of seawater, while more inland the δ 37 Cl values reflect chemical and physical processes of which ion-exchange and diffusion [67].
Water 2023, 15, x FOR PEER REVIEW 15 of 23 δ 37 Cl variations in the aquifer of Sousse ( Figure 12 and Table 1) are large, with δ 37 Cl values ranging from −0.6 to +0.3‰ vs. SMOC. However, most samples show δ 37 Cl values that are within the range of primary halite. δ 37 Cl variations in the study area depend largely on the well location. Only a limited δ 37 Cl range is found in the eastern section of the field area with values between −0.3 and 0‰. More inland δ 37 Cl values have a much larger range. Here, the δ 37 Cl variability is about equal to the full range that has been measured in this study. The data suggests that the δ 37 Cl of groundwater in the vicinity of the coastline is in part influenced by the presence of seawater, while more inland the δ 37 Cl values reflect chemical and physical processes of which ion-exchange and diffusion [67].

Tritium
Tritium distribution in precipitation is logged from the closest stations of the IAEA/WMO network at the Global Network of Isotopes in Precipitation (GNIP) [56]. For practical purposes, it is possible to make an estimation of the 3 H rain-out at any location by reference to the data (monthly average data for 1968-1997 available from Tunis-Carthage rain station, about 140 km from the study area, and for 1992-2022 from Sfax rain station, located 90 km from the Oued Laya region and at almost the same altitude) [6,56]. Tritium, as a radioactive isotope, is one of the best tools for the estimation of modern groundwater recharge.
In the Sousse region, salinization is likely not from the present-day, since most water samples have 3 H compositions near zero as well as fairly variable chloride contents.
Due to the presence of a clay fraction in the aquifer matrix and the fact that most samples have low tritium contents, we understand that the residence time of the recharge water is rather relevant and favorable to ion exchange process and dissolution of evaporite minerals.
The correlation between 3 H-Cl and 3 H- 18 O demonstrates that recent saltwater intrusion is the cause of the increased mineralization (Figure 13a,b) Consolidating the evaluation of these three variables enabled us to develop the hypothesis that the fundamental process causing the increase in salinization in the groundwater system is the dissolution of evaporite minerals. The stable isotope content in the Tritium distribution in precipitation is logged from the closest stations of the IAEA/WMO network at the Global Network of Isotopes in Precipitation (GNIP) [56]. For practical purposes, it is possible to make an estimation of the 3 H rain-out at any location by reference to the data (monthly average data for 1968-1997 available from Tunis-Carthage rain station, about 140 km from the study area, and for 1992-2022 from Sfax rain station, located 90 km from the Oued Laya region and at almost the same altitude) [6,56]. Tritium, as a radioactive isotope, is one of the best tools for the estimation of modern groundwater recharge.
In the Sousse region, salinization is likely not from the present-day, since most water samples have 3 H compositions near zero as well as fairly variable chloride contents.
Due to the presence of a clay fraction in the aquifer matrix and the fact that most samples have low tritium contents, we understand that the residence time of the recharge water is rather relevant and favorable to ion exchange process and dissolution of evaporite minerals.
The correlation between 3 H-Cl and 3 H- 18 O demonstrates that recent saltwater intrusion is the cause of the increased mineralization (Figure 13a,b) shallow well No. 6 is indicative of recent water mineralization and of an isotopic enrichment higher than that of seawater (Table 1 and Figure 13b).  Table 1 and expressed in percent modern carbon (pmC). 14 C activities in the groundwater of Oued Laya vary between 9 and 56 pmC in the deep aquifer and between 20 and 99 pmC in shallow aquifer.
The δ 13 C values vary between −10.02‰ and −14.6‰ with an average of −11.44‰. The estimation of the "real" groundwater residence times with 14 C consists of calculating the initial activity A0, corrected after the evaluation of the contribution of different carbon sources and isotopic exchange process [68][69][70].
The difficulty with 14 C dating lies in determining the initial activity A0 of the carbon. Several authors have proposed correction models for the initial activity A0, most of which are based on the different mixing processes and on the fact that the Dissolved Inorganic Carbon (DIC) of the waters originates from two carbon sources: CO2 gas and solid carbonates. Groundwater ages have been calculated and compared using eight correction Consolidating the evaluation of these three variables enabled us to develop the hypothesis that the fundamental process causing the increase in salinization in the groundwater system is the dissolution of evaporite minerals. The stable isotope content in the shallow well No. 6 is indicative of recent water mineralization and of an isotopic enrichment higher than that of seawater (Table 1 and Figure 13b).

Carbon 14 and δ 13 C
Carbon-14 ( 14 C) is known as an excellent natural tracer of groundwater flow in aquifers. For this aquifer, 23 dug wells and boreholes were sampled for 14 C and δ 13 C. The results are shown in Table 1 and expressed in percent modern carbon (pmC). 14 C activities in the groundwater of Oued Laya vary between 9 and 56 pmC in the deep aquifer and between 20 and 99 pmC in shallow aquifer.
The δ 13 C values vary between −10.02‰ and −14.6‰ with an average of −11.44‰. The estimation of the "real" groundwater residence times with 14 C consists of calculating the initial activity A 0 , corrected after the evaluation of the contribution of different carbon sources and isotopic exchange process [68][69][70].
The difficulty with 14 C dating lies in determining the initial activity A 0 of the carbon. Several authors have proposed correction models for the initial activity A 0 , most of which are based on the different mixing processes and on the fact that the Dissolved Inorganic Carbon (DIC) of the waters originates from two carbon sources: CO 2 gas and solid carbonates. Groundwater ages have been calculated and compared using eight correction models: Tammers [71], Ingerson and Pearson [72], Mook [73], Fontes and Garnier equilibrium [74], IAEA [75], Evans [76], and Eichinger [77]. The isotopic composition of different carbon used in the calculations are: (1) δ 13 C content of the soil gas: δ 13 C = −21‰ vs. PDB, the average value in δ 13 C obtained allowing measurements carried out on samples of plants in C3 and C4 [78], A 14 C 100 pmC; (2) δ 13 C content in solid carbonates: δ 13 C 0‰ vs. PDB, A 14 C activity of rock carbonates (dead carbon): A 14 C = 0 pmC [79].
The results of the different models have shown that: Tammers, Evans, Pearson, Fontes and Garnier and Eichinger models give negative initial activity. Values lower than the measured activities result in negative ages. However, both Mook and Fontes and Garnier equilibrium models do not differ too much from the Fontes and Garnier model and give logical results. Unfortunately, they present several negative corrected age values. The "IAEA" model presents much more logical results which agree with the other isotope data, and this model was adopted for the correction of the " 14 C ages". The corrected 14 C ages obtained using all models are shown in Table 2. The respective 14 C ages in aquifers could be used for the estimation of groundwater recharge rate (R). Vogel [80] and Cook and Herczeg [81] provided summaries of this approach. Presuming that an aquifer with uniform thickness has a constant replenishment in its unconfined area and has a restriction to not receive any recharge in its confined part of the aquifer, the Vogel model [80] may be used to calculate the recharge rate R from the groundwater ages.
The recharge estimation equation is: where H is the aquifer's overall thickness and ε is its mean porosity. Sample depth Z is the distance from the water table. The age of the groundwater at depth Z is given by t.
The mean aquifer thickness used in the calculation of the recharge rate was 30 m for the Mio-Pliocene aquifer. A mean aquifer porosity of 0.12 was used as a representative value for the porosity. Applying the aforementioned equation to the data revealed the recharge estimated to be between 2 and 12 mm/a.
The youngest ages were observed in the coastal areas, according to the geographical distribution of 14 C corrected ages, confirming the presence of recent recharge of this aquifer ( Figure 13). This map help to clarify conditions of recharge of the aquifer (current location of recharge areas and its evolution vs. time) as well as the general characteristics of groundwater flow (direction and flow velocity) ( Figure 14). The relatively young ages argue for a current significant recharge of rainwater to the water table through the Oued Laya.
The mean aquifer thickness used in the calculation of the recharge rate was 30 m for the Mio-Pliocene aquifer. A mean aquifer porosity of 0.12 was used as a representative value for the porosity. Applying the aforementioned equation to the data revealed the recharge estimated to be between 2 and 12 mm/a.
The youngest ages were observed in the coastal areas, according to the geographical distribution of 14 C corrected ages, confirming the presence of recent recharge of this aquifer ( Figure 13). This map help to clarify conditions of recharge of the aquifer (current location of recharge areas and its evolution vs. time) as well as the general characteristics of groundwater flow (direction and flow velocity) ( Figure 14). The relatively young ages argue for a current significant recharge of rainwater to the water table through the Oued Laya.
In the southeast, the 14 C activities decreased sharply. This is especially true for samples 4, 15, and 8 (Figure 14), where the 14 C activities of the water decreased to very low values presumably reflecting a reserve of old water. These low 14 C activities are not accompanied by changes in the chemical composition of water, confirming the limit of water recharge from precipitation at these sites. The same 14 C results are in good agreement with conclusions reached from stable isotopes and tritium, which strongly suggested that the recharge of Oued Laya is recent and renewed through rainwater. Finally, this map shows that the coastal zone is a new and renewable water area while the Southwest area consists of older waters and their ages exceed 5000 years.

Conclusions
In most aquifer systems located in the coastline of the Mediterranean region, the increase in the groundwater's salinity is usually mainly ascribed to overexploitation. However, this seems not to be the case for the studied groundwaters in the Oued Laya aquifer, based on the results obtained through the multidisciplinary (isotopic and geochemical) approach used in this study. Contrariwise, the dissolution of evaporitic minerals such as halite and gypsum can explain the water salinity of the studied aquifer systems, although In the southeast, the 14 C activities decreased sharply. This is especially true for samples 4, 15, and 8 ( Figure 14), where the 14 C activities of the water decreased to very low values presumably reflecting a reserve of old water. These low 14 C activities are not accompanied by changes in the chemical composition of water, confirming the limit of water recharge from precipitation at these sites. The same 14 C results are in good agreement with conclusions reached from stable isotopes and tritium, which strongly suggested that the recharge of Oued Laya is recent and renewed through rainwater. Finally, this map shows that the coastal zone is a new and renewable water area while the Southwest area consists of older waters and their ages exceed 5000 years.

Conclusions
In most aquifer systems located in the coastline of the Mediterranean region, the increase in the groundwater's salinity is usually mainly ascribed to overexploitation. However, this seems not to be the case for the studied groundwaters in the Oued Laya aquifer, based on the results obtained through the multidisciplinary (isotopic and geochemical) approach used in this study. Contrariwise, the dissolution of evaporitic minerals such as halite and gypsum can explain the water salinity of the studied aquifer systems, although the groundwater mineralization is also strongly controlled by cation exchange mechanisms between clays and groundwaters. The hypothesis of mixing between the shallow groundwaters of the Mio-Pliocene aquifer system with seawater (contributing to the groundwater contamination) seems to be supported by stable isotopic data. However, correlations between Cl − vs. δ 18 O (or vs. δ 2 H) are not considerable, excluding the hypothesis of seawater intrusion as the most important cause of groundwater salinity. The above-mentioned rather low correlations are usually explained by water-rock interaction processes leading to evaporite mineral (mainly halite and gypsum) dissolution, corroborated by higher Na + and Ca 2+ concentrations in groundwaters. In the study region, such processes are related to the presence of saline basins ("sebkhas") and to the type of rocks found below the ground surface, namely interbedded gypsum lenses in the sandstone geological formations. It seems that the recharge water's residence time is rather high, considering (i) the low 3 H concentration in most groundwater samples, (ii) aquifer matrix clay fraction, fostering the dissolution of minerals in evaporite rocks and cation exchange between clays and groundwaters. This research confirms the importance of a multidisciplinary approach to the study of groundwaters in shallow aquifers in semi-arid areas. The decrease in groundwater quality is usually ascribed to seawater intrusion. However, that was not the case in the study area. This indicates that measures to improve the groundwater quality in the case of the Oued Laya aquifers must be different from the measures that are needed in other areas where, for example, seawater intrusion is a major cause of groundwater deterioration. This does not only indicate that studies need to be separately done in each area, but also that measures to improve groundwater quality in a region need to be tailored to the needs of each specific situation.