Soil Erosion Estimates in Arid Region: A Case Study of the Koutine Catchment, Southeastern Tunisia

Abstract

Soil erosion remains one of the principal environmental problems in arid regions. This study aims to assess and quantify the variability of soil erosion in the Koutine catchment using the RUSLE (Revised Universal Soil Loss Equation) model. The Koutine catchment is located in an arid area in southeastern Tunisia and is characterized by an annual mean precipitation of less than 200 mm. The model was used to examine the influence of topography, extreme rainstorm intensity and soil texture on soil loss. The data used for model validation were obtained from field measurements by monitoring deposited sediment in settlement basins of 25 cisterns (a traditional water harvesting and storage technique) over 4 years, from 2015 to 2018. Results showed that slope is the most controlling factor of soil loss. The average annual soil loss in monitoring sites varies between 0.01 and 12.5 t/ha/y. The storm events inducing the largest soil losses occurred in the upstream part of the Koutine catchment with a maximum value of 7.3 t/ha per event. Soil erosion is highly affected by initial and preceding soil conditions. The RUSLE model reasonably reproduced (R² = 0.81) the spatiotemporal variability of measured soil losses in the study catchment during the observation period. This study revealed the importance of using the cisterns in the data-scarce dry areas as a substitute for the classic soil erosion monitoring fields. Besides, combining modeling of outputs and field measurements could improve our physical understanding of soil erosion processes and their controlling factors in an arid catchment. The study results are beneficial for decision-makers to evaluate the existing soil conservation and water management plans, which can be further adjusted using appropriate soil erosion mitigation options based on scientific evidence.

1. Introduction

Soil erosion is considered one of the most severe environmental problems. Soil erosion degrades soil fertility, soil conservation techniques and urban infrastructure, affecting the land and ecosystems sustainability [1,2]. The situation is being worsened by the increasing occurrence of extreme rainfall events due to the recently varying weather conditions and rising effects of climate change [3]. Accurate quantification of soil loss in a catchment, which is crucial for efficient water and soil management, depends on the knowledge of its controlling factors such as human activities, climate variability, topography and soil properties [4,5]. The soil erosion phenomenon can affect water availability (quantity and quality), agricultural production and sustainable development [6]. Many research activities showed a relationship between the impact of soil erosion and degradation on water quality [7] and the delivery of ecosystems goods and services [8].

In the arid region where water and soil resources are limited, rainfall is known for its high spatiotemporal variability and torrentiality. The driving factors of soil erosion also included population growth, human activities, vegetation degradation, rainfall and climate change, topography and soil erodibility and tectonics movement [9]. Both climate variability and human activity affect runoff and sediment load changes. They are essential factors in explaining the observed variation. However, as the effect of both are integrated at the catchment scale, accurate determination of exactly which part of the observed changes were due to climate variability and which to human activity remains extremely difficult [10].

Soil erosion and land degradation constitute major concerns of different countries in North Africa, such as Tunisia. In this region, the climate is influenced by the Mediterranean Sea to the North and the Sahara desert to the South. The arid, semi-arid and desert bioclimates cover more than two-thirds of the country [11]. The rainfall regime is known for its scarcity, high spatiotemporal variability, torrentiality and poor distribution [12]. Water erosion degrades annually around 15,000 ha of fertile arable soil in Tunisia [13].

The Tunisian arid zones located in the south part of the country cover over 6,290,000 ha, classified as 12% very degraded, 40% medium degraded, 17% marginally degraded and 31% not degraded areas [14]. Jebbari et al. [15] evaluated erosion risk in semi-arid regions based on extreme rainfall characteristics. These kinds of investigations are crucial for decision-makers for better management and designing intervention priorities against soil erosion risks at the catchment scale. For instance, bench terraces have been recommended as effective soil erosion mitigation options by reducing runoff, erosion and sediment transport on the hillslope [16,17].

Soil erosion due to extreme rainfall events is an exceptionally dynamic, distributed and complicated process. It is affected by many factors, among them the topography, vegetation and soil type. Furthermore, soil erosion processes during extreme rainfall events are strongly influenced by the spatiotemporal variability of rainfall characteristics [18]. An experimental study conducted by Escadafal et al. [19] in the south of Tunisia explained that soil erosion is more vigorous when rainfall is very intense, but it does not follow a linear relationship with rainfall rate and area exposed to erosion. Under laboratory flume-controlled conditions, however, it was found that soil loss is proportional to the area exposed during the steady-state equilibrium and when interrill is the dominant erosion mechanism. At the same time, at the beginning of the erosive event, the relationship is less clear due to the pronounced influence of soil’s antecedent conditions [20,21]. Therefore, during heavy storm events, the concentrated rill erosion becomes the dominant soil erosion process, especially under compacted and bare topsoil conditions, resulting in a substantial soil loss. For contour-plowed soil, however, only an exceptional rainfall on wet topsoil conditions causes rill and interrill erosions [19].

