The study area encompasses the Merguellil watershed, a significant watershed in semi-arid Tunisia, located west of the Kairouan governorate in central Tunisia (Figure 1). The Merguellil Wadi emanates from the elevated regions of Kesra and Makthar, traversing in a northwestern to southeastern direction towards its natural terminus, the El Kalbia Sebkha, adjacent to the northwestern periphery of Kairouan city (Kingumbi, 2006). The studied watershed has two distinct regions: a mountainous upstream area and a flat downstream plain, connected by the El Houareb dam, built in 1989 (Lazard, 2013). The Merguellil watershed features a mix of hilly terrain, including Jebel Kesra, Barbrou, and Trozza, along with expansive plains like El Alaa in its upper reaches. In contrast, the lower section is a broad alluvial plain primarily used for agriculture (Saadi, 2018; Jerbi, 2018).

The Merguellil watershed experiences a semi-arid climate (Henia, 2003) characterized by substantial year-to-year variations in rainfall, notable temperature disparities, and pronounced summer dry spells (Lacombe et al., 2008; Ben Ammar et al., 2006). Between 1980 and 2020, the region’s average annual rainfall was approximately 300 mm, varying from 265 mm in the plains to 515 mm in the highlands. The average annual temperature was 19.5°C, with daily temperatures typically falling below 12°C in January and exceeding 28°C during July and August. The local weather patterns are predominantly influenced by the convergence of disparate air masses emanating from the temperate Northern Tellian mountains and the arid southern regions (Saadi, 2018). During the winter months, the prevailing winds emanate from the north and northwest (Bouzaiane and Lafforgue, 1986), whereas during the summer months, they predominantly originate from the south and southwest. In the study area, in the year of 2020, the annual soil erosion rate was 16 t/ha/year and the sediment yield was 8.95 t/ha/year (Hermassi et al., 2023).

Using SWAT model, the first step of simulation process was the division of the study area into sub-basins (Neitsch et al., 2009). Next step was the elaboration of HRUs, which were established using combinations of slope, LULC and soil type datasets, with a threshold value of 10% applied to each dataset for consistency (Arnold et al., 2012). Moreover, all the necessary climatic variables were precisely prepared in a appropriate format and supplied as inputs to the model (Winchell et al., 2007). For this case, a monthly simulation has taken place, from 2000 to 2020, with an initial one-year warm-up period incorporated into the analysis.

The Sequential Uncertainty Fitting-2 (SUFI-2) algorithm was employed for model calibration and validation (Abbaspour et al., 2007). Sensitivity analysis and calibration, which are closely related, are crucial steps in effectively configuring the SWAT model (Abbaspour, 2015; Kouchi et al., 2017). SWAT-CUP’s primary objective is to identify the most suitable combination of pertinent parameter values following the initial sensitivity analysis. This selection aims to accomplish a strong alignment between the observed data and the simulated values throughout both processes. In this study, model calibration was conducted for the period from 2000 to 2012, while validation covered 2013 to 2020. Sensitivity analysis, a crucial phase in identifying the SWAT input parameters that most affect model output variability (Baker and Miller, 2013), was performed using the SUFI-2 algorithm. This analysis evaluated each parameter using two key indicators: the t-stat and the p-value (Abbaspour, 2015). The parameter with the highest absolute t-stat and the smallest p-value was considered the most sensitive (Abbaspour et al., 2007). In this specific context, a total of 7 parameters were selected for sensitivity analysis. These parameters were explored by taking into account the physical properties defined by the guidelines provided by Arnold et al. (2012). After 1,000 iterations, parameter sensitivity analysis was completed. In fact, 1,000 iterations in SWAT model calibration provide a balance between thorough parameter space exploration, reliable model performance, and reasonable computational effort, making it a widely accepted default in hydrological modeling. The sensitivity analysis was carried out on a monthly basis, spanning the period from 2000 to 2020. The initial year was a warm-up period to initialize unknown variables.

