Identification of zones at risk of erosion by the combination of a digital model and the method of multi-criteria analysis in the arid regions: case of the Bechar Wadi watershed

Abstract

Erosion and solid transport is a tricky and complex problem that negatively affects natural and urban environments. In Algeria, the effects of this phenomenon are apparent; their impact is no less devastating in the long term than the other spectacular catastrophic phenomena that can be observed. Sixty-five large dams in Algeria are threatened by the reduction of 62% of their storage capacity because of the siltation problem (ANBT) (National Agency for Dams and Water Transfers). The main objective of this work is the evaluation of the impact of the erosion phenomenon on Bechar watershed which is in an area characterized by an arid climate. The universal soil loss equation was used. This model is based on the combination of the five factors (erosivity, erodibility, topography, vegetation cover and support practices) that directly influence this phenomenon. Analytical hierarchy process is used to give a weighting value of each factor according to its degree of influence on the phenomenon. The sediment delivery ratio is calculated to determine the amount of soil that will arrive at the outlet of the watershed and contribute to the storage structures siltation. The obtained results will undoubtedly help decision makers to understand the threat of erosion degree in the study area in order to better take the necessary measures to face this issue.

Introduction

Soil erosion is a complex natural phenomenon that threatens soil stability. Water and wind are estimated as the main agents that arouse the appearance of this phenomenon, especially in arid zones (Balasubramanian 2017). It has become a visual problem, subject to a combination of factors such as rainfall intensity, soil type, vegetation cover and other parameters that regulate the intensity of soil erosion (Carvalho et al. 2015). The consequences of this phenomenon are numerous; among them, main consequences are the loss of soil, the reduction of the quality/quantity of the water and the flooding risk increase (Panagos et al. 2015), which classify the phenomenon of soil erosion among the world’s most important environmental problems (Pimentel 2006). In addition, researchers, in the environmental field in Algeria, consider erosion as the first factor affecting the water reserves of storage facilities in the country (Koussa and Bouziane 2018; Touahir et al. 2018).

In order to properly determine the impact of the erosion phenomenon and allow decision makers to intervene with appropriate solutions to reduce these negative effects on agriculture, infrastructure, water quality, etc. (Abdi et al. 2013), the computer tool is essential in order to draw up soil erosion risk assessment maps to distinguish the most vulnerable zones (Deepanshu et al. 2016) by using the geographic information systems (GIS) and digital terrain models (Rahman et al. 2009).

The Bechar watershed is vulnerable to soil erosion because of its location in an arid zone characterized by rainfall irregularity and dry climatic behavior (Boufeldja 2013). The soil erosion is a natural phenomenon, which has consequences for the agricultural sector and on the area water resources. Work that addresses this topic of estimating the soil erosion rate in the region is very limited. Belkendil (2014) and Bouzouina et al. (2014) presented contribution studies, estimating the erosion phenomenon in the Guir and Zousfana watersheds, respectively. The Bechar watershed is located between these two basins. The aim of this work is to determine the influence of soil erosion phenomenon on Bechar watershed by quantifying the long-term average annual land losses in tons per acre per year (A) and sediment yield (SY) using the USLE–SDR coupled model, with the multi-criterion decision analysis (MCDA) methods: analytical hierarchy process (AHP) (Thomas et al. 2018 and Algarín et al. 2017) to calculate the weighting of various factors of the erosion phenomenon according to the grade of importance of each factor. AHP is a powerful and systematic method for decision making, where a hierarchical structure is maintained among the objective, decision-making criteria and alternatives by successive levels (Saaty 1980, 1990).

Presentation of the study area

The Bechar watershed covers a 6357-km² area and is drained by Oued Bechar (Fig. 1). Located at the foot of the southern slope of the Saharan Atlas, it is limited to the north by the mountain chain of Jebel Grouz, to the southwest by the Ougarta Mountains, to the southeast by the Grand Erg Occidental, to the west by the Guir Hamada and to the east by the Oued Zousfana Hamada. It originates in the Jebel Grouz, at an altitude of 1590 m; it travels from northeast to southwest about 220 km.

Fig. 1

The location of the watershed

The watershed is in the arid to semiarid climate zone dominated by drought. The average annual temperature varies from 7 to 50 °C. Average annual rainfall ranges from less than 40 mm to over 100 mm in the northeastern part of the region (Boufeldja 2013). The topography of the region is generally flat except the eastern part which has a mountainous and rugged character (Djebel Bechar).

There are three major soil classes in this area. These are silts, clays and silty mud covering approximately 79%, 16% and 5%, respectively (HWSD 2012) (the land use pattern of the region includes 3.03% of forest land, 4.48% of urban land and 92.49% of bare land) (Fig. 3c).

Methodology

Annual estimate of soil loss by model (USLE)

The methodology is based on the universal soil loss equation (Eq. 1) that was established by Wischmeier and Smith (1978) to be applicable worldwide (Laflen et al. 2003). USLE is used with a computer program (Remortel et al. 2001). This model provides the best results for predicting erosion rates in ungauged watersheds, using knowledge of watershed characteristics and local hydro-climatic conditions (Deepanshu et al. 2016).

$$A = R \times K \times {\text{LS}} \times C \times P$$

(1)