In order to attenuate the soil erosion problem, numerous soil and water management options are increasingly implemented. However, the changing climate conditions, including seasonality and spatiotemporal variability, reduce these mitigation measures’ effectiveness. Decision-makers need a detailed understanding of soil erosion drivers, states and impacts in order to better guide mitigation intervention that require thorough spatial soil erosion investigation at the sub-catchment level.

Soil erosion modeling has shown to be a beneficial tool for predicting and describing soil erosion processes and evaluating soil losses and their controlling factors [1]. In arid regions, lack of systematic monitoring and data sharing are commonly limiting factors for modeling use. For instance, installing permanent monitoring devices for soil erosion in a typical dry area is a major challenge due to weather conditions, vandalism and technical problems (such as lack of power, appropriate maintenance and presence of dust). Therefore, numerous studies have considered rainfall simulation and/or existing natural or human-made erosion collection structures (such as cisterns, buffer strips, etc.) as appropriate approaches to investigate soil erosion processes and yields. One of the first models used was the USLE (Universal Soil Loss Equation) model introduced by the US Department of Agriculture (USDA) [22]. Later, it was revised [23] and retitled to RUSLE (Revised Universal Soil Loss Equation), which produces a spatial distribution of soil erosion and could be applied by decision-makers for conservation and management processes [24]. Northern regions of Tunisia have benefited from some studies such as [15,25,26,27,28,29,30]. RUSLE model, adapted to the semi-arid conditions in central Tunisia, provided good estimation results of measured soil erosion delivery [28,31]. However, using this model, few studies are conducted in the arid regions characterized by the short and violent character of their rainfall events and have very remarkable erosion yields.

This study is one step in that direction, where the spatial variability of soil loss in the Koutine catchment has been intensively investigated for 4 years. This study presents a soil erosion assessment methodology in data-scarce and dry areas combining field measurement and modeling. The objectives of this study are to (i) assess the spatial variability of soil erosion in the Koutine catchment and its controlling factors, (ii) quantify the soil losses using intensive field campaigns and laboratory analysis and (iii) to evaluate the RUSLE soil erosion model using the in situ data. The validated model provides the basis upon which decision-makers can test the efficiency of different soil erosion mitigation measures using a scientifically based tool.

2. Study Site

The study site (Wadi Koutine) is located in the Jeffara area, northeast of Medenine city, which is part of the southeastern Tunisia region (Figure 1). The study site covers an area of 279 km², drained by Wadi Koutine, which starts from the mountain range of Béni Khedache and reaches the city of Koutine. The approximate coordinates of the central point of the Koutine catchment are 33°16′ N and 10°08′ E. The highest elevation of the Koutine catchment is about 650 m that corresponds to the peak of mountain Jbel Moggar. The climate in the study site is arid. The coldest months are December–February. June–August is the warmest period of the year, during which the temperatures could reach as high as 48 °C. The temperature is affected by both the proximity to the sea and the altitude. The rainfall is characterized by low mean annual values, high spatiotemporal variability and torrential features [32]. Annually, the Koutine catchment receives mean rainfall varying between 234 and 162 mm in the mountains and eastern parts, respectively [33].

The annual average runoff of Wadi Koutine has been estimated to be 8.74 Mm³/yr [34].

Due to the low annual rainfall rates and their spatiotemporal variability, together with the relatively high elevation of the Koutine catchment, different surface water harvesting techniques were implemented (i.e., Jessour and Tabia). The Jessour are mainly constructed in the mountainous areas of the catchment. The Jessour is an ancient water-harvesting technique that is constructed in the inter-mountain and hill watercourses to intercept runoff and sediments [35]. This hydraulic unit is composed of three main components: a dike in the form of a small earth embankment with a spillway made of stones, a terrace that represents the cropping area and an impluvium, which is the runoff catchment area [36]. Tabias are basically situated in the middle of the catchment on moderate slopes. A Tabia is formed by a principal embankment of 50 to 150 m long situated along the contour with lateral bunds of about 30 m long at the ends. The Tabia gains its water directly from its impluvium or by the diversion of Wadi runoff. Water is captured until it reaches a height of 20 to 30 cm, after which it is diverted, either by a spillway or at the upper ends of the lateral bunds [37].