The effectiveness of a hydrologic model is gauged by how well it aligns its simulated outputs with observed data. In fact, according to Abbaspour (2015), the assessment of hydrologic models is significantly contingent upon two crucial steps: sensitivity analysis and parameter calibration.

In this study, the SWAT Calibration and Uncertainty Procedures (SWAT-CUP) software was employed to facilitate automatic calibration process (Abbaspour et al., 2007). This process was supported by the SUFI-2 algorithm (Abbaspour, 2012), a widely recognized tool for calibrating watershed hydrology models (Yang et al., 2008). SUFI-2 is adept at accounting for numerous sources of uncertainty, encompassing parameter uncertainty, conceptual model uncertainty, and input data uncertainty (Gupta et al., 2009). Two key performance indicators, recommended by SUFI-2, were used to evaluate the model’s performance during calibration and validation. These indicators include Nash–Sutcliffe Efficiency (NSE), Coefficient of Determination (R²), the ratio of the root mean squared error to the standard deviation of measured data (RSR) and percent bias (PBIAS).

These metrics collectively assess the model’s ability to replicate observed hydrologic behavior, indicating its accuracy and reliability.

To develop an accurate watershed model, initiating a calibration process is essential. Given that model parameters represent various natural processes and that many parameters collectively affect watershed outcomes, conducting a sensitivity analysis before calibration is crucial. The sensitivity analysis determines which parameters most affect key variables like runoff and refines the calibration process by concentrating on the most influential parameters while excluding less sensitive ones. A global sensitivity analysis was conducted utilizing the SWAT-CUP model. According to the statistical ranking of the model outputs, six parameters were identified as sensitive. The assessment of parameter sensitivity primarily relied on two statistical criteria: the p-value and the t-stat. Among these seven parameters, six emerged as the most sensitive, meeting the criteria of possessing absolute t-stat values equal to or surpassing 2, coupled with p-values less than 0.05. The key parameters identified included the curve number (CN2), SOL_AWC and the ESCO factor. Parameters delineated in Table 1 play a key role in influencing the attitude of runoff within the watershed.

Calibration and validation processes were employed to assess the impact of changes in the SWC area on runoff, soil erosion and groundwater recharge. In this study, a concept of multi-objective calibration in SWAT was applied and revolved around the simultaneous calibration of SWAT model across the watershed and its sub-watersheds using multiple objectives criteria. To assure a multi-objective calibration study, we use multiple criteria such as streamflow, sediment yield and groundwater recharge. Moreover, an objective function which are the performance indicators (NSE and R²) are often used. On the other hand, a multi-watershed calibration has taken place where the SWAT model is calibrated for the entire watershed and some sub-watersheds (SW) simultaneously. In fact, this approach search to ensure that the model performs well across various metrics, not just one.

The calibration was carried out for the period spanning from 2000 to 2012, while the validation phase encompassed data from 2013 to 2020. The study assessed the reliability of parameters that had been calibrated and validated using two different approaches. To evaluate streamflow consistency, observed and simulated monthly flows were compared graphically (Moriasi et al., 2015) during calibration and validation (Figure 2). The findings demonstrated a strong agreement between the simulated and observed values during calibration and validation phases in the outlet of the studied watershed. In terms of evaluating the model’s performance, various statistical metrics (Moriasi et al., 2015) were employed for both calibration and validation stages, using flow data. These indicators included NSE and R², RSR and PBIAS (Table 2). During the calibration phase, the performance indicators produced values of 0.82 for NSE 0.9 for R², 0.19 for RSR and 11.62% for PBIAS. As for the validation process, there was a strong correlation, with NSE reaching 0.81 and R2 at 0.9. Moreover, for the validation process RSR was equal to 0.22 and PBIAS equal to 10.96%. Consequently, the statistical indices showed a very good agreement between the observed and simulated streamflow for the calibration and validation processes.