where A is the long-term average annual land losses in tons per acre per year. This value can then be compared to the “tolerable soil loss” limits, R is the erosivity factor, K is the soil erodibility factor, LS is the topography factor, C is the plant cover factor and P is the factor of anti-erosion practices.

Erosivity factor R

The R factor represents the effect of rain and runoff on the soil, and it is the result of a product of the precipitation energy (E) and the maximum intensity for 30 min (I₃₀) (Nearing et al. 2017). Following his work in Morocco, Arnoldus partitioned Africa and the Middle East into climatic zones based on the ratio of annual rainfall to potential evapotranspiration. He then used the modified Fournier index to create an iso-erodent map in metric units for Africa north of the equator and the Middle East (FAO/UNEP/UNESCO 1979; Arnoldus 1980). In humid arias where there were not enough stations with calculated values, the relationships used in the less humid parts were extrapolated. Similarly, a relationship between the Sebou basin in Morocco has been extrapolated to the driest regions of Africa (Arnoldus 1980). Fournier uses monthly and annual average precipitation data according to the following regression:

$$\begin{aligned} & R = 0.264 * F^{1.50} \\ & F = \mathop \sum \limits_{i = 1}^{12} \frac{{p_{i}^{2} }}{{\bar{p}}} \\ \end{aligned}$$

(2)

where $p_{i}$ is the monthly rainfall and $\bar{p}$ is the annual precipitation.

In this study, we used 45 years (1955–2000) of data from nine weather stations located around/in the study area provided by Worldclim (FAO¹). According to the results obtained (Table 1), the values of R are between 4.43 and 27.5 MJ mm ha⁻¹ h⁻¹ year⁻¹ (Fig. 2a), the maximum value was recorded in the Bouarfa station and the minimum value in the station of Igli, and the average erosivity value in the studied watershed is 9.32 MJ mm ha⁻¹ h⁻¹ year⁻¹.²

Table 1 Erosivity values and location of meteorological stations used

Fig. 2

USLE parameters a R factor map. b K factor map. c LS factor map. d C factor map. e P factor map

Soil erodibility factor (K)

The K factor represents the degree of soil sensitivity to erosion. It measures the cohesion of the soil on an exploited field. Wischmeier and Smith (1978) used a 22.1-m-long pilot field with a 9% slope continuously maintained in summer fallow (Pham et al. 2018). The K value depends on the soil physical and chemical properties, such as texture, shear strength, permeability, grain size and organic matter content (Williams 1995), giving an equation for estimating k_USLE values

$$K_{\text{USLE}} = K_{W} = F_{\text{csand}} \times F_{{{\text{cl}} - {\text{si}}}} \times F_{\text{orgc}} \times F_{\text{hisand}}$$

where F_csand is a factor that lowers the K indicator in soils with high coarse-sand content and higher for soils with little sand; F_cl–si gives low soil erodibility factors for soils with high clay-to-silt ratios; F_orgc reduces K values in soils with high organic carbon content; and F_hisand lowers K values for soils with extremely high sand content:

$$F_{\text{csand}} = \left( {0.2 + 0.3 \times \exp \left[ { - 0.256 \times m_{\text{s}} \times \left( {1 - \frac{{m_{\text{silt}} }}{100}} \right)} \right]} \right)$$

$$F_{{{\text{cl}} - {\text{si}}}} = \left( {\frac{{m_{\text{silt}} }}{{m_{\text{c}} + m_{\text{silt}} }}} \right)^{0.3}$$

$$F_{\text{orgc}} = \left( {1 - \frac{{0.25{\text{oegC}}}}{{{\text{orgC}} + \exp \left[ {3.72 - 2.95 \times {\text{orgC}}} \right]}}} \right)$$

$$F_{\text{hisand}} = \left( {1 - \frac{{0.7 \times \left( {1 - \frac{{m_{\text{s}} }}{100}} \right)}}{{\left( {1 - \frac{{m_{\text{s}} }}{100}} \right) + exp\left[ { - 5.51 + 22.9\left( {1 - \frac{{m_{\text{s}} }}{100}} \right)} \right]}}} \right)$$

