Modeling of Soil Water Erosion in Tunisia Using Geospatial Data and Integrated Approach of RUSLE and GIS

1. Introduction

Large-scale soil erosion was known to be one of the most severe problems that can lead to environmental damage [1,2] and consequences affecting the political, social and economic aspects of countries [3], particularly the developing countries. The most dominant agent of the soil erosion is the water [4], removing the soil surface materials and including detachment, transportation and deposition of the particles by rainfall and runoff. Soil erosion caused by water is defined as the breakdown of soil structure and detachment of soil particles due to falling raindrops and water flow exceeding a critical threshold [5]. The raindrops, with which soil particles (sediments) are detached, hit the topsoil of ground surface by splash during rainfall. The detached sediments can be transported to rills. The rills were gradually joined together to form large channels leading to gully erosion. The natural process of the removal of soil particles by water runoff and their redistribution downslope represents one of the main types of soil degradation [6]. It can cause negative consequences on soil productivity, human health and the earth’s environment, as well as on habitat which can have gradual degradation [7,8]. Many researchers have focused their studies on soil erosion by water to assess soil erosion losses in different regions of the world using remote sensing and GIS technology [5,9,10,11,12,13]. The Land Use and Land Cover (LULC) can influence the soil erosion process in a positive or negative way [14]. Human activities can accelerate the process of erosion, transport, and sedimentation such as urbanization, mining and construction of roads, highways and dams [8].

More than 80 soil erosion models [15,16,17,18,19,20,21,22,23], with varying degrees of complexity, have been developed to evaluate potential soil loss for different spatial and temporal scales [24,25], such as the Universal Soil Loss Equation (USLE), which was developed by Wischmeier and Smith in 1965 from measurements on elementary experimental plots respecting exact dimensions [26]. However, this model had certain limitations since it only took into account sheet erosion processes at the scale of the plot [27,28,29]. To remove these limitations, the model has been improved, modified and revised in several versions by integrating runoff factors to adapt it to the scale of a rainfall event [30]. Indeed, the USLE model has moved to MUSLE model "Modified Universal Soil Loss Equation", taking into account the topographic complexity through the use of the Digital Elevation Model (DEM) and anti-erosive practices. Finally, the RUSLE “Revised Universal Soil Loss Equation” improves the determination of the different factors of soil water erosion [31]. This new Revised USLE equation maintains the basic structure of USLE but uses new algorithms to calculate and estimate the individual factors. RUSLE model takes the form of a mathematical equation that uses erosion factors as inputs to estimate the mean annual soil losses resulting from sheet and rill erosion [32]. It is an empirical model that brings together the factors affecting the rate of water erosion, namely the kinetic energy of intense rainfall, soil properties, terrain characteristics, soil protection by vegetation cover and anthropogenic practices. This model does not consider erosion processes such as detachment, transport and deposition to estimate soil loss. The revised Wischmeier equation is combined with GIS techniques to assess soil loss rate and its spatial distribution across different land covers, taking into account the advantages offered by GIS and geospatial data [33,34,35] that can lead to more efficient, more precise and less time-consuming results [36], thus helping to guide decision-making. RUSLE model is not used to estimate the amount of sediment migrating from a specific watershed, but the amount of soil lost from any area [37]. For this, the model retains the same form as the equation used in the USLE model.

The RUSLE model involves the spatial combination of the different factors contributing to soil erosion. It was used to calculate soil losses (A) which is a multiplicative function taking into account the erosivity and the aggressiveness of rainfall and runoff (R factor) [megajoules millimeter hour−1 hectare−1 year−1], the soil erodibility (K factor) [ton hour megajoules− 1 millimeter− 1] and the resistance of the environment (C, P and LS factors) [dimensionless]. LS in RUSLE equation is generally the combination of L and S, representing the effect of the topography on erosion rates [38]. The RUSLE model has been used worldwide and adapted according to the climatic, topographic and soil context [9,10,39,40]. It is still used widely to estimate soil erosion all over the world [41] and can be used to estimate soil erosion for large areas, up to country level [42].