In the next step, results of multi-watershed calibration were analyzed in two representative sub-watersheds (SW 8 and SW 16). Graphical comparison presents good agreement between simulated and observed values in SW 8 with a NSE equal to 0.75 and R² equal to 0.85, and in SW 16 with a NSE of 0.53 and R² of 0.57, for the period 2002–2008 (Figure 3). The strong R² value in SW 8 (0.85) confirms a high correlation between observed and simulated streamflow, indicating that the model effectively captures the variability and trends in this sub-watershed. The combination of a high NSE (0.75) and R² shows that both the magnitude and timing of peak flows, as well as baseflow conditions, are well simulated. This can likely be attributed to homogeneous watershed characteristics and fewer anthropogenic impacts in SW 8. In contrast, SW 16 exhibits moderate model performance, with a lower NSE (0.53) and R² (0.57) (Figure 3). These values suggest that while the model captures some aspects of streamflow variability, it struggles to accurately represent the dynamics of peak flows and low-flow periods. The discrepancy may stem from complex watershed characteristics in SW 16, such as mixed land use, varying soil properties and concentration of SWC techniques. Despite the weaker performance in SW 16, the calibration across both sub-watersheds demonstrates the model’s ability to handle diverse conditions within the multi-watershed framework.

For the studied area, the outlet was the Houareb reservoir, where only one bathymetric survey was conducted in 2011, providing a sediment yield value of 9.39 t/ha/year. Subsequently, a multi-watershed analysis was performed for sediment yield across several representative sub-watersheds (SW) in the region. Due to limited data availability, this investigation was restricted to four small dams located in SW 3, SW 7, SW 8, and SW 24, for the year 2011. The simulated and observed sediment yield values for these sub-watersheds are as follows (Table 3).

A comparison of simulated and observed sediment yield values reveals acceptable model performance. In fact, the model accurately simulates sediment yield in SW 7 and SW 8 however, it shows slight overestimations in SW 3 and SW 24. The observed value of 9.39 t/ha/year from the Houareb reservoir bathymetry falls within the range of simulated values, confirming that the model captures the overall trend of sediment yield, although local adjustments may be needed for more accurate predictions.

The modeling of groundwater recharge is of paramount importance in understanding and managing water resources effectively. Indeed, SWAT model is a powerful tool that provides quantitative insights into the dynamics of water movement within the watershed.

This study reveals a significant correlation between the simulated values generated by the SWAT model and the observed groundwater recharge values. The strength of this correlation is exemplified by a coefficient of determination reaching (0.8) (Figure 4).

The surface runoff distribution in the Merguellil watershed exhibits a clear lack of uniformity, indicating a non-uniform distribution across its subbasins. To assess the effect of these techniques, the SWAT model was run with and without SWC managements. Figure 5 displays the percentage of difference in runoff values between the treated watershed and the untreated condition.

Figure 3 provides a comprehensive overview of the runoff percentage differences observed in a watershed that has undergone SWC techniques compared to an untreated condition, yielding valuable insights into the impact of these practices. Negative percentage values in all subbasins signify a consistent reduction in surface runoff within the treated watershed.

The variations in the magnitude of reduction across subbasins, ranging from 0% to (−23%), reveal a diverse response to SWC implementation. Significant reductions in values, particularly in the form of larger negative values, were found to be concentrated within the central subbasins of the watershed, which are number 4, 5, 11, 13, and 15. These reductions surpassed 15% when compared to the untreated condition, indicating a considerable impact on the environmental parameters under consideration (Figure 5). For instance, subbasin 3 stands out with a significant difference, showing a 15% decrease in runoff due to the implementation of SWC techniques.

The distinctive feature of these central areas lies in their new cultivation approach, marked by the implementation of olive trees, orchards and wheat crops. In fact, the vegetation cover provided by these crops serves as a natural barrier, reducing the speed of water flow across the soil surface. This helps to minimize the risk of soil erosion and surface water runoff. Furthermore, it’s noteworthy that these zones have undergone thorough treatment employing SWC techniques. This comprehensive approach to SWC signifies a proactive effort to manage soil and water resources effectively. Indeed, the implementation of SWC techniques is crucial for sustaining the health and productivity of the cultivated areas (Figure 5). Overall, the combination of agricultural areas and strategic conservation measures highlights a holistic approach to environmental protection.

