Abstract

Soil loss triggered by water erosion constitutes a major issue that causes several environmental and socioeconomic concerns. The Moroccan Western High Atlas is the most vulnerable area in the High Atlas Mountains, due to the existence of different forms of landslides, and evidences of erosion are widely observed. This study aims at estimating and quantifying the amount of soil loss as well as highlighting potential areas to erosion risk, using the Revised Universal Soil Loss Equation (RUSLE) combined with GIS and remote sensing. The RUSLE model provides a possibility of computing erosion susceptibility for each pixel on the basis of the controlling factors which are rainfall aggressivity, topography, vegetation cover, soil erodibility, and support practices. In this study, results show that the erosion rate varies between 0 and 227.67 t/ha/year, with an average annual soil loss of 40.38 t/ha/year, and the Beni Mohand River basin is subject to very high rates of erosion which can be irreversible since it exceeds the tolerable standard rate which is 1 t/ha/year. These findings will provide land use planners baseline for land use and risk management and will provide data within the Moroccan Western High Atlas Mountains.

1. Introduction

Soil erosion is a widespread phenomenon that affects the environment and the economy all over the world [1]. The average rate of soil loss in the world is estimated to be approximately 12 to 15 t/ha/year; this implies a loss of 0.90–0.95 mm every year of topsoil [2], which is one of the most affected parts of the soil by the erosion process [3, 4]. It has been estimated that 11 million km2 is affected by erosion caused by water erosion [5]. Furthermore, climate change and intense agricultural practices are among the problems that cause soil erosion [6]. As a consequence, these problems contribute to fertility loss, quality degradation of soil resources [7], and dam siltation [8], which eventually endanger food production, leading to degradation of water quality as well as eutrophication in water bodies [9].

Soil loss has always existed, but it has rapidly grown in recent decades. Soil erosion causes extensive damage every year globally and is the reason why numerous researches have been conducted in different areas worldwide [1013] to improve methods and techniques used for the study of land degradation. On a basin scale, various models have been developed for the assessment of soil loss risk, namely, the EUROpean Soil Erosion Model (EUROSEM) [14], Limburg Soil Erosion Model (LISEM) [15], Soil and Water Assessment Tool (SWAT) model [16], and Water Erosion Prediction Project (WEPP) model [17].

The use of GIS and remote sensing for erosion assessment has proved to be a reliable tool [18], when combined with empirical/semiempirical models [19, 20]. The USLE model of Wischmeier and Smith [21] is the most commonly used empirical model, thanks to its simple function [22, 23]. This model is calculated on the basis of the Universal soil loss Equation (USLE), adapted to the Moroccan conditions by Arnoldus [24], and becomes the Revised Universal Soil Loss Equation (RUSLE). This equation uses five factors for estimating soil erosion, which are rainfall erosivity, soil erodibility, topography, vegetation cover, and soil conservation practices [25]. These factors have been computed to provide a soil loss map, where the erosion is calculated for each pixel. The resulting map delivers a spatial distribution of eroded areas as well as quantification of soil loss with higher accuracy [26].

In the Moroccan Western High Atlas, especially in the Beni Mohand River basin, it has been observed that extensive land degradation has been caused by landslides, soil erosion on the mountain slopes and along river banks. This is due to friable lithology, the existence of blocks of rocks, the absence of vegetation, and high-intensity rains. The aim of this study is to estimate the soil loss using the RUSLE model in order to highlight and identify areas of high erosion risk, which will help decision-makers to implement corrective measures.

2. Study Area

The watershed of Beni Mohand River is part of Taroudant-Oulad Teima piedmont in the Moroccan south Atlas. It covers 34,894 ha, and it is located in the southwest of Morocco between Taroudant and Agadir city (Figure 1), exactly between Issen River on the west and El Ouaar River on the east.

A morphostructural study investigated by Ait Hssaine [27] shows two distinguished units. The first is the ancient massif made up of very resistant Georgian limestones that arm narrow anticlines and Ordovician shales bent into synclines. These rocks rest on the tuffs and lavas of the Infracambrian and the whole is strongly tectonized. A second is the Sub-Atlas Zone, consisting mainly of Lower and Middle Cretaceous and Eocene limestones and marls, topped by pink sandstones and upper Eocene clays.

