Hydrogeochemical and isotopic investigation of the multilayer aquifer system: a case study of the Quaternary and the Campanian aquifers of the Mateur plain, Northern Tunisia

In recent decades, the world has been facing major challenges in providing water resources for the growing population (Dai et al. 2004). Water resources have become scarcer due to climate change and the increasing demands for water for soil irrigation, domestic and industrial uses. In addition, groundwater is becoming an increasingly crucial natural resource given the exponential growth of the human population. In Tunisia, nowadays, water availability is around 450 m³/person/year. According to actual previsions, this ratio would be only 315 m³/person/year in 2030 as the number of Tunisian inhabitants is expected to reach 13 million (Ben Ahmed & Ben Rouina 2019).

The Mateur multi-aquifer system, consisting of a Quaternary alluvial aquifer and a Campanian limestone aquifer, is one of the most important aquifers in northern Tunisia as it provides domestic and agricultural water supply for the Mateur region. Nonetheless, so far no precise certainty regarding the interconnection of these aquifers has been definitely established, at least from a geochemical and isotopic point of view. Previous studies have largely focused on groundwater hydrogeochemistry of the alluvial Quaternary aquifer (Tlili-Zrelli et al. 2016), on the spatial-temporal variation of its water quality (Tlili-Zrelli et al. 2018) and on groundwater nitrate contamination (Jendoubi & Bouhlila 2001; Nasri et al. 2014).

Despite the economic importance of groundwater in Mateur plain, little is known about the groundwater geochemistry of the limestone aquifer and inter-aquifer connectivity. Therefore, a holistic research approach is necessary to better ascertain the interactions between both aquifers. To achieve such an endeavor, geochemical and isotopic approaches can be used. Groundwater geochemistry is a valuable tool to determine inter-aquifer mixing and flow regimes, which can provide insights to pursue a flexible and responsive groundwater resources management (Dogramaci & Herczeg 2002; Edmunds et al. 2003; Raiber et al. 2009; Cartwright et al. 2010).

Stable isotopes can be used as a reliable and precious tool to determine groundwater origins and study groundwater evolution couples to hydrochemical, hydrodynamic and statistical methods (Nayak et al. 2017; Cerar et al. 2018). It can offer an evaluation of physical processes affecting water masses, such as evaporation and mixing. In addition, the application of isotope-based methods has become well established for water-resource assessment, development and management in the hydrological sciences, and is an integral part of many water quality and environmental studies (Clark & Fritz 1997; Cook & Herczeg 2000). Furthermore, the deuterium excess calculated as D-excess = δ²H−8δ¹⁸O (Dansgaard 1964) has been widely used in hydrological studies.

Tritium isotope (T) is always introduced into the hydrological cycle and dated with the fallout from atmospheric nuclear weapon tests conducted mainly during the early 1960s. It can be indirectly used to evaluate the rate of groundwater circulation and renewal rate (Clark & Fritz 1997). In deep aquifers, tritium (³H) appears to be a reliable tracer of modern waters. Most of this radioactive H isotope reached the groundwater system during the period of thermonuclear bomb testing (Clark & Fritz 1997).

In general, water with ³H content of <1 TU is regarded as having a pre-1952 age, a date which represents the peak of the artificial release of tritium through nuclear (atomic bomb) tests. Such waters are characterized to have been affected by little or no secondary processes such as evaporation prior to infiltration or isotopic exchange with the aquifer materials (Mazor 1991). However, ³H concentrations above the detection limit (1 TU) indicate recent water infiltration.

The hydrogeochemical and isotopic methods are applied as powerful tools for resolving different groundwater-related problems (Appelo & Postma 2005; Ma et al. 2009), and it has been generally applied to the geochemical investigation in multilayer aquifer for Tunisia (Farid et al. 2015; Mokadem et al. 2016; Hamdi et al. 2018) and in many aquifers around the world (Stimson et al. 2001; Hofmann & Cartwright 2013; Ma et al. 2018; Keesari et al. 2021).