On the other hand, the subbasins, namely numbers 1, 2, 8, and 27, exhibit a relatively lower change in comparison, with the recorded change being around 2%. Upon investigating the LULC maps and SWC maps, it becomes apparent that these areas are characterized by low vegetation cover, primarily consisting of rangelands and cereals. Additionally, there is an absence of any form of SWC techniques in these regions.

The combination of these factors provides a realistic explanation for the minimal runoff changes detected in these areas. The low vegetation cover suggests less effective protection against runoff and soil erosion. Without the presence of SWC techniques, there are limited mechanisms in place to mitigate the impact of water runoff. Moreover, the low percentage of runoff changes in these subbasins may also be attributed to a decrease in rainfall values over the years. The combination of low vegetation cover, lack of SWC techniques and reduced rainfall collectively contributes to the observed minimal change in these specific subbasins. Finally, the overall negative values underscore the general mitigating effect of SWC techniques on surface runoff.

The Merguellil watershed demonstrates negative percentage differences in soil erosion, and this favorable trend can be attributed to the effective implementation of SWC practices. These practices have contributed to mitigating soil erosion, leading to a notable improvement in the overall soil conservation status within the watershed.

The observed negative percentage differences reveal varying magnitudes of reduction, ranging from 0% to (−56%). This variability underscores the diverse responses of different subbasins within the Merguellil watershed to the implemented SWC measures. Each subbasin responds differently to conservation practices, shaped by factors like topography, land use, and specific SWC techniques (Figure 6).

The substantial reduction of up to (−56%) in soil erosion in some subbasins highlights the highly effective impact of SWC practices, particularly in areas with vulnerable soil conditions or high erosion risks. Conversely, subbasins showing a 0% reduction may indicate that existing SWC measures are effectively maintaining current soil erosion levels, reflecting a stable conservation state (Figure 6). Overall, the reduction in soil erosion rates, indicated by the negative percentage differences, demonstrates the effectiveness of SWC techniques in promoting sustainable land use practices and safeguarding the Merguellil watershed.

The presence of larger negative values, particularly in subbasins 4, 5, 11, 13, 14, and 15, indicates substantial decreases in soil erosion within these areas. These negative values reflect significant improvements in soil conservation, underscoring the effectiveness of the measures implemented to combat erosion (Figure 6).

The considerable reduction in soil erosion can be attributed to factors such as robust SWC practices, changes in land use patterns, and other conservation initiatives. These combined efforts have contributed to a more resilient and sustainable environment in the mentioned subbasins. This observation is key to understanding the positive impact of conservation measures on soil health and erosion control.

It implies that the implemented strategies in central subbasins have been successful in minimizing the adverse effects of erosion, promoting soil stability, and fostering a healthier ecological balance. Moreover, the relatively minor change in soil erosion is particularly evident in the northwest areas, encompassing subbasins number 1, 2, and 8, as well as in the southern zone represented by subbasin number 27 (Figure 6). Multiple factors contribute to this scenario, including the patterns of LULC, topographical characteristics, and the presence of SWC techniques in these subbasins. Notably, all these areas were classified as untreated zones, indicating the absence of specific interventions.

This outcome emphasizes the crucial role of implemented measures in these specific areas. The limited variation in soil erosion values suggests that SWC measures have been instrumental in maintaining soil stability and effectively controlling erosion.

To enhance the efficacy of future conservation strategies, it is imperative to elucidate the role of land use land cover, topography, and the application of SWC techniques in untreated zones on soil erosion processes. This recognition underscores the importance of continuous efforts to implement effective soil conservation measures to the specific conditions of each geographical region, ensuring sustained environmental health and resilience.