Therefore, modeling of soil water erosion in Tunisia using geospatial data and integrated approach of RUSLE and GIS aimed to estimate the spatial distribution of soil erosion of the overall study area, which can be useful to contribute to a better understanding of the soil degradation of large-scale areas. The spatial soil erosion map produced, can serve as a useful input for deriving land planning and management strategies and provide an opportunity to develop a decision plan for soil erosion prevention for better control and protection of both natural and man-made resources available in Tunisia.

2. Materials and Methods

2.1. Description of the study area

The republic of Tunisia, located in the North of African continent, is situated between 30°13’ and 37°20’ north latitude and 7°32’ and 11°36’ east longitude (Figure 1). It is located also in the south part of the Mediterranean Sea and bordered by Lybia at the South-East and by Algeria at the West. The study area covers 155084 km², including the two main islands (Jerba and Kerkennah) but without the other small islands, such as Kneiss, Zembra, Kuriat, Galite and Chikly. Tunisia is the smallest country in North Africa, with a coastline on the North and East having around 1,300 km.

Tunisia has different types of landscapes, with mountainous areas in the North-West (Figure 2), upper and lower Steppe areas in the center and wide plains in the East. The depressions of the great Chotts mark the beginning of the Sahara in the South, with the mountains of Dhahar and the plains of Jeffara. The lowest altitude is 27 m below sea level, relating to the chott areas, and the highest is 1549 m above sea level, in Jbel el Chaambi, with an average elevation of 256 m for the whole study area.

Tunisia has ten different bioclimatic zones (Figure 3), characterized by three factors: the mean annual precipitation which decreases from North to South, the highest mean temperature of the hottest month (August) and the mean lowest temperature of the coldest month (January) [43]. The winter in Tunisia is soft and humid with temperatures between 8°C and 15°C and the summer is hot and dry with temperatures between 22°C and 35°C which can even exceed 40°C in August. Its northern part has a humid higher, humid lower and sub-humid bioclimate with an average annual rainfall, for the previous 31 years (from 1990 to 2020), of 650-736 mm. The central regions have a semi-arid higher, semi-arid middle, semi-arid lower and arid-higher bioclimate with 250-650 mm annual rainfall. The southern part has an arid lower, Saharan higher and Saharan lower bioclimate with 37-250 mm annual rainfall.

2.2. Methodology and data source processing

The cartographic documents are generally produced at local scales. However, these data are often without recent updating, incomplete and sometimes non-existent, such as geographical reference data on topography, precipitation, land cover, etc... Similarly, georeferenced digital data on a small scale are mostly scarce or obsolete. On the other hand, the use of online databases on a global scale has become available these last years. Furthermore, the spatial analysis and modeling are facilitated mainly because of the existence of satellite data which has become very easy to have continuously. Indeed, the raster mode with a simple data structure and a regular shape of the grid is very convenient to keep geometric properties common to all layers of information, thus facilitating the combination of different layers, since all numeric values are dependent on the same base unit, the pixel.

In addition, these data need to be tested in several case studies on a global scale. That is why, one of the objectives of this study is to assess the potential of existing digital data for spatialized modeling in a GIS of soil water erosion in the current Tunisian context. It can be in some cases considered as an alternative to overcome the problem of lack of data.

The choice of RUSLE model (Equation 1) in this study was motivated by the fact that this model does not differ conceptually from the original soil loss equation (USLE) of Wischmeier and Smith [26]. On the contrary, it improves the quality of the environmental parameters which define the role of each factor. Also, the input data for this model is easier to access compared to other more recent models, which require more sophisticated data.

The RUSLE is written as: A=R*K*LS*C*P

(1)

Where:

A is the soil loss per unit area, expressed in the units of K and the period selected for R;

R is the rainfall and runoff erosivity factor;

K is the soil erodibility factor;

LS is the slope length and slope steepness factor;

C is the crop management factor;

P is the support practice factor.

RUSLE model was used to estimate the yearly average of soil loss due to water erosion, by taking into account the five factors indicated in equation (1). The latter has been applied using GIS and geospatial data composed of: (1) Climatic data including a 31-year average annual precipitation downloaded from the NASA-POWER meteorological parameters, taken from NASA's MERRA-2 assimilation model and necessary for computing the rainfall and Runoff Erosivity factor (R). (2) Soil map of Tunisia extracted from FAO Digital Soil Map of the Word (DSMW), essential for calculating the soil erodibility factor (K). (3) Global Digital Elevation Model (ASTER GDEM), which is used to calculate the topographic factor (LS). (4) LULC map of Tunisia extracted from ESRI 2020 Global Land Cover and used to calculate and map the vegetative cover factor (C) as well as the support and management practice factor (P).