According to the geological studies undertaken in the plain of Souss [28], but particularly the part of Beni Mohand River [29], it is characterized by Georgian schist-limestones, Acadian sandstone, quartzite, and Ordovician shale. The dominant materials are limestone pebbles of different sizes especially at the alluvial cones which occur in the region, where some cones between Issen River and Beni Mohand River can reach 5–6 km in length. The little amount of soil existing between the blocks allows some sparse seasonal vegetation to get established, which disappears by mid-June, leaving the soil surface bare.

The general climate is of the arid to semiarid type, with the annual rainfall ranging from 250 to 300 mm/year in the lower elevations to 500–600 mm/year in mountainous areas [30]. The average annual temperature is between 14°C and 18°C but can reach 20°C in low elevations [31]. The rainfall season is from November to March, and the aridity period can be up to 7 months (April to October). In Taroudant region, the aridity extends from mid-January to the end of October. The study area is characterized by a mild winter [29]; the first winter rains are violent and cause erosion by runoff, leading to tiny elements and exposing the bare skeleton of the ground [29].

3. Materials and Methods

In the current study, the cartographic method used is the USLE model of Wischmeier and Smith [21], revised and adapted to the Moroccan conditions by Arnoldus [24], and then called RUSLE. In previous studies, this model has shown remarkable flexibility for the available data and efficiency of the resulting map.

The RUSLE model was calculated in a GIS Environment [32] based on the universal soil loss equation, as a product of five factors that control soil erosion (equation (1)). These factors are calculated separately using GIS and remote sensing and combined to compute an annual soil loss map that predicts the rate of erosion. It also helps to estimate the amount of deposition in the entire area [33].

The RUSLE model is expressed by the following equation [21]:where A is the average annual soil loss rate (t ha−1 yr−1), R is the rainfall erosivity factor (MJ mm ha−1 h−1, yr−1), K is the soil erodibility factor (t h MJ−1 mm−1), LS is the topographic factor (L in m, S in %), C is the land cover factor, and P is the conservation practice factor.

Figure 2 shows the adopted methodology in this study.

3.1. Rainfall Erosivity Factor (R)

The R factor or the rainfall erosivity factor describes the aggressiveness of rain to cause erosion in a given location [34]. It reflects the sensitivity to erosion based on the amount and intensity of rainfall. The rainfall erosivity factor in RUSLE uses the kinetic energy over 30 minutes at the maximum intensity [21, 35]. However, due to the lack of datasets for Morocco, we opted for using available datasets for monthly and annual rainfall records [24] for over 40 years (1970–2010). This alternative approach has been applied to the 19 meteorological stations in and around the study area (Figure 3). The R factor has been estimated by the following equation developed by Arnoldus [24]:where R is the average rainfall erosivity (MJ.mm/ha.H.year), MRi is the average monthly precipitation in (mm), and AR is the average annual precipitation in (mm).

The point data of rainfall erosivity are interpolated using the Inverse Distance Weighted (IDW) method, as a deterministic interpolation method computed based on the influence of the distance for the climate station locations to generate a rainfall erosivity map in the catchment. The IDW has been tested and used to calculate the erosivity factor and has been successfully useful [36]. The result of the R factor represents a spatial distribution of rainfall aggressivity for the entire area.

3.2. Soil Erodibility Factor (K)

Soil erodibility refers to the susceptibility of soil particles to the detachment under the impact of raindrops, runoff, or both, taking into account the morphological, chemical, physical, and mineralogical characteristics of the soils [37]. The K factor describes quantitatively the inherent erodibility of each soil type. Soil erodibility can be influenced by soil structure, soil permeability, soil texture, soil profile, and organic matter. It refers to the amount of eroded soil per unit of erosive energy of rainfall, with 9% slope and 22.13 m long and considering a plot of clean bare soil [30].

To calculate the soil erodibility K, several formulas have been developed to meet this need. However, because of the absence/inaccessibility to some required datasets-soil permeability and soil structure datasets as mentioned in the original formula [21] or even the inadequacy of other formulas to the local conditions of the study area, we opted for the equation developed by Sharpley and Williams [38] (equation (3)) which has been successfully used in a similar environment [39].