In the areas behind Tabia and Jessour, water-harvesting techniques have high effectiveness in combatting water erosion. However, in the rangelands, current conservation measures have low effectiveness, and the rate of land degradation is increasing [38]. The Koutine catchment is used for cropland, including the production of cereals and grazing in various types of rangeland. The land use is dominated by sparsely covered, degraded steppes. Cropped sites, mainly for growing olives, are found on terraces behind water-harvesting structures [35].

3. Materials and Methods

3.1. Application of the RUSLE Model for Soil Erosion

Soil erosion modeling has shown to be a beneficial tool to reasonably describe and predict measured data under appropriate assumptions [1]. Empirical modeling, in particular, has proven to represent soil loss variabilities at the catchment scale considering simplified cause–effect relationships. Even though empirical modeling is a data-demanding approach and is usually implemented in rich in situ data studies, it is regarded as a beneficial approach when parametric models cannot be constructed and study site observations are available [39]. Empirical modeling allows the user to understand the relationships between the drivers and state variables that belong to a particular system or process based on different methodologies and approaches. In soil erosion, empirical models can support decision-makers to take the best solutions and select suitable scenarios for development plans.

The RUSLE soil erosion quantification and assessment model is an empirical spatially distributed model based on multiplying five independent factors (Figure 2). As a result, the total annual soil loss rate (A) is calculated:

$$A=R\times K\times C\times P\times LS$$

(1)

where A is the annual rate of soil loss (t/ha/yr), R is rainfall erosivity factor, K is soil erodibility factor, C is land cover management factor, P is practice support factor and LS is topography factor.

The use of GIS and erosion model has been integrated to estimate the degree and spatial distribution of erosion in the study area through mapping of the five different soil loss factors (R, LS, C, K and P).

3.2. Rainfall Data and Erosivity (R)

Rainfall erosivity (MJ mm/ha/h/yr) is one of the major factors that control soil erosion in Jeffara, south of Tunisia. Several indices have been developed in the past to quantify rainfall erosivity. In this study, we examined the performance of numerous indices with different temporal-resolution of the required rainfall inputs data. In Koutine catchment, rainfall data were measured with monthly and daily rainfall data were used to calculate the Fournier Index (FI) and half-month erosivity index (Mi).

Fournier developed an erosivity index using monthly rainfall data [40]. The FI is usually used due to the availability of monthly rainfall (Equation (2)). This index was later modified and called the Modified Fournier Index (MFI). It is defined by Equation (3). Another modification form of the Fournier index has been described in Equation (4) [41].

$$FI=\frac{P{i}^2}{P}$$

(2)

$$MFI=\frac{{{\displaystyle \sum }}_1^{12}P{i}^2}{P}$$

(3)

$$MFI=\frac{P_{wettest}^2}{P}$$

(4)

where Pi is monthly rainfall, P is annual rainfall and P_wettest is the maximum monthly rainfall (wettest month).

In our study, the R factor was calculated using Equation (5), where a and b are regional coefficients. For Koutine catchment, a and b were taken as 1 and 0, respectively. For MFI calculation, Equation (3) was used due to the data availability of monthly and annual rainfall data [42].

$$R=a\times MFI+b$$

(5)

Given the arid climate, the rainfall in the study site is characterized by low annual mean, high irregularity (both in time and in space) and torrential form. The annual average rainfall during observed years (1969–2018) is between 234 mm in the mountains and 162 mm in the downstream.

The study area has been equipped with climatic and weather stations for several years (

Table 1. Weather and rain-gauge stations and their corresponding data collection periods. P refers to the mean annual rainfall values generated from daily data.

Station	X (m)	Y (m)	Period	P	Type	Source *
Béni Khedache	611,653	3,679,340	1969–2018	234.2	Weather station	IRA
Ksar Hallouf	606,756	3,682,374	1969–2017	178.9	Rain gauge	CRDA
Ksar Jedid	620,637	3,685,121	1969–2017	172.8	Rain gauge	CRDA
Toujane Dkhila	609,157	3,699,761	1969–2017	193.1	Rain gauge	CRDA
Alamat	623,390	3,695,498	1969–2017	162.2	Rain gauge	CRDA
Koutine	627,980	3,701,600	1969–2017	173.1	Rain gauge	CRDA
Medenine	638,440	3,689,600	1969–2018	169.9	Weather station	INM
Sidi Makhlouf	637,412	3,707,889	1969–2017	179.1	Rain gauge	CRDA

* IRA: Instiut des Régions Arides, INM: National Institute of Meteorology, CRDA: Regional department of the Ministry of Agriculture.

). These data acquired in situ are being updated, depending on the availability of existing data at different time steps, including daily, monthly or annual. Initial quality analysis of rainfall time series focused on their duration and quantity in mm. Rainfall data were collected for more than 48 years (1969–2018) with a variable time interval (daily, monthly and annual) from eight stations located in the study region (Table 1).