where m_s, is the sand fraction content (0.05–2.00 mm diameter (%), m_silt is the silt fraction content (0.005–0.05 mm diameter (%), m_e is the clay fraction content (< 0.002 mm diameter) (%), orgC_c is the organic carbon (SOC) content (%)

The assumptions used to calculate K_USLE value are presented in Table 2.

Table 2 Method to calculate K_USLE value

The different soil fractions of the study area were obtained from the DSMW³ database. According to this database, the Bechar watershed is divided into three soil classes (clay, silt and silt), and the erodibility values vary between 0.020875 and 0.02276 t ha h ha⁻¹ MJ⁻¹ mm⁻¹ (Fig. 2b). The higher value was observed in the northern part of the Bechar watershed which has a rocky mountainous character, and the maximum value was observed in the southern part near the watershed outlet where there are strong sediments deposits which are easily erodible.

LS topography factor

The LS factor shows the influence of the length and the slope inclination on the erosion phenomenon. Fauck (1956) and Fournier (1967) claimed that a very low slope, in the order of 2%, can trigger the water erosion phenomenon. Wischmeier and Smith (1978) demonstrate that the distance of the slope is equal to the trigger distance of the flow until the beginning of the settling phase. Figure 5a shows that slope values in Bechar watershed vary between 0° and 58°, while the value of the flow accumulation ranges from 0 to 52 639. Moore and Wilson (1992) presented a modified equation, based on the original Wischmeier and Smith equation (1978); the objective of this modification is to be able to use this equation on the specialized software ArcMap (10.3) by using the numerical modeling of terrain of 90 m resolution. This value represents the same value of the cell size in the following equation:

$${\text{LS}} = \left[ {{\text{flow}}\;{\text{accumulation}}*{\text{cell}}\;{\text{size/}}22.13} \right] 0.4* \left[ {\sin \left( {{\text{slope}}*3.14 /180} \right) /0.0896} \right]^{1.3}$$

(5)

The result of applying this equation, using the raster calculator option, with the ArcMap software, has shown that LS values vary between 0 and 122.60 (Fig. 2c). The LS values’ distribution map shows that there is a dominance of LS low values, especially in the central zone and the southwestern part of the watershed; in return, the maximum values of LS were recorded in the northern part of the watershed (Djebel Antar, Djebel Horreit), and the average value of LS is 0.23.

Plant cover factor C

The plant cover factor (C) is estimated as the most important factor in the erosion phenomenon (Weiwei et al. 2011). It represents the positive effects of the vegetation cover on the soil particles stability and thus the soil losses reduction, by their actions characterized in the kinetic energy absorption of raindrops and the decrease in runoff. The effects of factor C vary with time and nature (Wenwu et al. 2013). We try to determine the value C using high-resolution Landsat-7 ETM⁺ satellite images. The image processing operation was done using ArcMap software using the supervised classification options for satellite images. The classification result shows that there are three classes (bare land, urban area and lean vegetation) (Fig. 2d). The values of the factor for these three classes vary between 0.01 and 0.35, and the average value of the factor C is 0.32.

Factor of anti-erosion practices P

The factor P represents the human intervention using the necessary facilities for the purpose of reducing the erosion rate, by the adjustment of the flow, the slope, the direction of runoff and the lowering, and therefore reducing the amount of sediment transported (Wischmeier and Smith 1978; Renard and Foster 1983). P is the ratio of land loss associated with conservation practice to land loss associated with line farming in the direction of slope. The most commonly used conservation practices are tillage against slopes, contour cultivation, The values of P vary from 0 to 1, the value 0 showing that the anti-erosion arrangement is perfect, however the value of 1 (Ashaq et al. 2011). To calculate the P-factor value in the Bechar watershed, we used the P-values based on the crop types and slopes given by Shin and Pesaran (1999) (Table 3), the topographic map and the watershed soil use (Fig. 2e).

Table 3 Supporting practice factor values by crop type and slope (Shin and Pesaran 1999)

The P values (Fig. 2e) vary between 0.55 and 1. The average value is 0.59 with a standard deviation value close to 0.12.

Sediment delivery ratio (SDR)

The sedimentation delivery ratio (SDR) is a fraction of the total erosion transported from a given area in a given time interval. This is the amount of sediment actually transported from the erosion zones to the watershed outlet. The SDR value in a given watershed indicates the basin’s ability to store and transport eroded soil through the compensation of sediment deposits that become increasingly important with the extension of the watershed surface and therefore determines the relative importance of sediment sources and their contribution (Sewnet 2016).

The SDR estimation methods are diverse. They have been developed on the basis of the watershed variable physical characteristics (Wu et al. 2017). In this paper, three SDR relationships are used: that of Maner (1958), that of Vanoni (1975) and that of Williams and Berndt (1972) (Table 4). Thus, we used the EPM⁴ and the ratio (R_u) to choose the appropriate model based on the comparative methods, the standard error (SE) and the coefficient of variation (CV) (Rostami et al. 2001).

Table 4 SDR results of the models used

The difference in the results obtained for the SDR (Table 4) lies in the variability of the factors used in each method. For comparison, we calculated the ratio of the sediment input (R_u) of the EPM model given by:

$$R_{u} = 4\frac{{(P \times D)^{0.5} }}{L + 10}$$

(6)

where $R_{u}$ is the sedimentation coefficient in the watershed; L is the length of the straight line joining the two ends of the watershed; P is the perimeter of the catchment area in km; D is the difference between the mean and minimum altitudes of the watershed, which is given as follows

$$D = D_{\text{ar}} - D_{0}$$

(7)

where D₀ is the altitude at the outlet in km and D_ar is the average altitude of the catchment area in km.

Thus, the sediment input ratio R_u = 0.369 (Table 5).

Table 5 The watershed estimated sedimentation rate by the EPM model

For the analysis and selection of the appropriate model for the study area, different statistical tests, including adaptive comparisons, standard error (SE) and coefficient of variation (CV) (Table 6) with respect to type, nature and relevant data, are used.

Table 6 Comparison of the results of the models evaluated against the basic SRR (0.27) and their ranks

$${\text{CV}} = \frac{\text{SD}}{{X_{0} }} \times 100$$

(8)

$${\text{SE}} = \sqrt {\frac{{\left( {X_{\text{e}} - X_{0} } \right)^{2} }}{{X_{0} }}}$$

(9)

where CV is the coefficient of variation; SE is the standard error; X₀ is the observed SDR (SDR₀); X_e is the estimated SDR (SDR_e).

Assessment of Soil Sensitivity to Erosion Index (SSEI) and analytical hierarchy process (AHP) application

The soil sensitivity to erosion index (SSEI) characterizes the influence of a multitude of environmental parameters on the phenomenon of soil erosion. This index evaluation interest is to determine the degree of the risk of this natural phenomenon as well as the detection of the erosion vulnerable zones.

Several authors have calculated this index. The choice and the selection of environmental parameters are done with the help of experts combined with an analysis of the studied area to know the parameters that can most influence the erosion phenomenon. Khatun (2017) assessed the SSEI in the Kushkarani basin using ten main parameters (drainage frequency, drainage density, slope, land use/soil cover, soil texture, hydraulic gradient, elevation, precipitation, NDVI⁵ and geology). Pradeep et al. (2014) used seven geo-environmental variables such as slope, relative relief, land use/land cover, landform, drainage density, drainage frequency and lineage frequency.

In this study, nine environmental parameters (such as annual mean precipitation, soil texture, land use/land cover, landform, terrain geomorphology, drainage density, drainage frequency, lineament frequency, slope and related terrain) were integrated into the analytical hierarchy process (AHP) using the geographic information systems (GISs) platform to generate the soil erosion map of the region vulnerable areas.

The annual precipitation parameter was used to show the importance of its role in the erosion phenomenon as the first driver of this phenomenon through the impact of raindrops on the soil surface during high-intensity storms increasing the soil particles loosening (Mohamadi and Kavian 2015). The Bechar watershed precipitation distribution map (Fig. 3a) is based on Worldclim data (1955–2000) from nine meteorological stations (Table 1). Inverse distance weighted interpolation (IDW) was used through the ArcMap toolbox software. According to the generated map, the annual rainfall values in the studied watershed range from 67.06 to 149.67 mm. The interval between these two values has been divided into four classes in the AHP comparison matrix where the class between 120 and 150 mm has the largest weight value, which is 48% (Table 10).

Fig. 3

soil Sensitivity to Erosion Index parameters. a Weighted average annual rainfall. map. b Soil texture map. c Soil use map/vegetal cover

Soil texture is the parameter that refers to the relative proportions of the soil (clay, silt and sand). Indeed, the texture influences erosion by the increase or decrease ratio of the soil components (Easton and Emily 2016). Soils with low clay content are less cohesive and inherently more unstable. These soils are more exposed to water and wind erosion. The FAO/DSMW database was used to develop the soil texture distribution map in the Bechar watershed (Fig. 3b). This map has three types of texture (coarse, medium and fine) with more than 79% of the area being of medium texture. This dominant class has a 24% weight in the AHP comparison matrix.

The change in vegetation cover (LUC) has important impacts on soil degradation, especially erosion (Sharma et al. 2011). Several components form the vegetation cover (forests, agricultural land, urban area, bare land, etc.), and each component has a specific value, basing its impact degree on the erosion phenomenon. To determine the vegetation cover distribution in the study area, a supervised classification operation was performed using Landsat TM three-band satellite imagery (red, green and blue) using the ArcMap classification tool.

The classification process revealed that the study area is divided into three classes (Fig. 3c): forests, urban area and bare land. This latter class dominates as the study area is subject to the arid climate regime. The use of satellite images has also contributed to the detection of geological structures (lineaments) that strongly influence terrain stability through increased permeability of the terrain and, consequently, soil moisture (Ali et al. 2018), which worsens the soil loss phenomenon in the region.

To determine the lineaments zones, in the Bechar watershed, a Landsat TM satellite image processing process, using the Geomatica software, reveals the existence of a high density of lineaments both in the northern part of the watershed (Djebel Antar) than in the eastern part (Djebel Bechar) and near the town of Kenadsa (Fig. 4a). The maximum value of the lineaments is 11 No/km².

Fig. 4

soil Sensitivity to Erosion Index parameters: a lineament frequency map. b Field forms map. c Relative relief map

The geomorphological position makes it possible to determine the infiltration (Cerdà 1998), which has a direct impact on the runoff and subsequently on the erosion phenomenon. Chabala et al. (2013) described the mapping methodology of the geomorphological structure using the digital terrain model. This methodology application revealed the existence of five geomorphological classes in the study area (Fig. 4b), and the class of summits (slope of the hills) takes the maximum value of the weights in the matrix of the AHP (Table 10). The area relief and slope are the result of various geomorphological processes occurring within and on the Earth’s surface (Sharma et al. 2018).

The topographic factor plays a key role in the soil erosion modeling operation; the relative relief parameters and the slope were calculated using the 90 m elevation digital elevation model (DEM). In the study area, the slope varies between 0° and 58° (Fig. 5a), while the relative relief varies between 3.08 and 557.12 m/km² (Fig. 4c). The respectively average slope and the relative relief values are 2.85° and 46.85 m/km², respectively.

Fig. 5

soil Sensitivity to Erosion Index parameters a slopes map. b Drainage frequency map. c Drainage density map

Drainage density is estimated as the total length of streams per unit area of the watershed and depends on factors such as lithology, permeability and vegetation (Moeini et al. 2015). Density of drainage has a direct effect on land stability, especially in mountainous areas due to the large contribution of groundwater recharge, which causes landslides (Pradeep et al. 2014). A drainage density map has been established for the study area (Fig. 5c) where drainage density values vary between 0 and 9.27 km/km². These values were divided into four classes to form the AHP comparison matrix (Table 10).

The frequency of drainage is the parameter that indicates the number of streams per unit area (km²). According to investigative work, large drainage frequency values are found in non-porous soil type zones. These areas are characterized by high slopes, intense rainfall, low vegetation cover and, subsequently, high degree of erosion (Kumar 2017).

According to the drainage frequency map of the Bechar watershed (Fig. 5b), we can see that the drainage frequency varies between 0 and 11 No/km²; the maximum class (> 9 No/km²) in the AHP comparison matrix corresponds to a relative weight of 54%.

Multi-criteria analysis methods help in decision making. The decision maker relies on criteria or factors that influence more or less, in a direct or indirect way, on the issue; here, the decision-maker mind reading is developed since 1960, including the methods of ELECTRE, PROMETHEE, AHP (Caillet 2003) and MACBETH (Bana e Costa and Beinat 2005; Hadji 2013). In this study, we chose the hierarchical analysis of AHP (analytical hierarchy process), a simple and easy method known for its contribution in several domains (Bhushan and Rai 2004), created and developed in 1980 by Saaty (Bernasconi et al. 2010).

The AHP relies on comparisons of the essential elements in a decision so that their ranking is a priority. It can be used in several areas (Sabaei et al. 2015). The advantages of this method are numerous; among them, the main advantages are the simplicity of application, the easy accessibility of the inputs and the comparison of several parameters at the same time (Saaty 1995).

The first step is to break down the complex problem into sub-criteria that react on the problem in a successive hierarchical ranking (Fig. 6).

Fig. 6

Soil sensitivity assessment (SSEI) methodology for erosion using AHP methods

The second step is to present the decision-maker role, who makes a pairwise comparison of these criteria based on his experience or information collected from expert work. The decisions evaluation will be linguistic and will be transformed into numerical values. The principle is based on a scale of absolute values ranging from 1 to 9 and their inverse, as given in Table 7.

Table 7 Recommended scale of comparison

The third step is to build the comparison matrix. Decisions to an intangible share differ from one expert to another depending on the action history or the experience. The values of the comparison will be organized in a matrix.

$$\begin{aligned} & M = \left[ {\begin{array}{*{20}c} 1 & {a_{12} \cdots } & {a_{1j} } \\ {\begin{array}{*{20}c} {a_{21} } \\ \vdots \\ \end{array} } & 1 & {a_{2j} } \\ {a_{i1} } & {a_{i2} \cdots } & 1 \\ \end{array} } \right] \\ & {\text{Number}}\;{\text{of}}\;{\text{comparisons }} = \frac{{n\left( {n - 1} \right)}}{2} \\ \end{aligned}$$

(10)

The value taken for one element corresponds to the degree of importance or the force of stress relative to the other (Saaty 2008). The determination of weights is done using the following formula (Pradeep et al. 2014).

$$a_{ij} = \frac{{w_{i} }}{{w_{j} }} = \frac{{{\text{attribute}}\;{\text{weight}}\; i }}{{{\text{attribute}}\;{\text{weight}}\;j}}$$

(11)

A consistency assessment is necessary after selecting the study parameters weights. The consistency ratio (CR) is the weights homogeneity measurement parameter assigned for the classes of each factor. For the matrix to be valid, the consistency ratio (CR) must be less than 0.1 (Saaty 1995) (Table 8). CR can be evaluated by the following formula:

$${\text{CR}} = \frac{{\text{CI}}}{{\text{RI}}}$$

(12)

where CI is the consistency index of the normalized matrix, calculated using the following formula:

$${\text{CI}} = \frac{{\lambda_{{\rm max} } - 1 }}{n - 1}$$

(13)

with λ_max being the largest matrix eigenvalue obtained from the priority matrix and n is the size of the comparison matrix.

Table 8 CI and CR values

RI is the average coherence index of the reproduced random comparisons (Saaty 1980) or simulations (Caillet 2003) (Table 9).

Table 9 Mean consistency index values for comparisons according to matrix n size

$\lambda_{\hbox{max} } = \left( {2.743 \times 0.32} \right) + \left( {5.767 \times 0.22} \right) + \left( {8.833 \times 0.14} \right) + \left( {19 \times 0.006} \right) + \left( {23 \times 0.04} \right) + \left( {22.5 \times 0.05} \right) + \left( {31 \times 0.03} \right) + \left( {13.58 \times 0.10} \right) + \left( {23 \times 0.04} \right) = 9,77$ The coherence index resulting value shows a high consistency ratio between the selected decisions, since it falls below the threshold value (0,1), so weighting values can be accepted.

The geo-environmental parameters used were integrated into the geographic information system (GIS) through the weighted overlay option in the ArcMap software by assigning the resulting weight values to the comparison matrix of AHP of each parameter.

After the conclusion of the process of weighting for each criterion (Table 11), these weights are used to correct the importance degree of each criterion by multiplying it by the map layers for each criterion using the weighted overlay option on ArcMap software (Kahsay et al. 2018), in order to extract the map of the distribution of erosion coefficient values following the equation:

$$\begin{aligned} {\text{SSEI}} & = \left( {{\text{Mean}}\;{\text{annual}}\;{\text{rainfall }} \times 0.32} \right) + \left( {{\text{Soil}}\;{\text{texture }} \times 0.22} \right) + \left( {{\text{Land}}\;{\text{use/ land}}\;{\text{cover}} \times 0.14} \right) \\ & \quad + \left( {{\text{Lanforms }} \times 0.04} \right) + \left( {{\text{Drainage}}\;{\text{density }} \times 0.04} \right) + ({\text{Drainage}}\;{\text{frequency}} \times 0.05) \\ & \quad + \left( {{\text{Lineament}}\; {\text{frequency}} \times 0.03} \right) + \left( {{\text{Slope}} \times 0.10} \right) + \left( {{\text{Relative}}\;{\text{relief}} \times 0.04} \right) \\ \end{aligned}$$

Results and discussion

A process of diagnosis and characterization of the impact of the erosion phenomenon in the Bechar watershed was carried out using recognized methods of spatial analysis. The purpose of using the USLE model is to evaluate the rate of soil loss in the region. The hierarchical analysis method (AHP) contributed to the determination of the areas most sensitive to erosion through the use of the nine environmental parameters (annual average rainfall, soil texture, land use/soil cover, reliefs, drainage density, drainage frequency, lineament frequency, slope and relative relief). Resorting to the use of AHP was due to its advantage in being able to use and handle many data and inputs at the same time in order to obtain a more accurate and clear result. Despite these positives, some researchers mentioned some of limitations of this model.Velasquez and Patrick (2013) mentioned that among the limitations of this model is self-assessment bias, which affects the internal validity of the results. Hartwich (1999) showed in detail a set of limitations for the AHP model, among which we mention that AHP does not give any constructive guidance to the structuring of the problem, because a different structure may lead to a different final ranking, bur a solution proposed by Saaty (1991) which is based on arranging the elements of the comparison in clusters to avoid extreme different between the structures of each decision maker. Some researchers criticized the judgment scale of Saaty which is based on verbal comparisons. Despite its appeal and ease of use, they are not clear in some cases, especially when you make some complex comparisons (Ishizaka and Labib 2009). Barzilai (2005) said that it is difficult to represent preferences using the ratio scale, because there is no absolute zero in such cases as is the case in comparisons related to temperature and electrical tension.

Calculation of the long-term annual average soil loss (A) and sediment delivery ratio (SDR)

The results of the application of the Arnoldus equation based on the monthly and yearly precipitation data for Bechar basin showed that the values or erosivity ranges between 4.43 and 27.5 MJ mm ha⁻¹ h⁻¹ year⁻¹. The highest value for R has been recorded in the northern part of the watershed. According to Fig. (2a), we note that the change in the value of R is changed gradually toward the southern part of the watershed where the smallest value is recorded for R. The mean value of erosivity of Bechar’s watershed is 9.32 MJ mm ha⁻¹ h⁻¹ year⁻¹. This value is considered close to the values recorded in both the Zousfana watershed 12.05 MJ mm ha⁻¹ h⁻¹ year⁻¹ (Bouzouina 2014) and the Guir watershed which is located in the west of Bechar’s watershed which is equal to 6.96 MJ mm ha⁻¹ h⁻¹ year⁻¹.

According the DSMW database, the distribution of soil types in the watershed of Bechar is divided into three types (clay, silt and silty clay) so that the clay dominates the formation of the soil by more than 70% of the basin area which is concentrated in the center and in the top part. On this basis, the value of k factor ranges between 0.020875 and 0.02276 t ha h ha⁻¹ MJ⁻¹ mm⁻¹. According to Fig. (2b), we can note that the biggest value has been recorded in north of the basin, this region contains specifically the mountains of Antar and the Horreit. These are considered to be rocky structures. The minimum value was recorded in the basin outlet. For a number of factors, perhaps the most important of them is the frequent positioning of small soil particles, which creates a soil layer with a fragile structure, easy to be eroded.

The study of the topography of a Bechar watershed showed that slope values range from 0° to 58°, more than 80% of watershed area with a slope value of less than 5°. The results of applying the Wischmeier and Smith equation showed that the value of LS of Bechar watershed ranges between 0 and 122.60 (Figs. 2c). The maximum values of LS have been recorded value on the mountains chains mainly concentrated in the northern part as well as the eastern part of the watershed (Djebel Antar, Djebel Horreit).

The results of the analysis of satellite images Landsat 7 ETM + showed that there are three classes of land cover for the watershed of Bechar (bare land, urban land and vegetal cover) (Figs. 2d). In order to know the value of c factor for each class, We have used the classification proposed by Jung et al. (2004). According to this classification, the value of c factor ranges between 0.01 and 0.35. The smallest value (0.01) has been concentrated in the north part of the watershed, and it is characterized by intense vegetal cover which is estimated as the ideal protection against soil erosion. Depending on the geographical location of the study area, which is located in a dry climate, most of the basin area is a barren land where the risk of erosion is at its maximum value.

In the watershed of Bechar, we can notice the lack of erosion control practices, most farmers depend on cereal crops and rarely have plowing parallel to the lines of contour; for that reason, we have used the table proposed by Shin and Pesaran (1999) which is based on the slope percent and farming type. The P factor map (Figs. 2e.) showed that the values of P factor in Bechar’s watershed range between 0.55 and 1. The value 1 represents the existence of efficient anti-erosion practices. Throughout the watershed, they represent only 4.35% (276.71 km²) of the total catchment area which is 6357 km².

The application of the USLE model is based on the multiplication of the factors layers. Using the Raster calculator option on ArcMap software, the result has shown that the average annual soil erosion rate varies between 0 and 4.61 t ha⁻¹ yr⁻¹ (Fig. 7), with the average value of 0.016 t ha⁻¹ yr⁻¹. The vulnerable to erosion areas mapping (Fig. 8) shows five classes distributed from the least severe in the south to the most severe north of the watershed. From Table 12, we can see that the weak class dominates the area of the watershed by more than 56% with annual erosion rate ranging between 0 and 2.77 t/ha/yr; this class is characterized by low erosivity (equal to 9.17 MJ mm/ha h, yr) and low erodibility (equal to 0.02162902 t ha h/ha MJ mm), whereas nearly 8% of the watershed area is subjected to sever and very severe soil erosion with annual erosion rate ranging between 2.77 and 4.61 t/ha/yr; this class is concentrated mainly in the north of the watershed (Djebels Antar and Horreit), which is characterized by very steep slopes and weak land cover.

Fig. 7

Average annual soil erosion

Fig. 8

Vulnerable areas classification map

The sedimentation delivery ratio (SDR) module was calculated to determine the soil fraction that contributes to the sedimentation of storage areas within the watershed. Three models have been used and subjected to a statistical treatment through the EPM module whose aim is to choose the SDR model most adapted to the watershed studied. The results indicate that the most suitable model among the three used is Vanoni’s (1975) with a fraction equal to 0.149%. The sediment yield is the result of the multiplication of the annual rate of soil loss and the SDR (Swarnkar et al. 2018). The zonal statistics option under ArcMap has been used to determine that the soil loss rate for each class gives a total soil loss in the Bechar watershed of 585,252 t/yr, i.e., a sedimentation yield of 87,202.55 t/yr. In proportion to the annual erosion rate, it can be concluded that the SDR is influenced by the watershed geo-environmental parameters where it can be distinguished that large SDR values are concentrated in steep and low-cover areas. These results were compared with the results of studies conducted on the basins surrounding the watershed of Bechar. Bouzouina et al. (2014) assessed the mean annual erosion rate of 3.64 t/ha/year in the watershed of Zousfana which is located east of the Bechar watershed, so that this watershed is characterized by a weak vegetation cover, except for some agricultural lands distributed along the length of the Zousfana mainstream, and a soil structure characterized by sand dunes that are easy to be eroded. According to Belkendil (2014), the mean annual erosion rate in the watershed of Guir is estimated of 1.73 t/ha/year, where we can notice that it features the same characteristics of the Zousfana watershed with the exception of the presence of soil structure more stable, reducing the value mean annual soil erosion.

Assessment of soil erosion risk/potential

In order to validate the results obtained, the soil sensitivity index for erosion (SSEI) was calculated using the hierarchical multi-criteria analysis (AHP) method. Nine geo-environmental parameters were included in this step.

The large area of the watershed is subject to an arid climate that is characterized by spatial variation and irregular rainfall. The northern part is characterized by greater precipitation intensity than the southern part of the basin; according to comparison matrix (Table 10), the rain intensity, which exceeds 140 mm, contributes by 48% of the weight of the erosion effect. Then, following the middle part, where the intensity varies between 100 and 120 mm and takes 32% of the weight of the precipitation effect, we can consider that this zone is an interim zone of arid to hyperarid climate. In the third part, which extends from the center to the south, the intensity is very low and does not exceed 80 mm. In this area, the erosive effect is very low and only reaches 0.7% of the influence weight.

The second parameter, in order of importance according to the AHP matrix, is the soil texture with a value of 0.22. This result seems logical because of the obvious impact of the soil composition on soil loss; clayey soils with high organic matter content are more resistant to erosion than other types of soil. Coarse soils are concentrated mainly in the northern part of the watershed in the mountains of Grouz, Horriet and Antar. Specifically, in the region called Kodiet Haidora, this area has a Syncline shape which promotes the deposition of coarse soil. We can also say that this zone is the part where solid transport begins. Medium-diameter soils occupy the majority of the basin; in this part, there is a slight slope, medium to low precipitation, and this prevents the transport of medium diameter soils. The transport of fine soils is carried out in areas located to the south where the average precipitation and the slope are very low, so it can be said that this perimeter is a spreading zone which promotes the sedimentation of fine soils due to low speed.

The land-cover/land-use parameter has a value of 0.14 which is the third influence on erosion. The arid zones are known by their low density of plant cover. Despite this, there is visual variability in the vegetal cover species. The study area is composed of three subclasses (Table 10): there is a large concentration in the north and in some agricultural lands near the rivers. The vegetal cover takes 73% of the weight of land cover (forest class) and is characterized by a high resistance to erosion relative to the other subclasses. The urban areas are distributed along the watershed, despite its density being low by its surface, and erosion phenomenon in urban land is low and contributes by 20% of the weight. Barren land occupies most of the watershed, and they are more venerable to erosion. So soil protection is very weak and takes the weight 0.7% of weighting (Tables 10, 11, 12).

Table 10 Summary of pairwise comparison matrices used for risk mapping/soil erosion potential

Table 11 Pair comparison matrix

Table 12 Classes of areas vulnerable to erosion

The hydrographic network of the Bechar watershed is very dense and varies between 0 and 9.27 km/km². The large values of hydrographic density are concentrated in the southern part; on the other hand, the northern part is characterized by a low hydrographic dense. The high drainage density implies the high drainage frequency which varies between 0 and 11 no/km². The distribution of these two parameters in the watershed strongly influences the erosion rate.

The geological formation of the watershed varies from the Paleozoic to the quaternary. The movement of these formations over the ages generates the faults that are distributed in the study area. It is easy to distinguish the fault using Landsat satellite images. The Ksiksou fault is estimated to be the largest fault in the region; this fault divides the watershed into two parts and extends from north to southwest. The assembly of faults in the north is tighter than the south which gives a form of sheet to the south platform called Chabka. Over time, the faults transform into effluents and streams, which increases the rate of erosion in the watershed.

The slope parameter, with a value of 10%, was ranked in the fourth row; according to Fig. 5a, one can distinguish that the watershed topography is flat in the majority with the exception of the northern and northwestern parts which are characterized by mountain ranges (Djebel Bechar) where the slope values reach more than 58°. This latter has a weight of 0.48%, the largest in the matrix of comparison.

The diversity of reliefs in the center of the study area is not as clear as the northern part and the eastern part of the basin where you can easily notice big reliefs, which is because of positioning of the mountain like that of Antar and Grouz. The maximum slope in the basin reaches a value of 58%, and the major part of the study area is characterized by flat lands with a gradual variation topographic from north to south (the outlet).

Figure 9 illustrates the distribution of sensitivity index erosion in the studied watershed; values vary from 0.130133 to 0.0445591. The maximum values were recorded in the northern part and precisely on the mountain ranges of Antar and Horreit Djebels; however, the minimum values were observed in its central and southern parts.

Fig. 9

Soil erosion index map

Figures 7 and 8 show that there is a high vulnerability to erosion in Antar Mountain by a value of 0.13 and less values of erosion risk in the other adjacent mountains that are located in the same area which is the north of the basin, where rainfall is high, this can be argued by the weak plant cover which prevents erosion, and the medium-sized soils are dominant in this part, as well as the great heights of relief. In the same map, we can notice a moderate vulnerability class in the region of Ouakda, which is due to high density of vegetal cover because of the existence of agricultural areas, high relief deficiency and high frequency of drainage (Horiet and Djebel Bechar). The central part is characterized by low vulnerability of erosion, though the terrains are bare, the drainage density is high and the soil size is medium. This result relied heavily on weak amounts of precipitation, which is estimated as the first launcher parameter of erosion. The southern zone is characterized by a very low vulnerability, and this part is characterized by flat lands which helps to soil deposition.

Conclusion

This study presents a the impact spatial characterization of the erosion phenomenon in the Bechar watershed, located in southwestern Algeria, and subject to an arid climate regime. The methodology consisted of calculating the annual rate of erosion by the USLE model and then evaluating the sediment yield through the EPM module. The AHP multi-criteria analysis method contributed to the determination of erosion-sensitive areas in the watershed using geo-environmental parameters that influence the erosion phenomenon in the study area. According to the USLE model, the annual erosion rate varies between 0 and 4167 t ha⁻¹ yr⁻¹ with the average annual value of 0.016 t ha⁻¹ yr⁻¹.

The Vanoni EPM module (1975) indicates that 14.9% reach the watershed outlet with a sediment yield equal to 585 252.01 t/year. According to the classification of the erosion rate final map, based on the degree of erosion intensity, it can be seen that more than 56% of the study area is subject to low erosion intensity, a rate that varies between 0.92 and 1.84 t h⁻¹ yr⁻¹.

This zone is characterized by a low rainfall (80–100 mm/year) and a flat topography which favors the flow velocity dissipation and tends toward the of soil particles deposit. A small area of the watershed (1.30%) subject to very severe erosion is concentrated in the north of the watershed and particularly in the mountain ranges that dominate in this part (Djebel Antar and Djebel Horreit).

According to the AHP pairwise comparison table, we can see the geo-environmental parameters large values of mean annual rainfall (0.32), soil texture (0.22), land use/land cover (0.14) and the slope (0.10).

The other parameters used in the study (landform, drainage density, drainage frequency, lineament frequency and relative relief) did not significantly affect erosion in the study area.

All the results quoted above present an accurate diagnosis of the situation of the Bechar watershed on the erosion phenomenon. The erosion-sensitive areas maps were made to provide support for researchers and decision makers in this area to better intervene through appropriate tools to combat soil degradation in the study area.