The used soil data in this study has been derived from the universal soil database to create the erodibility map since the collection of soil data of K factor from the field is expensive and time-consuming:wherewhere SAN, SIL, and CLA are the percentage of sand, silt, and clay, respectively, C is the organic carbon content, and SN1 is the sand subtracted sand content of 1 and divided by 100.

3.3. Topographic Factor (LS)

The topographic factor called also Slope Length and Steepness Factor (LS factor) plays an important role in the erosion process and indicates the effects of topography on erosion, referring to the length (L) and steepness (S) of the slope that influence the surface runoff speed [40]. In other words, the LS factor considers that the slope gradient and slope length influence the production and transport of deposits [41]. The erosion process increases as the slope length increases [30]. Moreover, it can be increased promptly with the increase of the slope steepness. The vegetation cover and soil particle size are two parameters that influence the relationship of soil loss to the slope gradient.

Both slope length and slope gradient were extracted from the Alos DEM (12.5 m resolution) in a GIS Environment and then combined to generate a topographic factor map using the following relation (equation (5)) of Wischmeier and Smith [21], developed by Bizuwerk et al. [42]:where S is the slope gradient (%) and L is the length of the slope (m). L = flow accumulation × DEM spatial resolution and the value of “m” as given in Table 1.

3.4. Cover Management Factor (C)

After the topography, vegetation cover factor is the second most important parameter that triggers erosion. Vegetation cover intercepts the effect of rainfall and increases the infiltration and permeability and eventually the speed of rock alteration. According to Wischmeier and Smith [21], C factor differentiates between vegetation type and density. It also distinguishes between bare and covered lands.

In this study, the resulting C factor map is derived from a Normalized Difference Vegetation Index (NDVI) (equation (6)) of the Landsat 8 OLI image (resolution of 30 m). This image has been preprocessed radiometrically and atmospherically. The rasterized map of NDVI was then classified into five classes that represent five land uses [21], where the area of each type was calculated separately (Table 2). Results were validated with the help of the ground truth verification. For each type of land use, an adequate C value is assigned. In this study area, C factor values range from 0.05 to 0.7 as shown in Table 2. Higher values stand for areas with no vegetation cover while low values indicate areas with vegetation cover [43]:where the NIR is the near-infrared band, which refers to surface reflectance values of band 5 (Landsat8 OLI), and R is the red band, which refers to the surface reflectance values of band 4 (Landsat8 OLI).

3.5. Support Practice Factor (P)

The support practice or the soil conservation practice factor indicates the ratio of soil loss assigned to each cultivation method depending on its ability to conserve the soil from erosion, such as contour plowing, ridging, slope terracing, and alternating strip crops [44]. This factor is estimated based on the relationship between slope and agriculture practices. As these datasets are not available, we opted for a method developed by Shin [45]. This method has been conducted in a similar environment where each type of conservation practices has a relationship with a particular slope interval, which means that the P factor has been estimated depending on slope variations. The P factor values range between 0 and 1, the value of 0 corresponds to areas with anthropic erosion resistance, and the value of 1 has been designated to areas with steep slopes and indicates the absence of anthropic support practices (Table 3).

3.6. Potential Erosion Map

Within a GIS Environment, all the thematic maps of factors R, K, LS, C, and P were processed. The resulting maps of 30 m resolution were combined using the empirical equation of the RUSLE model, to produce a soil loss map of the area. The soil loss map will identify areas with a high rate of erosion as well as a quantitative description of vulnerable areas.

4. Results and Discussion

4.1. Rainfall Erosivity Factor (R)

The thematic map of rainfall erosivity factor was computed based on meteorological datasets, located in and around the watershed. These datasets contain the annual and monthly average rain for the 19 stations as described in Figure 3. The erosivity map (Figure 4) has shown a spatial distribution of the R factor all along the watershed. Low to moderate values occupy the center of the watershed. However, high to very high values are observed in the watershed borders, which means that the rainfall aggressivity increases from the middle area to the edges.