3.3. Soil Characteristics and Erodibility (K)

Since the 1930s, the majority of studies have defined soil erodibility according to inherent soil properties. They focused on the role of soil texture, chemical properties such as soil organic matter and soil profile descriptors such as structure and permeability [43]. The factor K expresses the vulnerability of the soil to be eroded by rain. It depends on the physical and chemical properties of the soil (such as particle size distribution, aggregation, structural stability, porosity, organic matter content and soil moisture).

The soil texture is determined in different sites of the catchment by the mass percentage of sand, silt and clay particles in selected soil samples. The soil particle size distribution determines whether it is sand, silt (loam) or clay. Sand particles have a diameter of 0.05 mm to 2 mm, the diameter of a silt particle ranges from 0.002 mm to 0.05 mm and the diameter of clay particles is smaller than 0.002 mm. Wischmeier and Smith [22] and Renard et al. [23] recommend an algebraic estimation of the soil permeability that includes five soil parameters (texture, organic matter, coarse fragments, structure and permeability) in the below equation (Equation (6)):

$$K=2.173\times (2.1\times {10}^{-4 }{M}^{1.14}\left(12-OM\right)+3.25\left(s-2\right)+2.5\left(p-3\right))/100$$

(6)

where K is erodibility (t ha h/ha/MJ/mm), M is textural factor, which is equal to (m_silt + m_vfs) × (100 − m_c), m_c (%) is clay fraction content (<0.002 mm), m_silt (%) is silt fraction content (0.002–0.05 mm), m_vfs (%) is very fine sand fraction content (0.05–0.1 mm), OM (%) is organic matter content, s is soil structure class (s = 1 for very fine granular, s = 2 for fine granular, s = 3 for medium or coarse granular and s = 4 for blocky, platy or massive), p is permeability class (p = 1 for very rapid and p = 6 for very slow). Six classes are distinguished [44]: (1) rapid, (2) moderate to rapid, (3) moderate, (4) moderate to slow, (5) slow and (6) very slow.

In this study, the soil map (Figure 3), based on the French soil classification system (CPCS) [45], was used to locate and collect soil samples from Wadi Koutine catchment (Figure 4). The physicochemical analyses were determined in the soil laboratory of Arid Regions Institute (IRA) of Médenine. For each sample, the particle size analysis was performed in three classes (sand, silt and clay) using Robinson’s pipette method [46]. Organic carbon was quantified by the Walkley–Black wet dichromate oxidation [47] and converted to organic matter by multiplying it by 1.724. Soil permeability (hydraulic conductivity) was determined using the flow rate per unit area, measured perpendicular to the direction of flow, at a constant gradient of hydraulic pressure and soil texture was obtained using Soil Textural Pyramid [48].

According to the Taamallah map [49], Regosols and Rendzinas (Figure 3) cover most parts of Koutine catchment. Soil types in the mountain area are mainly Regosols. These soils are generally thin to fully infiltrate surface runoff and therefore are easily eroded. Soil characteristics of the Koutine catchment are distinguished by three surface formations [49] (Figure 3):

• The coherent formations that include the calcareous crust and gypsum encrustation frequently exist in the Tunisian Jeffara;
• The remaining formations consist of loess and sandy loams that occur in the western mountain area;
• The formations of red silts cover underlying formations, particularly the calcareous crust is encountered mainly in the littoral zone.

The physical soil characteristics (organic matter and soil texture) are acquired by laboratory analyses of sixty-two (62) samples collected in the field from different soil types and with homogeny spatial distribution at the scale of Koutine catchment.

3.4. Land Cover Management Factor (C)