The overall aims of this study are (i) to identify the origin and the chemical behavior of groundwater solutes from the Quaternary alluvial and the Campanian limestone aquifers; (ii) to determine the most relevant processes which control the groundwater geochemistry for each aquifer by using geochemical and isotopic tools and (iii) to appreciate the interconnection between the two aquifer systems which contribute to groundwater management in the Mateur region.

The Mateur plain covers different stratigraphic units ranging from Triassic to Quaternary (Melki et al. 2011). The stratigraphic sequence primarily consists of carbonate formations which characterize the Tertiary deposits of the Ichkeul Mountains. The Santonian outcrops to the southwest of Mateur also form the two limestone hills of Mateur-Ras El Ain (Fournet 2001). It is made up of grey marls, with intercalations of small limestone beds.

Miocene deposits are mainly developed in the East of Mateur plain and the North of Ichkeul lake and are composed of gypsum and marl. The Quaternary and recent alluvium cover the major part of the plain and are composed of silt and detrital formations.

The limestone aquifer is the main productive level which corresponds to Abiod Formation, dated as Campanian–Maastrichtian in age. The hydraulic characteristics of this aquifer are: a transmissivity varying from 33 to 254×10⁻⁴ m/s; a water infiltration speed in the vertical direction less than 10⁻⁷ m/s in the most unfavorable cases and greater than 5×10⁻⁷ m/s in the most favorable cases; a storage coefficient of around 8×10⁻² for the limestone structures of Ras El Ain and 12.5×10⁻² for those of Mateur ville; an average drainance of the order of 6×10⁻⁹.

The limestone aquifer is fed by direct infiltration of rainwater, in areas of fracture-related permeability (limestone hills), or by drainage alluvium, especially during the low water period, with a flow rate of 7 l/s. The drainage is horizontal in limestone outcrops of Mateur and Ras El Ain, and vertical between them (Ennabli 1980).

The upward and diffuse movement of the deep waters and their elimination by the superficial capillary fringe constitute the main factor of discharge of the deep limestone aquifer. This aquifer is exploited by boreholes of the National Water Supply and Distribution Company (SONEDE) to meet the region's drinking water needs (DGRE 2010). Groundwater from the alluvial aquifer is used for domestic and agricultural activities with 49 and 51%, respectively, whereas the limestone aquifer is mostly used for local domestic water supply (97.3%) owing to its good groundwater quality (DGRE 2005). Most soils in the Mateur plain are fluvisoils or hydromorphic and halomorphic soils with high salinity (Supplementary file, Figure 2) (CRDA 2007).

A geochemical and isotopic sampling campaign was conducted during September 2018. A total of 27 samples were collected from the Mateur aquifer and are distributed as follows (Figure 1):

• - Twenty-two samples from pumping wells and boreholes exploiting the Quaternary aquifer.

• - Five samples from boreholes exploiting the Campanian aquifer.

pH, electrical conductivity (EC), dissolved oxygen and temperature were measured in the field. The samples were kept at 4 °C for their subsequent chemical analyses. Afterwards, the samples were filtered through a 0.45 μm Millipore filter. Chloride was determined by standard AgNO₃ (Rodier 1984). Sulfate (SO₄²⁻) content was measured by the gravimetric method using BaCl₂. Sodium (Na⁺) and potassium (K⁺) were measured by flame photometry and calcium (Ca²⁺) and magnesium (Mg²⁺) with atomic absorption. Alkalinity was determined by titration with HCl. The quality of chemical analysis was checked by making an ionic mass balance, accepting an error rate lower than 5%. The obtained results were applied to a statistical study using Andad software (CVRM 2000). Surfer 8 (Surfer^® (Golden Software, LLC)) was used to generate TDS spatial distribution map in the study area. The groundwater samples were classified into water types in the Piper diagram based on the ionic composition of different water samples, where ion concentrations are expressed in milliequivalent per liter (Piper 1944).

In order to investigate thermodynamic controls on mineral–water interactions, the geochemical computer code PHREEQC (Parkhurst & Appelo 2013) was used to calculate mineral saturation indices which were used to find out which minerals could possibly dissolve or precipitate in the aquifer system. Positive values indicate precipitation or a stable condition for minerals; zero value indicates that the solution is in equilibrium with a mineral; and negative values indicate dissolution of a mineral (Lee & Gilkes 2005).