Each of the factors was derived separately in raster format based on rainfall pattern, soil type, topography and LULC data in the context of soil water erosion modeling, in order to generate global scale maps. However, the current study, that deal with the estimation of soil erosion by water, did not consider development stage of the crops. Furthermore, the season variability effects were not taken into account.

The geospatial data with different types and origins are listed in Table 1. Figure 4 shows the flowchart of the method used for the estimation of each factor and for the analysis of soil loss based on the RUSLE model and GIS technique.

The following parts of this article make a brief description of the various factors showed in equation (1) and present the results obtained after modeling and analysis of geospatial data.

2.2.1. Rainfall and runoff erosivity factor (R-factor)

The estimation of the R factor requires knowledge of the kinetic energies and the average intensity over 30 min of the raindrops that fall each time, over a long period of up to 30 years [44]. But this direct method requires precipitation records at high resolutions to calculate R-factor. Other authors have developed alternative formulas when these data are not available. The equations most used to calculate the R factor using only annual precipitation are those of Renard and Freimund [45] (Equations 2 and 3), whose expressions are:

R = 0.0483. Pi^1.61               if 𝑃_𝑖 < 850 𝑚𝑚

(2)

R = 587,8 -1,219 Pi +0,004105. Pi ²        if 𝑃_𝑖 > 850 𝑚𝑚

(3)

Where: R is the measure of rainfall erosivity and Pi is the average annual precipitation (mm).

The estimation of rainfall and runoff erosivity using rainfall data with long-time intervals has been conducted by many authors for different regions of the world [36,46]. To calculate the average annual R-factor values in this study, a 31-year average annual data has been used and an interpolation of this data was applied to have a representative rainfall distribution map which is used as input for calculation of R-factor, knowing that the rainfall data available for the study area is not homogenous. This parameter, expressed in (MJ.mm.ha−1. h−1.year−1), takes into consideration the influence of climatic aggressiveness on soil loss [47].

To generate the R-factor for the whole Tunisia, a 31-year average annual precipitation from 1990 to 2020 required for this study, was downloaded freely, as a Comma-Separated Values (CSV) file, from the National Aeronautics and Space Administration (NASA) Prediction of Worldwide Energy Resources known as POWER project which was initiated in 2003. It provides access to solar and meteorological data sets for the entire globe and particularly over regions where surface measurements are sparse or nonexistent. The meteorological parameters are based upon the MERRA-2 assimilation model. The monthly and annual precipitation data are derived from NASA's Global Precipitation Measurement (GPM). and an interpolation of this data was applied to have a representative rainfall distribution map which was used as input for calculation of R-factor. Since the mean annual precipitation in Tunisia does not exceed 850 mm, equation 3 was used to calculate the rainfall-runoff erosivity factor, taking into account the bioclimatic map of Tunisia which has been developed for climatically homogeneous areas. The Inverse Distance Interpolation (IDW) method was used to minimize the calculation error and above all to optimize the processing time on the ArcGIS software. It is a method that consists in creating, from the aggregated points and the corresponding information, a value grid with a certain continuity.

2.2.2. Soil erodibility (K-factor)

The soil erodibility K, expressed in [t. h. MJ−1.mm−1], determines the resistance of different types of soil to erosion, knowing that soils are more or less sensitive to water erosion. This K factor is determined according to the characteristics of the soil: infiltration capacity, retention texture and susceptibility to particle removal. The infiltration and the high cohesion of the materials increase the soil resistance to removal and gullying of particles. The high rate of sand stabilizes the structure of the soil and makes it less sensitive to climatic aggression. Similarly, the organic matter improves the physical and chemical properties (cohesion, structural stability, porosity) of the soil, increasing the ability to retain water, and strengthens resistance to erosion [48]. Thus, the higher the percentage of sand, the more the soil is permeable, which implies a low value of the K-factor and vice versa. Bolline and Rousseau (1978) [49] established the classification, in Table 2 interpreting the soil susceptibility index.

Table 1. Description of the data sources.

Table 1. Description of the data sources.