The rainfall erosivity factor values range between 35.65 and 44.56 MJ mm/ha.H.year. The highest values are located in the South-Western and North-Eastern parts, representing downstream and upstream areas of the studied catchment, respectively.

For a better representation of the R factor’s intensity, we use the standard deviation classification. Values between 39.18 and 44.56 MJ mm/ha H year represent high values of erosivity, situated in the upstream and downstream parts of the watershed. This indicates that 75% of the area is at risk caused by rain, 18% suffers from a moderate rate of erosion, and only 7% of the area is exposed to low erosion. These results have shown that the studied area is subjected to a significant rainfall erosivity.

4.2. Soil Erodibility Factor (K)

The soil erodibility varies according to climatic variability, soil type, and cultivated methods [46]. In this study area, the K factor varies from 0.0176 to 0.0227 t ha H/ha MJ mm (Figure 5), the highest value representing the class of 0.0227 t ha H/ha MJ mm; this class stands for soils with high vulnerability to erosion and it occupies 13.74% of the entire watershed and occurs in the downstream. The value of 0.0183 represents areas with medium susceptibility to erosion risk, occupying 25% of the total area and located in the upstream of the watershed. The areas with low erodibility are in the center of the basin, this category represents 62%, and it is regarded as almost the dominated category.

4.3. Topographic Factor (LS)

The topographic factor constitutes one of the major components of erosion risk studies from the slope length/slope gradient point of view. Other parameters can also influence this factor such as soil type and lithology [47]: hard rock often protects steep slopes and resistant floors while soft rock gives fragile soils on more gentle slopes; this explains that erosion often arises where the slopes are moderate.

According to the calculation of the LS factor, the values obtained vary between 0 and 664.57 (Figure 6). These results are grouped into four classes using standard deviation classification, where 94% of the watersheds were characterized by low slopes, the majority of which are located in the southern part of the basin. Moderate slopes occupy 2.11% and 3.61% of the entire area represented by steep slopes. The moderate to very high classes are corresponding to medium to high elevation and occur in the upper part of the basin area. This means that the upper part (upstream) is subject to a high risk of erosion in comparison to the downstream areas.

4.4. Cover Management Factor (C)

At Beni Mohand River basin, four categories of dominant land uses are bare lands, cultivated lands, lands with moderate vegetation, and lands with dense vegetation. Land use with high C factor value contributes significantly to soil erosion. The cultivated lands are limited to few areas, particularly in the downstream areas where the elevation is lower. The agricultural lands help to fix the ground through cultivation to a certain extent; meanwhile, it enhances erosion by disturbing the texture and organic matter of the soil. The moderate and dense vegetation areas are spread in the upper part of the watershed. However, the bare lands are scattered throughout the Beni Mohand basin.

The C factor values of Beni Mohand River basin vary between 0.05 and 0.7 (Figure 7). The values of 0.05, 0.10, 0.6, and 0.7 correspond to the following land uses, respectively: dense vegetation, moderate vegetation, cultivated land, and bare land. The dominant land use class is bare land covering more than three-quarters of the entire watershed, followed by lands with moderate vegetation cover (12.98%), lands with dense vegetation (6.42%), and at last the cultivated lands (2.12%).

4.5. Support Practice Factor (P)

The support practice factor is the ratio of soil loss with specific support practice on agricultural land to the corresponding loss with parallel slope plowing [21]. In other terms, the P factor values vary according to the type of agriculture applied and slope. The conservation practices are extremely useful in reducing soil loss in areas with steep slopes and where the erosion rate is high. In this area, no significant conservation practices have been taken to reduce the soil loss rate apart from cultivated lands, where some species have been planted such as citrus, argan, and olive trees.

In this study, P values were estimated based on slope values. The soil conservation practice factor varies between 0.55 and 1 (Figure 8), where the higher the value, the lower the opportunity to the soils to be supported, and as a consequence a high risk of erosion. High values correspond to areas of high slopes and vice versa. Values varying between 0.55 and 0.6 represent areas with support practices, or at least areas with a low slope that do not need any support. The value of 0.8 corresponds to areas with moderate slopes. The values of 0.9 to 1 stand for areas without any support practices; these areas represent more than 60% of the surface of the Beni Mohand River basin. Moderate to high P values represent zones with high erosion susceptibility.