Hydrogen and oxygen isotope analysis were performed at the Laboratory of Radio-Analysis and Environment (LRAE) of the National School of Engineers of Sfax (ENIS) using the Laser Absorption Spectrometer LGR DLT 100 (Penna et al. 2010). The results were expressed as δ (‰) relative to the international standard V-SMOW (Vienna Standard Mean Ocean Water). Twelve samples were selected for tritium content analyses using electrolytic enrichment and liquid scintillation technique (Taylor 1976). ³H concentration is expressed in Tritium Units (TU). One TU is defined as the isotope ratio ³H/¹H = 10⁻¹⁸. Typical precisions are ±0.2, ±1.0‰ and ±0.3 TU for oxygen-18, deuterium and tritium, respectively.

Hydrochemical characteristics of groundwater in the study area are summarized in Table 1.

The temperature values of groundwater varied from 17.9 to 24.5 °C and from 18 to 20.0 °C in the alluvial and limestone aquifer, respectively.

The groundwater samples in the limestone aquifer are characterized by a moderate mineralization. TDS values of all groundwater samples vary between 1,518 and 1,932 mg/l and are lower than those recorded in the Quaternary alluvial aquifer. The highest TDS value is recorded in groundwater sample C1 located north of the hill of Ras El Ain, characterized by a thickening of alluvium (about 50 m thick), surmounting the limestone reservoir, which has an effect on the replenishing rate of water and consequently on the salt load. Low TDS value is recorded in the groundwater sample C5 situated near the Campanian hill of Ras El Ain, which constitutes a safe water supply area.

The sulfate contents varied between 146.9 and 5,025 mg/l and between 293.8 and 773.8 mg/l in the alluvial and in the limestone aquifers, respectively.

In groundwater samples from the Campanian limestone aquifer, we note the dominance of chlorides and sulfates among the anions with 60 and 28% of all anions, respectively. The sodium and calcium are dominant among the cations with 45 and 37% of all cations, respectively.

Regarding the limestone aquifer, groundwater samples have a mixed facies Ca–Na–SO₄–Cl to sodium chloride facies (NaCl).

The Quaternary aquifer samples’ isotopic compositions range from –5.2 to –1.0‰ and from –31.5 to –13.4‰ for δ¹⁸O and δ²H, respectively (Table 2). Those of the Campanian aquifer vary between –5.8 and –4.9‰ for oxygen-18 and from – 36.1 to –29.7‰ for deuterium.

In Tunisia, tritium content in rainfall is monitored by the International Atomic Energy Agency (IAEA 2015) at two stations (Tunis and Sfax). The annual natural ³H content for the Bizerte region was defined for the first time by Ben Ammar et al. (2020), from the Sfax chronicle, which showed a stabilization of the annual ³H signal over the period from 1992 to 2013 with an average of 3.87±0.65 TU (Ben Ammar et al. 2020).

In the study area, tritium activities range from 1.95 to 5.05 TU in the Quaternary aquifer, except for samples where tritium contents are below 1 (0.79 and 0.93, respectively), indicating a recharge prior to the era of thermonuclear bomb testing (Bajjali et al. 1997). Concerning the Campanian aquifer, tritium contents are 1.64 in groundwater sample C2 and 2.67 in groundwater sample C3.

The spatial variation of temperature is linked to the water recharge and to the infiltration transfer time, which in turn both depend on porosity, lithology and thickness of the unsaturated zone (Tlili-Zrelli et al. 2018). The spatial variation is marked by a decrease as the depth increases.

In fact, the water mineralization of the Campanian aquifer is probably not solely linked to the dissolution of calcite. This implies that another potential source of excess Ca and SO₄ may be linked to the dissolution of gypsum from the alluvial aquifer. Consequently, the dominance of Ca⁺⁺ and SO₄⁻ in groundwater samples from the limestone aquifer may be ascribed to vertical mixing with the shallow quaternary groundwater which led to the change of hydrochemical facies of the Campanian aquifer from Ca–HCO₃ to Na–Ca–SO₄–Cl.