Erodibility (K)	Type of soil
K < 0.10	Soil highly resistant to erosion
0.10 à 0.25	Soil fairly resistant to erosion
0.25 à 0.35	Soil moderately resistant to erosion
0.35 à 0.45	Soil with low erosion resistance
>0.45	Soil with very low resistance to erosion

However, soil erodibility is not constant, as it varies not only with soil type, but also with seasons and cultivation techniques. The soil map of Tunisia was prepared using FAO Digital Soil Map of the Word Shapefile (DSMW), where the soil data layer has been clipped, according to the country boundaries relative to the study area, in the ArcGIS Desktop environment. K-factor was calculated using Williams (1995) [50] formula (Equations 4) and FAO Digital Soil Map.

The raster dataset of the FAO/UNESCO Digital Soil Map of the World (DSMW) has a spatial resolution of 5 * 5 arc minutes and is in geographic projection.

$$K=f_{csand}\ast f_{cl-si}\ast f_{orgC}\ast f_{hisand}$$

(4)

Where:

${\mathrm{f}}_{\mathrm{c}\mathrm{s}\mathrm{a}\mathrm{n}\mathrm{d}}$ is a factor that gives low soil erodibility factor for soils with high coarse-sand content and high value for soil with little sand.

${\mathrm{f}}_{\mathbf{c}\mathbf{l}-\mathbf{s}\mathbf{i}}$is a factor that gives low soil erodibility factor for soils with high clay to silt ratios.

${\mathrm{f}}_{\mathrm{o}\mathrm{r}\mathrm{g}\mathrm{C}}$ is a factor that reduces soil erodibility for soils with extremely high sand content.

${\mathrm{f}}_{\mathrm{h}\mathrm{i}\mathrm{s}\mathrm{a}\mathrm{n}\mathrm{d}}$ is a factor that reduces soil erodibility factor with high coarse-sand contents and high value for soil with little sand.

The factors are calculated using these equations (Equations 5, 6, 7 and 8)

$$f_{csand}=0.2+0.3\mathrm{e}\mathrm{x}\mathrm{p}[-0.256{\mathrm{m}}_{\mathrm{s}}\left(1-\frac{{\mathrm{m}}_{\mathrm{s}\mathrm{i}\mathrm{l}\mathrm{t}}}{100}\right)]$$

(5)

$$f_{cl}-f_{si}={\left(\frac{{\mathrm{m}}_{\mathrm{s}\mathrm{i}\mathrm{l}\mathrm{t}}}{{\mathrm{m}}_{\mathrm{c}}-{\mathrm{m}}_{\mathrm{s}\mathrm{i}\mathrm{l}\mathrm{t}}}\right)}^{0.3}$$

(6)

$$f_{orgC}=1-\frac{0.25\mathrm{*}\mathrm{o}\mathrm{r}\mathrm{g}\mathrm{C}}{\mathrm{o}\mathrm{r}\mathrm{g}\mathrm{C}+\mathrm{e}\mathrm{x}\mathrm{p}[3.72-2.95\mathrm{o}\mathrm{r}\mathrm{g}\mathrm{C}]}$$

(7)

$$f_{hisand}=1-\frac{0.7\left(1-\frac{{\mathrm{m}}_{\mathrm{s}}}{100}\right)}{1-\frac{{\mathrm{m}}_{\mathrm{s}}}{100}+\mathrm{e}\mathrm{x}\mathrm{p}\left[-5.57+22.9\left(1-\frac{{\mathrm{m}}_{\mathrm{s}}}{100}\right)\right]}$$

(8)

Where, m_s is the percent sand content (0.05-2.0) mm diameter particles,

m_silt is the percent silt content (< 0.002) mm diameter particles,

m_c is the percent clay content (0.05-2.0) mm diameter particles,

orgC is the percent organic carbon of the layer (%).

2.2.3. Topographic factor (LS-factor)

