Received 6 December 2015; accepted 26 February 2016; published 29* February 2016
Soil erosion represents one of the most serious land degradation problems  . It is defined as the loosening, dissolving and removal of rock materials from all parts of the earth’s surface triggered by a complex interaction process of many factors: natural (climate, topography, soil, vegetation) and anthropogenic (tillage systems, soil conservation measures, overgrazing and deforestation)   .
In the north of Morocco, soil erosion affects negatively agricultural productivity, reduces water infiltration, underground water resources and water availability.
In Morocco, 40% of land is affected by water erosion  . In some parts of the Rif in northern Morocco, erosion rates sometimes reach 30 to 60 t∙ha−1∙y−1   . Therefore, dams lose their water capacity initial storage due to their siltation which is estimated to 0.5% by year  . The largest Moroccan dams receive each year approximately 50 million tones of sediment  , which affects their storage capacity and brings about an annual loss of almost 300 million Dirhams  .
Due to the intensification of agricultural practices leading to unsustainable farming practices (e.g. inappropriate tillage practices, straw exportation, overgrazing) and specific bioclimatic conditions (e.g. recurring and severe droughts), more than 15 million hectares of the Moroccan agricultural land is under serious threat. As reported by Namr and Mrabet  , it was estimated that out of 22.7 million hectares potentially exploitable in the Northern part environmental problem worldwide of Morocco, 77% were exposed to very high erosion risks  . The Global Assessment of Human Induced Soil Degradation (GLASOD) survey carried out during the 1980’s by the United Nations Environment Programme (UNEP) and the International Soil Reference and Information Centre (ISRIC) established that the severity of human induced degradation had been classified as severe and very severe for more than 20% of the Moroccan territory  .
In the Rif Chain (Northern Morocco), several local studies have been applied to evaluate and quantify the erosion risk, especially in the western and central Rif   .
Therefore, the estimation of erosion factors and exposed areas to soil erosion can be very helpful to identify the increment and the degree of the risks and, finally, to establish conservation measures and soil/water management plans.
In the current study, an effort to predict potential annual soil losses in Oued El Makhazine watershed in Morocco has been conducted using the Universal Soil Loss Equation (USLE) adopted in a GIS environment. It is used for the prediction of sheet erosion depending on the distribution of the aggressiveness of rainfall, the erodibility of soil, topography, land use practices and crop management.
2. Study Area
Located in the north of morocco, Oued El Makhazine watershed (Figure 1) covers an area of 2414 Km2, and stretches from north latitude 34˚45'15.3" to 35˚15'44.6" west longitude 5˚14'32.4" to 5˚51'9.1", its altitude is between 10 m and 1677 m.
Figure 1. Oued El Makhazine watershed location map.
The study area consists west of plains with a very pronounced topography, while to the east the relief becomes more hilly and mountainous. The altitude increases gradually eastward with the first reliefs of the Rif mountain range, this configuration of the relief associated with the presence of the Atlantic Ocean to the west is the main cause of significant rainfall measured in the watershed, the study area receives an annual average rainfall of 1100 mm  .
The predominance of impermeable soils in the watershed promotes runoff, which is increased by the slope effect as one progresses eastward.
The vegetation cover of the watershed is defined by the presence of matorral, typical Mediterranean landscapes, with a predominance of forests. Most of the forest cover is at the upper watershed. Downstream, the agriculture is the dominant economic activity in the watershed consists of cultivated plains.
3.1. Data Processing and USLE Factors Generation
In recent years, GIS and remote sensing are often used to assess and map water erosion effects. Truthfully, with these modern techniques, it is increasingly exposed the advantages of spatialization methods for assessing and mapping soil erosion over large areas and setting up scenarios for rehabilitation.
Many techniques and studies realized worldwide were done about the evaluation of soil loss. Most of them are using the Universal Soil Loss Equation (USLE) and its revised version (RUSLE)  . Others had modified part of the equation to adapt in every country’s situation. In this study, we tried to promote the process using the potentials of GIS.
The USLE equation is a product of five input factors (Figure 2) in raster data format: soil erodibility; rainfall erosivity; slope length and steepness; cover management; and support practice. Each factor varies over time and space and depend on other input data. Therefore, soil erosion within each pixel was calculated using the Universal Soil Loss Equation (USLE).
The USLE equation is described as:
where A is the computed spatial average of soil loss over a period selected for R, usually on yearly basis (t∙ha−1 ∙y−1); R is the rainfall erosivity factor (MJ∙mm∙ha−1∙h−1∙y−1); K is the soil erodibility factor (t∙ha∙h∙ha−1∙MJ−1∙mm−1);
Figure 2. Flowchart of the methodology.
LS is the slope length steepness factor (dimensionless); C is the cover management factor (dimensionless, ranging between 0 and 1); and P is practices factor (dimensionless, ranging between 0 and 1).
In this study, the Universal Soil Loss Equation (USLE) was combined with GIS technologies to estimate the potential soil loss from areas within the Oued El Makhazine watershed, generate soil erosion severity maps, and analyze areas of critical soil erosion conditions which claim urgent need for convenient conservation measures and land management.
3.1.1. Rainfall Erosivity Factor (R)
The erosivity factor (R) in the USLE equation is accounted for by the rainfall-runoff, it is considered as a driver of soil erosion processes. The R factor represent the effect of raindrop impact and also shows the amount and rate of runoff associated with precipitation events, it is defined as the product of kinetic energy and the maximum
30 minute intensity and shows the erosivity of rainfall events  .
Rainfall data of 10 years (2004-2014) was obtained from the Loukkos Hydraulic Basin Agency and imported into GIS environment since all the weather stations had geographical coordinates. Annual and monthly rainfall data of Oued El Makhazine watershed was used to calculate the R-factor in this study.
The rainfall erosivity values for the different stations were used to interpolate a rainfall erosivity surface using the Inverse Distance Weighted (IDW) interpolation method in ArcGis environment to generate a raster map for R factor. The IDW interpolation method was selected because rainfall erosivity sample points are weighted during interpolation such that the influence of rainfall erosivity is most significant at the measured point and decreases as distance increases away from the point. The IDW interpolation method is based on the assumption that the estimated value of a point is influenced more by nearby known points than those farther away   .
The Equation (2), below developed by Wischmeier and Smith  and modified by Arnoldus  was used in the computation:
where R is the rainfall erosivity factor (MJ∙mm∙ha−1∙h−1∙y−1), Pi is the monthly rainfall (mm), and P is the annual rainfall (mm).
The mean values of R factor range from 74.63 MJ∙mm∙ha−1∙h−1∙y−1 to 116.11 MJ∙mm∙ha−1∙h−1∙y−1 (Figure 3(a)), approximately 76.82% of the rainfall erosivity factor R ranges between 74.63 MJ∙mm∙ha−1∙h−1∙y−1 and 96.59 MJ∙mm∙ha−1∙h−1∙y−1.
3.1.2. Soil Erodibility Factor (K)
The soil erodibility factor (K) relate to the average long-term soil and soil profile response to the erosive power associated with rainfall and runoff. It is also considered to represent the rate of soil loss per unit of rainfall erosion index for a specific soil. The USLE model utilizes the technique proposed by Wischmeier  to measure the K factor of a soil type.
The specific methodology concerns many properties of soil such as gain size (silt, clay and sand), organic matter content, soil structure and soil permeability   .
The soil samples collected from Oued El Makhazine watershed were analyzed in soil laboratory of National Institute of Agricultural Research of Settat. For each sample the particle size analysis was performed in five classes (sand, very fine sand, silt, very fine silt, and clay) using Robinson’s pipette method  . Organic carbon was measured by the Walkley-Black wet dichromate oxidation  and converted to organic matter by multiplying it by 1.724. Soil permeability code was obtained from the National Soils Handbook No. 430 (USDA, 1983)  and Soil structure code was determined using Soil Textural Pyramid  .
To calculate the soil erodibility, the formula of Wischmeier and Smith  was used.
where K is the soil erodibility factor (t∙ha∙h∙ha−1∙MJ−1∙mm−1), M is the particles percentage (% of very fine sand + % of silt x (100% clay)), a is the organic matter content (% C x 1.724), b is the soil structure and c is soil permeability.
The analysis showed that 72 % of samples had clay texture, 18% Sandy clay loam and 10% Clay loam texture.
Erodibility is low for clay rich soils with a low shrink swell capacity, as clay particles mass together into larger aggregates that resist detachment and transportation  .
Organic matter content varies from the lowest of 1.37 % up to the highest of 4.57%. The permeability of Oued El Makhazine watershed is very low, a result that was observed in numerous publications  - that the presence of clay decreases the level of soil permeability.
The result (Figure 3(b)) shows that Oued El Makhazine watershed has soil erodibility, ranges from 0.24 to 0.85, the watershed has a moderate to severe soil erodibility with 71.26% of the area ranges between 0.68 t∙ha∙h∙ha−1∙MJ−1∙mm−1 and 0.85 t∙ha∙h∙ha−1∙MJ−1∙mm−1.
3.1.3. Slope Length and Steepness Factor (LS)
The L and S factors in USLE reflect the effect of topography on erosion. It has been demonstrated that increases in slope length and slope steepness can produce higher overland flow higher erosion  . Moreover, gross soil loss is considerably more sensitive to changes in slope steepness than to changes in slope length  .
Therefore one of the key factors in soil loss is topography, especially when the ground slope exceeds a critical angle. Topography factor plays a major role in soil erosion since it dominates the surface run off rate.
The combined LS factor of the watershed was derived from Digital elevation model by means of ArcGis environment. The computation of LS requires factors such as flow accumulation and slope steepness.
The flow accumulation and slope steepness were computed from the DEM using ArcGIS Spatial analyst extension.
The Equation (4), was adopted to calculate LS factor as proposed by Moore and Burch   .
where flow accumulation denotes the accumulated upslope contributing area for a given cell, LS = combined slope length and slope steepness factor, cell size = size of grid cell (for this study 30 m) and sin slope = slope degree value in sin.
The LS factor value in the study area varies from 0 to 22 (Figure 4(a)). The Majority of the study area (77.72 %) has LS value ranges between 0 and 6.29.
3.1.4. Cover Management Factor (C)
The coverage factor ranges from 0 (full soil coverage) to 1(no soil coverage)  .
C factor represents the effect of soil disturbing activities, plants, crop sequence and productivity level, soil cover and subsurface bio-mass on soil erosion  . Natural vegetation plays a predominant role in reducing water erosion  .
Vegetation indices such as the NDVI (Normalized Difference Vegetation Index) are quantitative measures, based on vegetation spectral properties that attempt to measure biomass or vegetative vigor  . As an indirect estimate of vegetative density, a Normalised Difference Vegetation Index (NDVI), which approximates chlorophyll density, was calculated for the study area using Landsat TM images with a spatial resolution of 30 meter.
Satellite imagery acquired during the rain season are more adapted for this application given that soil erosion is most active and vegetation cover is at its peak during this season.
where NIR is light intensity in the near-infrared, and RED is light intensity in the red band.
This index is an indicator of the energy reflected by the Earth related to various cover type conditions. When the measured spectral response of the earth surface is very similar to both bands, the NDVI values will approach zero. A large difference between the two bands results in NDVI values at the extremes of the data range.
Photosynthetically active vegetation presents a high reflectance in the near IR portion of the spectrum (0.76 - 0.90 μm), in comparison with the visible portion Red band (0.63 - 0.69 μm); therefore, NDVI values for photosynthetically active vegetation will be positive. Areas with or without low vegetative cover (such as bare soil, urban areas), as well as inactive vegetation (unhealthy plants) will usually display NDVI values fluctuating between −0.1 and +0.1. Clouds and water bodies will give negative or zero values.
The following formula was used to generate a C factor surface (Figure 4(b)) from NDVI values  - :
where a and b are unitless parameters that determine the shape of the curve relating to NDVI and C factor. According to Van der Knijff an a-value of 1 and a b-value of 2 seem to give reasonable results  .
3.1.5. Conservation Practice Factor (P)
The P factor is the soil-loss ratio with a specific support practice to the corresponding soil loss with up and down slope tillage  . Usually, the conservation practice factor corrects the USLE estimation for management and tillage practices that protect the soil from erosion.
In the present study the P factor map was derived from the land use/land cover and support factors. The values of P factor ranges from 0.55 to 1 (Figure 4(c)), in which the high value is assigned to areas with no conservation practices; the minimum values correspond to built-up-land and plantation area with strip and contour cropping.
4. Results and Discussion
Average Annual Soil Loss
The data layers (maps) extracted for K, LS, R, C, and P factors of the USLE model were integrated using Equation (1) within the raster calculator option of the ArcGIS spatial analyst in order to quantify, evaluate, and generate the maps of soil erosion risk and severity for Oued El Makhazine watershed.
Generally, if the estimated (A) value is high, it means a higher rate of sediment yield, while a lower value denotes a lower rate of sediment yield  .
The Oued El Makhazine watershed was classified into five soil erosion risk categories (Figure 4(d)). The area and proportion of soil erosion risk classes are illustrated in Table 1. The maximum and minimum losses are respectively about 0 t∙ha−1∙y−1 and 735 t∙ha−1∙y−1, about 65.26% (1575 km2) of the watershed ranges between 0 t∙ha−1∙y−1and 95 t∙ha−1∙y−1.
The results shows that the Oued El Makhazine watershed is exposed to a very high erosion risk, the highest
Figure 3. (a) Erosivity factor, (b) Soil erodibility factor.
(a) (b)(c) (d)
Figure 4. (a) Topography factor, (b) Land cover factor, (c) Conservation practice factor, (d) Annual soil loss.
Table 1. Classification of soil loss of Oued El Makhzine watershed according to area.
soil loss values are spatially correlated with the steepest slopes.
USLE is a straightforward and empirically based model that has the ability to predict long term average annual rate of soil erosion on slopes using data on rainfall pattern, soil type, topography, crop system and management practices. In the present research, annual soil erosion rate map is generated for Oued El Makhazine watershed, a mountainous area, which represents the first reliefs of the Rif mountain chain. Several data sources are used for the generation of USLE model input factors and are stored as raster GIS layers in ArcGis environment.
The combination of the model USLE with GIS techniques has shown divers advantages, specifically the resulats of each factor involved in the erosion process. GIS allows rational management of a multitude of data, with respect to the various factors responsible of land degradation; in fact, in this case, it allows us to conclude that the primary factor responsible of the degradation of Oued El makhazine watershed is the steep and the slope morphologies, the soil erodibility and the vegetation cover comes in the second place. GIS also simplifies and facilitates the enrichment continuous updating of the database.
The intensiveness of the rainfall coupled with steep gradient slopes causes severe erosion runoff in the study area. The result of this high runoff and soil detachment is considered as the main agent for the high rate of soil erosion at Oued El Makhzine watershed.
There is an urgency to take into consideration this soil loss by all possible ways so as to reduce the existing amount of soil loss and to raise watershed rehabilitation and productivity.
Finally, the present investigation has demonstrated that GIS techniques are low cost tools and simple for modeling and mapping soil erosion, with the purpose of evaluating erosion potential risk for Oued El Makhazine watershed.