Shallow groundwater salinization risks under climate change: FEFLOW model application in Cebala Borj-Touil irrigation zone, northern Tunisia

Effects of climate change (CC) on the quality and quantity of water resources are becoming more and more visible (Barica 1972; Leavesley 1994; Colombani et al. 2016; Mastrocicco & Colombani 2021). Arid and semi-arid regions are very vulnerable to these changes due to the high evaporation, which leads to a progressive depletion of aquifers, particularly shallow groundwater (Barica 1972). Coastal aquifers are highly vulnerable because of changes in groundwater recharge and seawater intrusion (Wang et al. 2023; Ansari et al. 2024). In the case of irrigated perimeters, the anthropogenic processes contributing to the groundwater salinization are mainly summarized in the rise of the groundwater table due to inappropriate irrigation, the failure of the drainage system, and the disturbance of the natural cycle of the soil leaching by the creation of wells (Foster & Chilton 2003). Several studies analyze the relationship between irrigation water quality and groundwater salinization in similar cases (El oumlouki et al. 2018; Sefiani et al. 2019). In addition, the choice of climate models is a major challenge. Information about CC at the regional scale is essential to obtaining accurate results about the environmental impacts of CC. Recently, scaling-down methods have enabled the extraction of results from global climate models at finer scales (Giorgi 2008). The simulations of the different climatic scenarios confirm the high coastal groundwater vulnerability to salinization (Barica 1972; Leavesley 1994; Colombani et al. 2016; Mastrocicco & Colombani 2021). The impact of the different combined processes on water salinization can be analyzed with numerical modeling (Li et al. 2019), even for complex cases or when observed data are scarce and poor (Baccouche et al. 2024). To analyze the salinization of the shallow groundwater, the perimeter of Cebala Borj-Touil was selected, which is irrigated by treated wastewater. In Cebala Borj-Touil, it is not confirmed that irrigation is the origin of the groundwater salinization. Only a few studies have been conducted in this region due to a lack of data on the salinity in the aquifer. However, all studies confirm the presence of extremely salty shallow groundwater (Hachicha et al. 1990; Hachicha & Job 1994; Hassine 2005; Dahmouni et al. 2017; Rzigui et al. 2022a, 2022b). To understand these developments, models are used to address data quality issues and create a reliable model of salinity that can be used for predictive simulations. The FEFLOW model was selected for a 2D simulation of flow and salt transport under complex initial and boundary conditions. The modeling approach is divided into two parts: the flow model and the salt transport model. The transport model depends on the calibrated and validated flow model and has several uncertainties. In the study of Bahri et al. (2015) and Rzigui et al. (2022a, 2022b), the numerical models MODFLOW and FEFLOW were applied for the simulation and optimization of the Cebala Borj-Touil perimeter drainage system. The results confirmed that the drainage system is not working anymore. Water quality and the impact of CC on groundwater salinization remain unknown. The objectives of this paper are the modeling of the variation of the salinization of the shallow groundwater and the estimation of its evolution under CC conditions. This research presents an innovative methodological framework that addresses issues related to data scarcity by utilizing only the available data and numerical models to investigate the impact of CC on a complex deltaic coastal shallow groundwater scenario.

Since 1990, irrigation and drainage systems have been installed. The efficiency of the drainage network has progressively deteriorated over the last three decades. The natural drainage system is composed of two rivers, the Khelidj and the El Maleh wadis. The climate of the study area is semi-arid, with an average annual rainfall of 450 mm and a potential evapotranspiration of 1,300 mm·year⁻¹. The Cebala soils are composed of loamy clay on Quaternary alluvial deposits. The soil salinity varies between 2.0 and 3.5 dS·m⁻¹ at the soil surface layer.

The shallow groundwater table ranges from 0.2 to 1.8 m above sea level (ASL), with a high electrical conductivity varying from 10 to 20 dS·m⁻¹.

The first step in setting up the groundwater model involves collecting geological and hydrological data to create a representative conceptual model. The database included climate information along with data on aquifer characteristics and boundaries, well extractions, irrigation doses, the depth of the water table, and groundwater quality. The data collection process involved sources such as research and regional institutes from the Tunisian Agricultural Ministry, the National Meteorology Institute (NMI), and the National Office of Mines, in addition to relevant research related to the study area. The NMI data include parameters like wind speed, air temperature, humidity, rainfall, and evaporation. These measurements refer to the climatic stations at Cherfech and Carthage and have formed a data series since 1970. Rainfall and evaporation measurements were used in the calculation of the groundwater net recharge. The synthesis of geological data was based on previous studies and focused on the characterization of the stratigraphic units of the site to define the geometry of the aquifer. Based on previous work and monitoring of the groundwater and soil parameters at the laboratories, the characteristics of the aquifer were defined: permeability, total porosity, root depth, and piezometric variations. Hernot's (1949) study was the basis for developing the reference state of the water table. In terms of quality, the investigation focused on the historical trends of groundwater and irrigation water salinity. Climate data are available on a daily basis, while piezometric and salinity measurements are seasonal and less consistent throughout the observation period from 1990 to 2007. Time series and spatial data analysis were performed using R 3.6.2 (R Core Team) and ArcGIS (Esri) (10.2.2).