The value of C factor depends on plants characteristics related to vegetation type, stage of growth and cover percentage [50]. For the Tunisian context, the C value has been defined for different plant covers [27]. In RUSLE model, C factor means the ability of vegetation cover to reduce soil erosion. It ranges from zero (total protection) to one (no protection). Values are assigned from RUSLE tables or others based on plots experimental studies. In this study, we performed a supervised classification of the 30 m resolution Landsat image (downloadable from USGS Open Access Hub (https://earthexplorer.usgs.gov/, last accessed on 28 May 2021) to extract a land-use map. C values were then assigned according to the tables used in regional and international studies [29,30].

3.5. Practice Support Factor (P)

Practice support factor represents the role of different conservation practices in reducing soil erosion. Its values vary from zero (total protection) to one (no protection). In our study area, we proceeded to generate the percentage slope map on which we assigned the values, by slope class, of the conservation practices most used in the Koutine catchment, namely, Contour, Jessour and Tabias. We applied the values cited in SHIN [51] and used them for predicting water erosion in arid Tunisian lands [52].

3.6. Topographic Factor (LS)

Several studies are interested in the effect of topography on erosion. Borst et al. [53] and Wischmeier et al. [22] exposed that topography acts on soil loss by two compounds such as length and degree of the slope. The topography also intervenes by the slope profile shape, depending on whether the slope is concave or convex. On the convex slopes, erosion is very low upstream, but it increases rapidly downstream. On concave slopes, erosion is very high upstream but decrease gradually downstream [54]. The slope can similarly act on erosion through infiltration parameters. With the increasing slope angle, Wang and al. [55] observed a decrease in the runoff rate attributed to the destruction of the crust of the crop under the kinetic energy action of runoff, which increases with the slope.

LS factor means the role of the slope length (L) and its steepness (S) in soil erosion. To obtain the LS value map, we chose to use the Unit Stream Power Erosion and Deposition (USPED) method for calculating the LS factor because it was obvious that it could be completed with the spatial analysis tools included in a standard ArcMap installation [56]. The USPED is a physically established model that integrates a spatial module. In RUSLE, $L$ is dependent on linear distance ${\lambda }_i$, which is the horizontal length from the start of sediment transport to point $i$ on the slope. Thus, they are inherently a single-dimensional function. The USPED uses the area of upland as an alternative, contributing flow at distance $i$. The area is switched in place of the former slope length. The L calculation for point $i$ on a slope is shown in Equation (7).

$$L=\left(m+1\right){\left(\frac{{\lambda }_A}{22.1}\right)}^m$$

(7)

where $L$ is the slope length factor at some points on the study site landscape, ${\lambda }_A$ is the region of upland flow, $m$ is an adaptable value depending on the soil’s vulnerability to erosion and $22.1$ is the unit plot length.

The S calculation is shown in Equation (8).

$$S={\left(\frac{\sin \left(0.01745\times {\theta }_{deg}\right)}{0.09}\right)}^n$$

(8)

where $\theta$ is the slope value in degrees, $0.09$ is the value of the slope gradient constant and $n$ is an adaptable value reliant on the soil’s vulnerability to erosion.

In our case study, we used 12 m digital elevation model, which is optimal for this kind of work [52].

3.7. Field Monitoring

A total of 25 cisterns (6 Fesguia and 19 Majel) have been selected for erosion measurement by monitoring solid deposits after rainfall events (Figure 5). In the context of RUSLE simulation, calibration and validation were conducted during the assessment period. Erosion was estimated by measurement of deposited sediment in selected cisterns sediment settlement basins in the Koutine catchment. Cisterns are often built to catch and store rainwater. Fesguia and Majel are two common types of cisterns. Each unit consists of three main parts: the impluvium, the sediment settlement basin and the storage reservoir (Figure 6). They are distinguished according to the form of their reservoirs. Generally, Fesguia has a larger storage capacity, and therefore its construction cost is relatively higher than Majel. In our case study, the measurement concerned only deposited sediment in the settlement basin of each infrastructure (Fesguia and Majel), which depends only on the topography, soil and land cover of the impluvium. Therefore, the differences between Fesguia and Majel have no influence on the measured deposited sediment. Majel is commonly implemented in the upstream area (mountains), while Fesguia is generally constructed in the downstream area of the catchment characterized by a gentle slope (plain).

The measurement sites were defined by three Cartesian coordinates (x, y, z) in UTM projection Zone 32N (Figure 4). During field campaigns, we measured the thickness and the weight of the solid deposit in the sediment settlement basin. Before cleaning sediments for the next observation time, a sample is taken for laboratory analysis (soil humidity, organic matter). Sediments data for the entire network are collected and reported in the measurement sheet. These ratios account for every cistern and characterize the sediments. The eroded soils are transported towards the cisterns sediment settlement basins. According to field measurements, it is important to mention that erosion values were significant only for precipitation events that exceed 20 mm [57].

Sediments in cisterns settlement basins represent an important record of relatives’ sub-catchments erosion. In fact, they give information about the main exceptional rainfall events responsible for erosion delivery. From this point of view, the measurements were used to estimate average and maximum erosion rates for each sub-catchment. The average erosion rate was calculated at the end of the year using measurements after rainfall events of solid deposits. From these data, an annual average erosion (t/ha) was determined. These calculations were executed for all twenty-five (25) cisterns sediment settlement basins.

4. Results and Discussion

4.1. Rainfall Data and Erosivity (R)

The analysis and evaluation of the rainfall data trends were conducted based on the spatial and inter-annual variation coefficients. The selection of different stations was determined based on the quality and availability of rainfall data. R factor calculated for the eight rain-gauge and weather stations around Koutine catchment shows, after interpolation, a gradient of values evolving between 24.97 and 33.15 MJ · mm/ha/h/yr (Figure 7a).

During the monitoring period, the spatial distribution of the rainfall erosivity in the year 2017 shows that the Koutine catchment is subjected to moderate rainfall intensity upstream. The MFI was equal to 109.0, 106.8 and 113.0, respectively, in Béni Khedache, Zammour and Allamet stations. The R factor values used in our study were adopted from Kefi et al. [42], considering the different soil properties and rainfall characteristics of both studies. This results in different R factor ranges (R is in [24–33] in our study <[31–66] used by Kefi et al. [42]).

4.2. Soil Erodibility (K)

The K factor varies in the Koutine catchment between 0.031 and 0.075 (Figure 7b). According to the USDA triangle, soil types can be classified into 12. The result shows that 45.2% of samples had Sandy Loam (SL) texture, 33.9% had Loam Sandy (LS) texture, 19.4% are Sandy (S) and 1.6% are Sandy Clay Loam (SCL). The study area had a high sand percentage that varied from 54.7% up to 94.4%. In contrast, the clay content ranged from 3.7% up to 32.2% and silt content from 0.0% up to 29.7%. Detailed information on soil properties of the Koutine catchment and impluvium of different cisterns are given as Supplementary Materials.

Erodibility is low for soils with a high percentage of clay because clay particles aggregate together and resist detachment and transport processes [58]. According to Equation (6), it was found that erodibility is negatively correlated with clay content (Figure 8a) but is positively correlated with sand content. This result corroborates with the finding of Duiker et al. [59].

The results of this analysis showed that organic matter content varied from the minimum of 0.18% up to a maximum of 2.02%, with 60% of samples ranging from 1.01% up to 2.02%. Organic matter content is considered a moderate level that can prevent soil from detachment by the rainfall kinetic energy and provide a shallow threshold of soil erodibility.

According to USDA classification, 96.0% of Wadi Koutine catchment have soils with moderate to rapid permeability. It was observed in numerous studies [60] that the presence of clay decreases the level of soil permeability.

4.3. Topographic Factor (LS)

The average elevation of Béni Khedache Mountain is 650 m, where the highest point is at Jbel Moggar, with an elevation of 700 m (Figure 7c). High precipitation occurs in the upstream and middle parts of the Koutine catchment. Due to the topography of the region, the direction flow is West–East. The average slopes in sub-catchments of selected cisterns were more than 10% in the upstream part, with a maximum of about 30%, while for the middle and downstream areas, the average slopes are 6 and 2%, respectively.

The LS factor in the study site shows values in the range 0–227.3. A total of 30% of the catchment has LS values less than 1. About 14% of the lands are between 1 and 5. Land with LS values between 5 and 10 represents 53%. Only 3% of the catchment shows values exceed 10 (Figure 7d).

4.4. Land Cover Management Factor (C)

In the Koutine catchment, the C factor shows four classes of values: Rangeland (C = 0.15) represents 39% of the study site. The Bare land (C = 1) represents 22% of the catchment. The land covered by olive trees and fruit trees behind Jessour and Tabias (C = 0.9) occupies 13% of the total surface. Intensively cultivated plots (cereals, vegetables, drop by drop irrigation, legumes, etc.) or well-preserved tufts of natural vegetation (C = 0.4) cover 26% (

Table 2. Land cover management values for different C factor classes [30].

Vegetation Cover	C Factor
Bare soil, fallow	1.0
Fruit trees	0.9
Cereal rotation	0.4
Ordinary pasture	0.15

, Figure 7e).

4.5. Practice Support Factor (P)

Different water and soil conservation (WSC) practices were essentially implemented in the study site, such as Jessour, Tabias, Check dams and Contour benches. The efficiency of these practices varies according to the slope value (

Table 3. P factor values for different slope classes in Koutine catchment [51].

Slope Classes	Contouring	Terracing
0.0–7.0	0.55	0.12
7.0–11.3	0.60	0.14
11.3–17.6	0.80	0.14
17.6–26.8	0.90	0.16
>26.8	1.00	0.18

, Figure 7f). The WSC practices are entirely useless on slopes higher than 26% (p = 1), representing only 3.5% of the Koutine catchment. They become less efficient as the slope declines.

4.6. Simulated Annual Soil Loss (t/ha)

The simulated annual soil loss results of the RUSLE erosion model in the Koutine catchment vary between 0 and 70.9 (t/ha/yr) with an average annual rate of 1 t/ha/yr (Figure 9) and total losses of 26,487 t/yr. The distribution under RUSLE classes of soil loss shows a predominance of loss lower than 0.5 t/ha/yr, which concerns the major part of the Koutine catchment (70%). Then comes the class 1–5 t/ha/yr, which characterizes 25% of the catchment and corresponds approximately to the piedmonts of the reliefs and the moderate slopes. Rates above 5 t/ha/yr represent only 5% of the study site. This means that 95% of the land in the catchment has less than 5 t/ha/yr (Figure 9).

In order to validate simulated results, in situ soil loss measurements were quantified based on observation of deposited sediment in the sediment settlement basin of cisterns mentioned above.

Table 4. Cisterns characteristics and yearly soil loss during the sampling periods (2015–2018). The fi and mi referred to Fesguia and Majel, respectively. Six erosive events (put in bold) are detailed in Figure 10 below.

Sampling Number	Code	Sub-Catchment Area (m2)	Slope (%)	Observation Years	Number of Storm Events	Mean Soil Loss (t/ha/yr)
1	f1	32,500.00	2	2015, 2016, 2017, 2018	9	0.016
2	f2	2528.82	2	2015, 2016, 2017, 2018	13	0.822
3	f3	5923.79	2	2015, 2016	6	0.023
4	f4	832.56	2	2015, 2016, 2017, 2018	9	0.258
5	f5	1092.55	2	2017, 2018	6	0.329
6	f6	963.54	3	2017, 2018	6	0.156
7	m1	368.94	3	2015, 2016, 2017, 2018	9	0.065
8	m2	111.51	6	2015, 2016, 2017, 2018	13	0.459
9	m3	304.67	7	2015, 2016, 2017, 2018	13	0.559
10	m4	744.34	10	2015, 2016, 2017, 2018	13	0.277
11	m5	233.68	11	2015, 2016, 2017, 2018	13	1.348
12	m6	233.72	14	2015, 2016, 2017, 2018	13	1.754
13	m7	163.83	15	2015, 2016, 2017, 2018	13	0.286
14	m8	1062.53	3	2017, 2018	6	0.566
15	m9	2970.05	5	2017, 2018	6	0.277
16	m10	4416.08	5	2017, 2018	6	3.071
17	m11	504.36	6	2017, 2018	6	0.031
18	m12	312.68	3	2015, 2016, 2017, 2018	9	0.158
19	m13	1213.62	5	2015, 2016, 2017, 2018	9	0.031
20	m14	232.41	30	2015, 2016, 2017, 2018	13	2.711
21	m15	740.00	5	2015, 2016, 2017, 2018	9	0.732
22	m16	1591.33	5	2015, 2016,	5	0.227
23	m17	400.22	2	2015, 2016, 2017, 2018	9	0.049
24	m18	530.00	17	2015, 2016, 2017, 2018	13	12.504
25	m19	947.00	22	2015, 2016, 2017, 2018	13	9.008

shows a summary of cisterns characteristics and mean annual soil loss during the monitoring period. Measurement campaigns started in 2015 and continued for a 4-year period. The number of observations for each cistern stretches between 5 and 13 depending on the characteristic of the storm event and the importance of the erosive event. It is worth mentioning that the deposited sediment was not monitored for all storm events. The measurements from storm events characterized by very short duration and negligible eroded sediment were omitted during the study. Therefore, the number of observations depended on cisterns location and storm events; a total of 240 measurements were obtained (Table 4).

The annual soil loss in the Koutine catchment during the monitoring period stretches between 0.01 and 12.5 t/ha/yr, slightly exceeding the permissible limit of 12 t/ha/yr [22]. Fersi et al. [57] also recorded similar results during an observed period of 5 years (1972–1977) with a minimum erosion of 0.3 t/ha/yr and a maximum of 12.3 t/ha/yr. A maximum of 7.3 t/ha per event was observed after a specific erosive event occurred in the upstream part of the Koutine catchment. Simulation results obtained by Schiettecatte et al. [61] using modified Sediment Transport Model (2D-version) indicate that values stretch between 1 and 5 t/ha/yr.

The maximum soil loss rate was observed in the two successive relatively wet hydrological years (2016/2017 and 2017/2018) for specific cisterns. However, for some cisterns, there were only data for single- or two-measurement campaigns due to spatial rainfall distribution. This was the case for cisterns f3, m9, m13, m15 and m16, for which the maximum solid deposits in sediment settlement basins could not be determined.

At the end of every year during the monitoring period (2015–2018), the annual soil erosion was evaluated at selected measurements points. The erosive rainfall events that characterize the different erosion periods are significant in 2017 and 2018. Therefore, the dates related to the above measurement periods (for average and maximum soil loss) were used to identify the exceptional rainfall events that occurred during these periods.

The spatial distribution of erosive rainfall together with topography soil texture and organic matter were analyzed for 25 sub-catchments. Since there were climatic (spatial rainfall distribution), physical (topography, slope), soil characteristics (texture, organic matter) and geographical conditions, the catchments were divided into three groups. The first group is located in the catchment upstream, represented by 10 cisterns (9 Majel and 1 Fesguia) and characterized by a mean annual rainfall of about 200 mm, and the mean slope exceeds 15%. The second group is located in the middle of the Koutine catchment and is represented by eight cisterns (6 Majel and 2 Fesguia), characterized by a moderate slope (4%) and a mean annual rainfall of 183 mm. The last group contains seven cisterns (4 Majel and 3 Fesguia), and the mean slope is 3% and the annual mean rainfall is 168 mm.

The mean annual soil loss was 3.8 t/ha in the upstream and 1.3 t/ha in the middle and downstream regions. Figure 10 emphasizes the relationship between soil loss and slope. Results showed that the permissible limit of 12 t/ha/yr [22] was achieved in the upstream under the slope of 22%. In the downstream area, values stretch only between 50 kg/ha and 227 kg/ha (Figure 10).

Figure 11 represents the relationship between specific events and maximum soil loss during a wet year (2017). Soil erosion loss was more important for specific events with high precipitation values (over 40 mm). The minimum soil erosion was 0.2 t/ha for events with a depth between 20 and 38 mm.

The maximum (3.7 t/ha) was recorded for a rainfall event with an intensity depth of only 64.6 mm. Thus, erosion does not depend only on precipitation value but also on its intensity and soil properties. The event with a rain depth of 173 mm that happens at the beginning of the dry season results in less soil loss than a rainfall event of 64.6 mm that occurred in the wet season (Figure 11), reflecting the effects of antecedent conditions on soil-erosion delivery.

At the end of this work, and to verify the accuracy of the model results, we made a comparison between values of erosion predicted by the RUSLE model and observed results (Figure 12).

The RUSLE model could reasonably reproduce the measured soil losses (Figure 12), reflected by the square of the correlation coefficient (R²) of 0.81. Our results showed lower soil loss in the arid region of Tunisia (average of about 1 t/ha/yr). Other studies revealed higher soil erosion in the semi-arid compared to the arid regions (such as [15,42]). In the arid region, the topography and the erosivity are the determinant factors for soil-erosion delivery. While in the semi-arid region, the interaction between rainfall and land use becomes an important factor in either accelerating or attenuating the soil detachment during the erosive events [62].

Comparing the spatially distributed RUSLE model simulation and observed values contributes to a better understanding of the mechanisms and processes leading to catchment soil loss. The validated model using appropriate meteorological and hydrological data can be used as a predictive tool for better describing and predicting soil erosion processes and deliveries under different scenarios. In addition, the RUSLE model can serve as a decision-making tool for testing the efficiency of different soil erosion mitigation options using a scenario-based approach. Gained insights from combining in situ data and appropriate modeling investigation are fcrucial for successfully implementing science-based management.

5. Conclusions

Under the current climate variability, the study site is characterized by arid conditions and extreme inter-annual and spatial variability in total rainfall. This study conducted in the Koutine catchment showed that the arid regions of Tunisia, despite their low and variability annual rainfall, could be subject to significant rates of water erosion. The use of finer DEM resolution (12 m) as a basis of this work allows an appropriate and even an optimal level of precision for the quantification of erosion at a local scale (pixel level). We opted for the use of the modified Fournier index to determine R factor values. The superposition of the different model layers allowed us to compute the different parameters and obtain an annual rate stretching between 0 and 70 t/ha/yr with an average rate of 1 t/ha/yr. Obtained results were validated by the field measurements. The average annual soil loss in selected sites in the Koutine catchment during the monitoring period varies between 0.01 and 12.5 t/ha/yr. A maximum of 7.3 t/ha per event was registered in the upstream after a specific event. The importance of soil loss depends on the slope, and it was more than the permissible limit of 12 t/ha/yr in the upstream and 0.227 t/ha in the downstream area. It was also shown that cisterns could be used as effective soil erosion monitoring sites in dry regions. Obtained results can be used for planning purposes to develop mitigation strategies for soil erosion management. In addition, the results can be used to design water-harvesting techniques, especially in the upstream area, and eventually to decrease soil erosion.

The developed RUSLE model is quite reliable and can be continuously updated based on future available data. Additionally, it can serve as a science-based decision supporting tool which can be extrapolated to the regional level for sustainable development of soil and water management plans.