The impact of SWC techniques on groundwater recharge is a critical aspect of watershed management and hydrological sustainability. SWC practices has a central role in influencing the quantity and quality of water that infiltrates the soil and reloads groundwater aquifers. By minimizing soil erosion and controlling runoff, SWC techniques contribute to increased water absorption and percolation into aquifers. The effectiveness of SWC techniques in groundwater recharge can vary based on factors including topography, soil type, climate variables and the implemented SWC techniques.

The modelling of groundwater recharge is of paramount importance in understanding and managing water resources effectively. Indeed, the SWAT model is a powerful tool that provides quantitative insights into the dynamics of water movement within the watershed.

In the Merguellil watershed, the range of groundwater recharge displays noticeable variability, extending from below 10 mm to surpassing 50 mm. The variation observed is likely affected by a combination of meteorological, pedological and anthropogenic factors (Land use practices), as well as the differential effectiveness of SWC techniques. The interplay of these elements contributes to the observed fluctuations in groundwater recharge across the studied area.

Figure 7 explains the percentage of differences in groundwater recharge across various subbasins, revealing distinctive responses influenced by both anthropogenic interventions and natural conditions. The observed variations underscore the importance of water management strategies, considering the unique characteristics of each subbasin. In fact, the differential responses of subbasins to environmental and anthropogenic pressures are influenced by factors such as land use land cover change, altered precipitation regimes, and the implementation of various SWC strategies. Ranging from less than 5% to more than 20%, the differences in groundwater recharge percentages across various subbasins exhibit a diverse pattern. Notably, several subbasins, including subbasins 4, 5, 11, 13 and 15, exhibit an important percentage of groundwater recharge difference (Figure 7). The enhanced infiltration rates into shallow aquifers in central zones can be ascribed to the adoption of innovative SWC practices. Thus, the present study highlights the important role of these techniques in conserving and recharging shallow aquifer in a semi-arid watershed.

Finally, the detection of groundwater recharge in subbasins with established SWC practices highlights the efficacy of these techniques in enhancing water retention and infiltration capacities. This aligns with the main objective of mitigating surface runoff and increasing soil moisture, contributing to sustainable water resource management.

To pinpoint erosion hotspots, we used sediment yield (SY) data simulated by the SWAT model. This data spatially mapped sediment concentration, highlighting areas with the highest erosion rates. To assess the severity of erosion, the annual sediment yield was computed and categorized into different levels, from low to very severe. In fact, the majority of sediment yield originates from cropland areas, with soil erosion severity categorized as follows: (0–5) low, (5–10) moderate, (10–15) high, (15–20) very high, and > 20 t/ha/year as extremely severe. However, the distribution of severity varies across the different sub-watersheds. Through spatial analysis and categorization of watershed areas, we delineated erosion hotspots. These regions, characterized by high erosion susceptibility, require immediate attention for the application of erosion control strategies to mitigate sediment loss and maintain landscape stability.

The result shows that the average annual sediment yield of the study area ranges from <1 to >20 t/ha/year (Figure 8), showing the spatial variations of sediment yield. A comparative analysis of runoff values across sub-basins reveals substantial differences in runoff generation. On the other hand, it is important to notice that the spatial patterns of runoff, soil erosion, and groundwater recharge are significantly influenced by the spatial distribution of soil and water conservation management practices.

Among the 27 sub-watersheds, SW 1, SW 6, SW 7, SW 8, SW 23 and SW 25 presents a very high and extremely severe classes of sediment yield which is >15 t/ha/year (Figure 8). Finally, all sub-watersheds which presented a category from moderate to severe sediment severity classes, are identified as soil erosion hotspot areas. Therefore, the integration of SWC measures becomes crucial in ensuring sustainable water infiltration and groundwater recharge, even in areas naturally predisposed to favorable conditions. Finally, this approach does not only reduce overall sediment yield, but also preserves soil fertility and prevents further land degradation in the watershed. Within the watershed.