The model boundaries have been modified and extended to reliable boundary conditions.

In this research, a high-resolution mesh was generated for the drainage network. The small size of the elements is necessary to achieve a more accurate solution at the nodes defined by drainage conditions (Rao 2011). Also, FEFLOW has the best performance to interpolate maps within the finite elements using three interpolation methods with their respective property settings: inverse distance, kriging, and Akima.

The results of the flow model simulation were presented in a previous paper that aimed to evaluate the drainage system performance and the scenarios for the rehabilitation of the perimeter (Rzigui et al. 2022a, 2022b). The flow model simulation shows that a water logging problem is present in the major parts of the perimeter due to the drainage system deficiency and the low slope. Several management scheme scenarios were simulated. The three principal scenarios include: the drainage system efficiency improvement by increasing the drain conductance from 10⁻⁵/s to 10⁻⁴/s (scenario 1), the reduction of the water irrigation amount by applying a revised net recharge for the simulations where the new recharge rates are determined based on a 70% decrease in irrigation (scenario 2), and the third scenario illustrated the combination of the two prior scenarios. The final scheme leads to resolving the problem, even in the lower parts.

The salt transport model is based on the flow model. Under transient conditions, the transport parameters are porosity, molecular diffusion, and longitudinal and transverse dispersivity. Salt transport in the groundwater is considered non-reactive. Boundary conditions are defined by a Dirichlet condition at the North boundaries (Medjerda wadi) and the east boundary (sea) with a salinity of 1.5 and 35 g/l, respectively. A homogenous initial concentration of 3.5 g/l has been assigned to the entire domain. This salinity corresponds to the salinity of the irrigation water. The aquifer porosity varies from 0.3 to 0.35. The α dispersion of the medium is defined by the ratio of dispersion to pore velocity. Two components of dispersivity are distinguished: longitudinal dispersivity (α_L: in the direction of flow) and transverse dispersivity (α_T: perpendicular to the flow direction). The ratio of transverse to longitudinal dispersivity is estimated to be 0.1. A longitudinal and transverse dispersivities of 100 and 10 m, respectively, have been used. Molecular diffusion coefficient values range from 10^–10 to 10⁻¹¹ m²/s for a non-reactive chemical species in clay placers.

Statistical parameters served as error criteria in the calibration and validation of the flow model, allowing for an evaluation of the fit between observed and simulated groundwater data. The calibration process of the flow model was achieved with a correlation coefficient of 0.98 and an root mean square (RMS) error of 0.51 m, indicating an appropriate match between observed and calculated heads in 1949. A transient simulation was first run in 1990. The obtained correlation coefficient of 0.85 shows a significant relationship between simulated and observed hydraulic levels all over 15 piezometers. After calibration, the flow model was verified using a transient simulation from 1990 to 2017 showing a good agreement between the computed and the observed heads of the transient models for the observation period 1990–2007 with a calculated mean RMS error of 0.87 m. Due to the convective salt transport in the Cebala Borj-Touil aquifer, a steady-state solute transport model calibration is not possible. The groundwater flux defined by the flow model is the main factor in the salinity distribution. The calibration of the solute model was conducted using a transient calibrated flow model based on the reference salinity map (Section 2.2).

The water table salinity within the perimeter varies between 6.7 and 15 g/l. The salinity tends to increase from west to east (boundary with the sea). The minimum salinity concentrations of 6.7 g/l are observed locally at the entrance of the Khelidj wadi.

CC has an impact on rainfall distribution and sea level rise. Both factors strongly affect the coastal irrigation perimeter of Cebala Borj-Touil. The impact of different CC scenarios on the evolution of groundwater salinization was simulated based on the calculated salinity map for 2018. The scenarios are based on the Intergovernmental Panel on Climate Change (IPCC)'s (2014) reports on CC in the Middle East-North Africa region. Three main representative pathways for the evolution of greenhouse gas concentration in the 21st century (RCP2.6, RCP4.5, and RCP8.5) are considered in the 5th IPCC assessment report. The most pessimistic scenario RCP8.5 was adopted. Climate scenarios were simulated in the short term, corresponding to 2050, due to the need for a fast forecast of the salinity evolution of the groundwater. The three main scenarios are: (1) a 15% decrease in precipitation (scenario 1), (2) a sea level rise of 0.2 m (scenario 2), and (3) a degradation of irrigation water quality by increasing its salinity to 5 g/l (scenario 3). Simulation scenarios resulting from the combination of scenarios 1 and 3, as well as scenarios 2 and 3 aim to evaluate the most severe influence of CC on groundwater salinization (scenario 4).

The precipitation decrease increases salinity over almost the entire water table, except in the northwest boundary of the perimeter where the average rise is about 1 g/l. However, the increase is greater in the remaining areas of the perimeter, with a maximum increase in the southeast part where the calculated salinity exceeds 15 g/l. This can be explained by the lack of a drainage network and the accumulation of water in those areas. This scenario shows a considerable impact of CC on the general salinization of the aquifer.