4.6. Potential Erosion Map

The erosion map has been prepared by the combination of RUSLE factors, which are climatic aggressivity R (or rainfall erosivity), soil erodibility K, the combined effect of gradient and length slope of LS factor, vegetation cover C, and soil conservation practices P. The integration and processing in a GIS Environment of these five factors led to the development of a map indicating the potential erosion at Beni Mohand River basin. This map allows us to have a spatial distribution of erosive risk due to the effect of these natural factors. The obtained potential soil erosion map (Figure 9) ranges from 0 to 227.67 t/ha/year, with an average of 40.38 and 1.4 million tons of losses every year.

The resulting map has been classified into six classes using FAO’s [48] classification to describe soil loss intensity [49] and visualize both susceptible areas to erosion and areas of deposition; these classes are as follows (Figure 10):(i)The first class groups areas with low potential erosion risk under 5 t/ha/year. It represents 84.5% of the surface of the studied area and it is spread along the watershed particularly in the southern part.(ii)The second one is a class of areas with moderate erosion risk, between 5 and 25 t/ha/year. This class represents 13% of the total surface.The third class gathers areas with medium potential erosion risk between 25 and 50 t/ha/year. It constitutes 2.3% of the studied basin surface. Moderate to medium classes are focused mainly on the upper part of the watershed.(iii)The fourth category varies between 50 and 200 t/ha/year. This class represents areas with high to very high potential erosion risk, and it covers 0.21% of the studied area.(iv)The last class represents a class of critic potential erosion risk of more than 200 t/ha/year. It constitutes only 0.06% of the total area.

The soil loss map of Beni Mohand River basin has successfully highlighted and identified the areas exposed to erosion, the majority of which are focused on areas with moderate to high gradient slopes, generally in the upstream area of the watershed. In addition to the topographic factor, areas of greater rainfall, bare land, and friable soil or even rock blocks involve more vulnerability to water erosion. These areas are mostly located in the upper part of the watershed, where the estimated soil erosion is considered to be moderate to high.

This result acquired from the RUSLE model reveals the current or potential alteration of eroded materials and informs about the amount of soil loss provided by erosion in the watershed. With the lack of datasets in the Moroccan Atlas Mountains, the RUSLE model remains effective and powerful as it gives decision-makers the possibility of anticipating erosion risk in order to combat or at least minimize this risk, particularly in priority areas where the risk is high. The validation of the obtained results is fundamentally required for model reliability, although it is challenging especially with the lack of datasets in the studied area. However, some methods remain useful as they are used as an alternative data source.

Indeed, the validation methods used in this study are as follows. The first method is the ground truth which is a field verification of sampled locations from the soil loss map crosschecked with ground truth data. 80% of the checked points matched the ground truth, which means that these eroded areas have been clearly identified and only 20% of these points are areas where erosion can be potentially induced as mentioned before. The second method is the comparison with studies done before in areas having almost the same characteristics, especially those of the Atlas chain [5052] which have shown the high performance using the RUSLE method despite its limitations; it constitutes a useful tool to identify priority areas to be preserved from the erosion risk.

According to the FAO’s report [2], soil erosion in Morocco affects up to 40% of its territory with the total annual soil loss evaluated at 100 million tons, and this is equivalent to 50 million m3 annual reduction in dam storage capacity. This cumulative annual erosion is reaching an alarming level; that is why it becomes a necessity to look for concrete solutions. Indeed, this study aims to provide a significant tool to mitigate not only the amount of erosion but also various problems and consequences related to soil loss that are as follows:(i)Natural risk and vulnerability: erosion destabilizes slopes and initiates landslides or collapsed blocks which can damage villages and infrastructure downstream.(ii)Environmental: erosion processes remove the fertile part of the soils and thus reduce the effective depth to be exploited by roots and the amount of water available to plants. This is considered a major constraint limiting productivity in Morocco [53]. In addition, soil erosion can change the vegetation distribution and the stability of watersheds.(iii) Hydrological: the region is considered to be semiarid and has soil erosion in causing dams siltation which consequently reduces their longevity. Recently, the Moroccan government policy is oriented towards the construction of hill dams for water retention and the supply of groundwater.(iv) Socioeconomic: in this part of Morocco, the rural exodus constitutes a scourge endangering the stability and development of neighboring urban areas due to soil loss for cultivated areas in such watersheds. This study suggests some practical solutions in order to manage these agricultural areas and at the same time minimize the impact of soil erosion, which will undoubtedly contribute to reduce the rural exodus and thus the pressure of the urban areas.