The two parameters that constitute the topographic factor are slope gradient and slope length factor were estimated, here in this study, through an ASTER Global Digital Elevation Model (ASTER GDEM) with a 30 m resolution downloaded freely from (https://search.earthdata.nasa.gov/downloads/) [73]. The DEM of Tunisia’s territory is ranged from -27 m to 1549 m (Figure 2). The latter was incorporated into GIS to determine accurately the slope gradient (S) and slope length (L), taking into account that the effect of slope length and degree of slope interaction should always be considered together [51]. So, to create the slope length and steepness (LS-factor), equivalent to topographic factor and relief factor, respectively, flow direction and flow accumulation maps were created after the filling process of DEM data by the use of the GIS extension Arc Hydro Tools. The flow accumulation was determined for each cell by the flow that flows through that cell. The greater the flow accumulation value, the easier the area will form runoff and vice versa. The slope map in degree, required to estimate the LS factor, shows that most of the slope of Tunisia is ranged from 0 to 10° which represent about 95,54 % from total area (Figure 9).

In 1985, Moore and Burch [52] developed the equation below (Equation 9) to compute the length-slope factor:

$$LS={\left(\frac{\mathrm{A}\mathrm{s}}{22.13}\right)}^{\mathrm{m}}\times {\left(\frac{\mathrm{s}\mathrm{i}\mathrm{n}\mathsf{\theta }}{0.0896}\right)}^{\mathrm{n}}$$

(9)

The equation (9) has been applied by many researchers [51,53,54]. The exponent value m can be taken equal to 0.4 while the value of n can be taken equal to 1.3 [53,54,55].

Where:

𝐿𝑆: is slope steepness–length factor.

𝐴𝑠: is the flow accumulation (in meters).

𝜃: is the slope angle (in radians).

𝑚 = 0.4 – 0.6 and n = 1.2 – 1.3.

In general, as the length of the slope increases, the total soil erosion and the soil erosion per unit area increase due to the gradual accumulation of runoff water in the downhill direction of the slope. As the steepness of the slope increases, the velocity and erosivity of the runoff increase.

2.2.4. Vegetative cover factor (C-factor)

It is a dimensionless factor presenting the effectiveness of the vegetation cover in relation to the susceptibility of the soil to erosion [56]. Vegetation cover and its spatial distribution play an important role in reducing the effects of runoff by amortizing the impact of rainwater on an area [15].

According to several studies, this RUSLE factor that can range from zero for water and a very well protected soil with very strong cover effects, to 1 for a surface that produces a lot of runoffs (bare soil) and leaves the soil very susceptible to water erosion [57]. Several authors consider that C-factor is around 0.01 (1/100) under dense forest, 0.05 (5/100) under grasslands and 0.24 (24/100) under crops (Table 3). Other researchers have adopted the calculation of C-factor by new simplified approaches: use of remote sensing techniques such as the classification of satellite images [58,59] and vegetation indices [60] or use of the Land Use/Land Cover (LULC) map and class assignment for each entity [61]. In fact, C-factor reflect the effect of LULC, cropping and management practices on the rate of soil erosion [62,63,64,65].

C-factor was determined, in this study based on the literature, by assigning a value for each type of LULC extracted from ESRI 2020 Global map of LULC, derived from ESA Sentinel-2 imagery at 10 m resolution and downloaded freely from https://livingatlas.arcgis.com/landcover/ for all land masses on the planet. Two Individual GeoTIFF scenes (32S_20200101-20210101 and 32R_20200101-20210101) covering the Tunisian territory, have been downloaded from Esri 2020 Land Cover Downloader application. These scenes were joined to form a mosaic image and the values of the C-factor (Table 3) for 7-class LULC types covering the study area were used in the data analysis to build the C-factor map for the year 2020 at 10 m resolution, performing the analysis in ArcGIS and then converted to a grid with 30 x 30 m spatial resolution.

2.2.5. Support and management practice factor (P-factor)

This factor, which is dimensionless, represents soil protection based on anti-erosion cultivation techniques that reduce runoff speed and thus reduce the risk of water erosion. It varies according to the landscaping carried out, namely cultivation on a level curve, in alternating strips or on terraces, reforestation in benches, ridging and ridging [15].

The P-factor is between 0 and 1, in which the value 0 represents a very good environment of resistance to erosion and the value 1 shows an absence of anti-erosion practice [72]. Wischmeier and Smith (1978) [15] established the classification according to the slope in percent (Table 4), showing that this factor can be distributed according to two zones: the cultivated zones and the other zones.

3. Results and discussion

In this part, the results of the processing as well as the calculations of each factor made for the spatialization of the soil loss rates in the study area will be presented in more details. After applying the RUSLE model, an assessment of the impacts of erosion at the regional scale of Tunisia will also be detailed in this section.

3.1. Spatial distributions of RUSLE factors

3.1.1. R-factor

The rainfall erosivity factor quantifies the effects of rainfall aggressiveness and therefore reflects the amount and rate of runoff associated with a rainfall event. It was calculated from annual average precipitation data recorded over a time interval of 31 years (from 1990 to 2020).

After calculation, the values of the R factor in Tunisia range from 24.29 to 1859.84 MJ.mm.h^-1. ha^-1. year^-1. However, the average value in the entire country is 473.73 MJ.mm.h^-1. ha^-1. year^-1. The maximum value is observed in the North-West part of Tunisia and the minimum value is observed in the South.

The maps in Figure 5 and Figure 6 show that precipitation and erosivity factor (R-factor) decrease gradually from the South towards the extreme North-West part of Tunisia. In addition, it is important to note that rainfall and runoff erosivity factor is a crucial factor in assessing soil erosion for future LULC and climate change.

3.1.2. K-factor

The physic-chemical properties of soils, such as organic matter, fine sand, clay and silt contents were obtained from the FAO-UNESCO Soil Map of the World. The percentage of organic matter (OM) depends on the presence and nature of the vegetation on the ground as well as on the use of the latter.

Determination of the K factor was performed using equations 4, 5, 6, 7 and 8 from Williams (1995) [50] and the soil map of Tunisia (Figure 7), extracted from the FAO-UNESCO Digital Soil Map of the Word (DSMW) shapefile [74]. The soil erodibility spatial distribution map (Figure 8) shows that the values of K-factor range from 0.10 to 0.19 t. h. MJ^-1.mm^-1. According to the classification of Bolline and Rosseau (1978) [49], the soils in Tunisia are therefore all classified as soil fairly resistant to erosion.

3.1.3. LS-factor

The slope map (Figure 9) is obtained by the Slope tool in Spatial analyst which was then used to calculate the steepness of the slopes. The mean value of the slope throughout the country is 2.24° with a standard deviation of 3.50° and values above 10° especially in the South-East, Center and North-West of the country.

Knowing that the LS factor is the derivative product of equation (9) of Moore and Burch (1985) [52], the equation of this product is formulated on Map Algebra's Raster Calculator following the syntax:

LS = Power ("Flow Accumulation" * Cell size / 22.13,0.4) * Power (Sin ("Slope in degree" * 0.0174533) / 0.0896,1.3)

Such that “Slope in degree” represents the slope map as represented by Figure 9 and the cell size represents the spatial resolution of the ASTER GDEM equal to 30 m. To convert the slope to radians, it was multiplied by 0.0174533, since 1 degree = 0.0174533 radians.

Thus, the following LS topographic factor map was obtained (Figure 10), showing that the value range of the LS factor is between 0 and 12.41. However, the average value in the whole country is 1.09.

3.1.4. C-factor

The vegetation cover factor map shows that the values vary from 0 to 1 over the whole country. The minimum values are found in the northern part of the country, where according to the LULC map (Figure 11), there are mainly dense forests and cultivated areas corresponding to well protected soil with very strong cover effects. On the other hand, the value of C-factor (Figure 12) is close to 1 going south, especially on bare soil that produces a lot of runoffs and leaves the soil very susceptible to water erosion.

3.1.5. P-factor

The support and management practice factor explains the effects of agricultural practices that can minimize the impact of rainwater and reduce the rate of runoff, and thus logically reduce soil loss (Wischmeier and Smith, 1978) [15].

The determination of the P-factor required the establishment of the map of agricultural and non-agricultural land (Figure 13) from the LULC map of Tunisia extracted from ESRI 2020 Global map, derived from ESA Sentinel-2 as well as the slope map in percent.

The values of the P-factor in agricultural areas decrease as the slope decreases and vice versa on steep terrain. These values vary from 0.1 on low slope land to 0.33 on very steep land. The other non-agricultural areas are all classified as areas without erosion control practices, which are assigned the value P = 1 (Figure 14).

The results of the calculations on each factor can be summarized as follows:

• The values of the erosivity factor of the rains vary from 24.29 to 1859.84 MJ.mm.h^-1. ha^-1. year^-1 with an average of 473.73 MJ.mm.h^-1. ha^-1. year^-1 over the whole country.
• The soil erodibility K-factor is classified as fairly resistant to erosion, with values varying from 0.10 to 0.19 t. h. MJ^-1.mm^-1.
• The value range of the topographic factor LS is between 0 and 12.41, with an average value of 1.09 over the whole country.
• The values of the vegetative cover C-factor vary from 0 to 1.
• The values of the support and management practice factor vary from 0.1 to 0.33 for agricultural land and are equal to 1 for non-agricultural land.

3.2. Calculation of soil loss rate and quantification of erosion

The modeling of the factors involved in water erosion of soils has enabled the quantification and spatialization of this phenomenon.

The RUSLE model gives an approximation of the average value of soil loss, expressed in tons per hectare per year, on a given scale of study. This approximation is based on the joint product of the five factors described above. In this study, the modeling is essentially carried out on ArcGIS. To estimate and map the spatial distribution of soil losses for the entire Tunisian territory, the maps of the five factors, seen in the previous sections, were projected according to the same coordinate system "Universal Transverse Mercator (UTM) zone 32 N using the Datum: WGS 1984" with a spatial resolution of 30 m x 30 m for each. We then proceeded to the calculation pixel by pixel by the RUSLE equation (1) with the "Raster calculator" tool and the "Map Algebra" function in ArcMap.

The application of this model allowed the establishment of the map of Figure 15, displaying the average values and the distribution of soil loss rates for the whole country.

By carrying out the calculation according to the empirical and spatialized model RUSLE on ArcGIS, the average value range of soil losses in Tunisia was obtained, varying from 0 to a maximum value of 619.44 t/ha/year. These results have been subdivided into 5 classes, according to Kefi et al. (2012) [9], represented in Table 5.

Despite a large domination of very low soil erosion by water, representing 56.34 % of the total area of Tunisian territory, 10.64 % of the total area are severely exposed to water erosion of the soil (High and very high soil erosion).

4. Conclusions

This study was the first to provide complete and detailed water soil erosion map for Tunisia at a country scale. The RUSLE model, an empirical model to assess soil losses per year, has been conducted to evaluate soil losses in Tunisia which were mapped at a 30 m cell size. The usage of the RUSLE equation is significantly facilitated by its deployment in Geographic Information Systems (GIS) and the use of geospatial data, provided over vast regions, such as the global territory of Tunisia. So, the use of this approach contributed to regional soil erosion risk assessment through the rapid derivation of appropriate indices.

Therefore, this study has found an adequate method that integrates geospatial data and GIS with the RUSLE model to estimate the spatial distribution of soil water erosion in whole Tunisia and provide spatial erosion map. The quantification of water soil erosion losses is an essential approach to spatialize the zones most sensitive to water erosion. Consequently, this approach can give, in practice, relevant results for the potential evaluation of soil losses at the regional scale and make an important contribution to regional soil erosion risk assessment through the derivation of appropriate factors, despite the few studies that have been carried out for the modeling of soil erosion at the scale of countries around the world.

Soil erosion by water is a natural phenomenon that may or may not be accelerated by various agents. The impact of soil erosion depends on several parameters. In this study, the rates of soil loss by water erosion were listed according to five different classes, including: very low, low, moderate, high and very high.

There is normally a tolerable amount of soil loss in the natural cycle of the earth's elements, so that soils can maintain their long-term productivity. Tolerable soil loss is the maximum annual amount of soil that can be removed without affecting the long-term natural productivity of soil covering the slope of an area. The limit value is set in this study at 10 t/ha/year.

The results indicated that 6.43% of the Tunisia’s surface, corresponding to 9989.66 km2, was affected by a very high soil loss rates (>30 t ha−1 y−1). 4.20% of the total area of the country, corresponding to 6535.02 km2, was affected by a high soil loss rates ranging from 20 to 30 (t ha−1 y−1). These results can be useful to identify the vulnerable areas of the country and to develop a less time-consuming decision plan more efficient for soil erosion prevention and control. This integrated approach, based on GIS and erosion model, can be useful to contribute to a better understanding of the soil degradation of large-scale areas in Tunisia, for better preserving and protecting both natural and man-made resources available in the country.