The absolute value of the increase depends on the distance of the observation point from the sea. Figure 7 shows the value of two selected piezometers (piez) 2 and 17, located, respectively, 11 and 6.5 km from the sea.

The closer the piezometer is to the sea, the greater the impact of rising sea levels with an increase in water salinity of 66 and 74% for the observation points piez 2 and piez 17, respectively.

Irrigation with salty water will increase the groundwater salinity of piez 2, piez 28, and piez 41 to 148, 128, and 53%, respectively (Figure 8). These results are alarming but remain not realistic because irrigation will stop as soon as the soil salinity exceeds the salt tolerance of the plants.

Simulations were performed on additive scenarios that merged the degradation of irrigation water (scenario 3) with the two previously analyzed scenarios 1 and 2.

In this paper, the salinization of the shallow groundwater in the irrigation perimeter of Cebala Borj-Touil was analyzed. In coastal areas, the salinization of aquifers is also linked to natural factors such as the impacts of CC (Greene et al. 2016). The dynamics of the groundwater salinity in Cebala Borj-Touil are complex. In all scenarios, an increase in the salinity of groundwater from 17 to 150% is observed. Simulation results show good agreement with the literature. In coastal areas, groundwater salinity would increase by 14% in 2050, due to CC impacts (Islam et al. 2020). This alarming development is caused by the coastal location, the semi-arid climate, and the increased anthropogenic factors. Several studies (Akbari et al. 2020; Castaño-Sánchez et al. 2020; Soundala & Saraphirom 2022) highlight that the repercussions of sea level rise are considerably more pronounced. The findings correspond with the outcomes observed in the Cebala case study. The impact of reduced rainfall will be evident. The salinity of Cebala groundwater exceeds 10 g/l. An increase in the frequency and intensity of rainfall enhances the downward movement of chemicals from the surface and vadose zone, leading to a greater influx of suspended and dissolved solids into aquifers (Dao et al. 2024). CC will lead to an increase in soil and water salinity, which will affect irrigation water quality. The decline in water irrigation quality highlights a crucial potential impact on the salinization of the groundwater in the Cebala soil. Scenario 3 suggests that salinity could reach a maximum of 20 g/l by 2050.

The combined scenarios (scenario 4) indicate that the maximum soil salinity reaches 25 g/l. Such a condition may result in unrealistic findings, as the negative impact on production would necessitate the cessation of irrigation. Such a situation is present in various deltaic areas across numerous semi-arid and arid areas. The Cebala case study exemplifies the complexities involved. Therefore, these results should be generalized with caution. The distinct nature of the region and the selected model may result in particular outcomes. Despite its limitations in illustrating salt chemistry in soils, the model is still a practical tool to deal with data scarcity. The FEFLOW model demonstrates a strong capability to simulate groundwater salinization in various climate scenarios (Lam et al. 2021; Yang et al. 2023). The validation of the model's results remains challenging because of the issues of data scarcity. Also, in all studies concerning the impacts of CC, there is always an inherent uncertainty in climate scenarios due to an insufficient understanding of greenhouse gas emissions (Xepapadeas 2024). According to Nouaceur & Murărescu (2016), current climate data from Morocco, Algeria, and Tunisia do not match the forecasts related to CC. Nevertheless, the findings from this research on the Cebala Borj-Touil case remain pertinent and advantageous for the investigation of irrigation perimeter rehabilitation.

The main objective of this work is the development of a transport model for the salinization of the Cebala Borj-Touil perimeter and its future evolution under the impact of CC. The validated flow model made it possible to map the spatial distribution of salinity in 2018, which confirms the persistence of widespread salinization throughout the perimeter. The results of the climate scenarios simulations show that the salinity of the groundwater tends to increase particularly under the effect of sea level rise. A decrease of 15% in precipitation leads to an increase in salinity across nearly the entire aquifer, averaging 10–15 g/l. The rise in sea level has a pronounced effect on the salinization process. Seawater intrusion leads to a water concentration between 16 and 20 g/l. The impact of irrigation is very pronounced. High salt concentrations in irrigation accelerate the salinity levels in groundwater, with increases exceeding 100%. All combined scenarios involving the use of irrigation water at a concentration of 5 g/l show a groundwater salinity level greater than 20 g/l by 2050 in the entire area. Although these results are alarming, they remain largely hypothetical, as irrigation is usually stopped when soil salinity begins to affect plant growth negatively. These results can serve as a basis for decision-making about the rehabilitation of the perimeter. However, it is recommended to conduct more sophisticated studies aimed at reducing uncertainties by examining additional factors and updating the groundwater salinity measurements.

We would like to acknowledge DHI-WASY for providing a FEFLOW license for the current study.

No funding was received for conducting this study.

R.H., M.H., and W.B. contributed to the study's conception and design. Material preparation, data collection, and analysis were performed by R.H. and H.G. The first draft of the manuscript was written by R.H. and all authors commented on previous versions of the manuscript. All authors read and approved the final manuscript.

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

The authors declare there is no conflict.