Another hypothesis that can be formulated is that the high sulfate in deep groundwater may be linked to the dissolution of marly intercalations in Campanian limestones.

The detailed analysis of stable isotope data (Figure 5) allows the distinction of two groups of waters:

• - The first group includes water samples representing limestone aquifer and the upstream part of the Quaternary aquifer which are lying widely on the GMWL. These samples are characterized by relatively depleted oxygen-18 and deuterium contents, indicating that they are not significantly affected by evaporation, which imply rapid infiltration of rainfall waters and a low amount of evaporation.

• - The second group comprises groundwater samples below the GMWL more enriched in δ¹⁸O and δ²H which are located for the most part in the downstream part of the Quaternary aquifer.

The enrichment in heavy isotopes is interpreted as an evaporative influence which is controlled by (i) thickness of the unsaturated zone which decreases globally from upstream to downstream part of the aquifer (Supplementary file, Figure 3), indeed, the rate of evaporation is the highest as the thickness is the lowest; (ii) the nature of the sediments of the unsaturated zone in which clay proportion also increases from upstream to downstream (Supplementary file, Figure 4), and therefore influencing the rate of infiltration and (iii) location of irrigated perimeters allowing irrigation return flow.

Furthermore, the insignificant difference between δ¹⁸O and δ²H values across the two aquifers confirms that the Mateur aquifer system is a hydrogeological entity formed by the superficial and deep levels communicating by drainage.

An attempt was made in this study to establish an improved understanding of the factors and mechanisms controlling the groundwater chemistry in Quaternary alluvial and Campanian limestone aquifers, and their mixing and evolution in the Mateur region. Hydrogeological, hydrochemical and isotope tools were used herein to identify such mechanisms.

In the Quaternary alluvial aquifer, the dominant geochemical processes influencing the groundwater mineralization are: (i) water origin (leaching of low mineralized water in the recharge zone and salt water into the discharge zone); (ii) dilution in the wet season and evaporation in the dry season; (iii) thickness, grain size and lithology of the unsaturated zone. It is noteworthy that the water–rock interaction does not contribute significantly to the solutes acquisition since minerals within the aquifer are sparingly soluble.

Regarding the Campanian limestone aquifer, groundwaters have Na–Cl or Na–Ca–Cl–SO₄ facies. Such facies types indicate that groundwater mineralization from the Campanian aquifer is not under the dominant influence of the water–rock interaction; it may be linked to vertical mixing with the shallow quaternary groundwater, more enriched in sulfate. This led to a change in the hydrochemical facies of the Campanian aquifer from Ca–HCO₃ to Na–Ca–SO₄–Cl.

Furthermore, the δ² H and δ¹⁸O values of groundwater from the two aquifers do not have a greater difference which confirms that the Mateur aquifer system is a hydrogeological entity formed by superficial and deep levels communicating through leakage.

The enrichment in heavy isotopes is interpreted as an evaporative influence which is controlled by (i) thickness of the unsaturated zone which decreases globally from upstream to downstream and where the evaporation rate is higher as the thickness is low; (ii) the nature of the sediments of the unsaturated zone in which clay proportion is increasing from upstream to downstream, and consequently influencing the rate of infiltration and (iii) location of irrigated perimeters allowing irrigation return flow.

The tritium contents for the Campanian aquifer are lowest than those of the Quaternary aquifer, indicating a relatively old age or mixture with recent waters. This confirms the hydraulic continuity between the Quaternary aquifer and the Campanian aquifer inducing the mixing of waters of the two aquifers. This is also consistent with the geochemical analytical results. Lastly, we believe that the present study will have significant implications for a better management of water resources in one of the most important aquifers of northern Tunisia.

The authors gratefully acknowledge the contributions of the members of the Bizerte Water Resources Division/Agriculture Ministry, for their help in the field campaign. Special thanks are due to Professor Kamel Zouari and the technical staff at the Laboratory of Radio-Analyses and Environment of the National Engineering School of Sfax (ENIS) for isotopic analyses. Finally, we wish to thank the editor and the anonymous reviewers for their insightful comments, which greatly helped to improve the paper.

All relevant data are included in the paper or its Supplementary Information.

The authors declare there is no conflict.