It is now obvious that soil loss represents a real danger for population, environment, and infrastructure for which all countries including Morocco must adopt strategies to ensure their stability. In addition, this document constitutes a road map for stakeholders for future projects as it provides (1) a detailed study with relevant results and (2) solutions to minimize the soil loss damage such as increasing levels of C in the topsoil using conservation measures [54], developing the vegetation cover, prohibiting deforestations, and managing watersheds using bleachers and streams deviation.

5. Conclusions

Soil erosion is a common phenomenon that harmfully impacts the ground. This study mainly deals with the prediction of potential soil erosion areas using the RUSLE method and GIS techniques. It provides a concrete estimation of the distribution of soil erosion at Beni Mohand River basin, throughout the combination of the influencing factors such as rainfall aggressivity, soil erodibility, topography, vegetation cover, and soil conservation practices.

The soil loss estimated at Beni Mohand River basin varies from 0 to 227.67 t/ha/year, with an annual loss of 1.4 million tons/year. Results have shown that 15% of the surface area displays an important rate of soil loss, and it represents zones with moderate to severe soil erosion. These areas are mainly situated in moderate to very high slope where the runoff is important. This loss is favored by the other factors of erosion, which are also combined to accelerate erosion, significant losses (93% of the total area has a very important rainfall erosivity), moderately erodible soils (49% of the soils show a K factor between 0.0183 and 0.0227 t ha H/ha MJ mm), 78% of the bare land, and 6% of the areas representing steep slopes. These statistics show that the Beni Mohand River basin is subject to a medium to high erosion risk and even severe to some specific areas.

The universal soil loss equation combined with GIS and remote sensing plays an important role in producing the soil loss map, starting from raw data to the final result that provides a detailed assessment of erosion risk. This method is widely used since creating a database of the soil erosion of the ground is costly difficult and time-consuming. Hence, the resulting map can be more reliable while using high-resolution spatial data.

It is difficult to stop the soil erosion risk. However, it could be reduced with a suitable land use management and adequate support practices which can fix the topsoil in the region. Yet, the managers and planners can use the soil loss map in order to conserve the areas of priority and preserve or manage potential areas to the erosion risk.

Data Availability

To evaluate soil loss in the Mohand River basin, the data used are as follows. (1) The ALOS (12.5 m) digital elevation model (DEM) was used to retrieve the slopes, slopes’ lengths, and slopes’ orientation. Every result of this DEM was converted to a 30 m resolution raster map as all datasets in this study and it is available at https://search.asf.alaska.edu. (2) Historical data of precipitation (1970–2010) were collected from Agence du Bassin Hydraulique de Souss Massa (ABHSM), which included a time series of the monthly precipitations of different climatological stations in and around the studied area. (3) The soil data were extracted from the Universal Database of the soil of FAO, and it is used to fill the lack of soil data in the area, particularly data of soil type, organic component, and texture. It is freely downloaded via http://www.fao.org/land-water/databases-and-software/en/. (4) The Landsat 8 satellite image was used to calculate the Normalized Difference Vegetation Index which helps to extract land use and vegetation cover map of the study area. This image was acquired on 19 July 2017. LC08 imagery was captured using a multispectral imaging sensor providing 11 bands, with a UTM/WGS84 projection (zone 29) and 30 m of resolution; it is a 16-bit encoding and freely downloaded from GloVis USGS at https://glovis.usgs.gov/app?fullscreen=1. (5) Ground truth data were based on High-Resolution Google Earth to validate areas representing a risk of erosion.

Conflicts of Interest

The authors declare that they have no conflicts of interest.