Water Erosion Risk Analysis in the Arribes del Duero Natural Park (Spain) Using RUSLE and GIS Techniques

Abstract

Nowadays, soil erosion is a global problem of great environmental and social concern, affecting natural resources, natural spaces and agricultural production. Therefore, it is necessary to carry out an erosion risk analysis to estimate the amount of soil lost, as well as to establish possible conservation practices to mitigate this loss. One way of doing this is through the integration of empirical equations such as RUSLE and GIS techniques, giving rise to a mapping of potential and actual erosion, considering the factors that make up this equation. The results obtained indicate that the areas with extreme erosion levels in Arribes del Duero, that is, with the greatest losses (greater than 200 Tm/ha/year), correspond to areas with steep slopes, poorly developed soils such as Leptosols and Regosols and vegetation with little or no vegetation cover. On the other hand, areas with stable levels of erosion (up to 10 Tm/ha/year) are found in flat areas, with more developed soils, such as Alisols and Luvisols, and vegetation with a higher density and herbaceous cover. Finally, it is concluded that the integration of GIS techniques with parametric equations constitutes a simple and economic tool for estimating these losses and, together with land use, allows different mitigation measures to be established, which, in our study area, focus on reducing the length and gradient of the slope, such as contour cultivation, construction of terraces and “bancales”.

1. Introduction

Soil loss through erosion is a global problem, affecting natural resources and agricultural production [1,2,3,4,5], and has increased significantly during the 20th century, becoming a global problem of major environmental and social concern [6]. Average soil erosion rates worldwide are estimated to be between 12 and 15 T ha⁻¹ yr⁻¹ [7], which means that each year, the earth’s surface loses about 0.90–0.95 mm of soil [8]. Among the climatic factors causing soil erosion, water is more influential than wind [2].

Water erosion is a natural process involving the separation, transport and sedimentation of materials [9,10,11,12]; the main triggers of which are, among others, precipitation, slope and land use changes [13,14]. It constitutes the major form of land degradation, serving as a precursor of irreversible effects on the soil, causing loss of fertility, slope instability and loss of surface horizons [15,16,17].

Human-induced erosion, on the other hand, is linked to deforestation, poor agricultural practices, overgrazing, forest fires and rapidly increasing urbanisation. The above-mentioned factors, together with other inappropriate land practices, are responsible for triggering erosion along with other inappropriate land management practices [18,19,20]. In general, these practices can reduce the productive potential of agricultural regions, generate slope instability and reduce soil porosity, which lead to loss of water retention, infiltration and percolation capacity. As a result, surface runoff, sediment transport, siltation and water pollution due to the transport of agrochemicals such as fertilisers and pesticides increase [13,18,19].

Soil erosion may become more severe in the near future due to climate change, further aggravated by increased population pressure, overexploitation of natural resources and poor land and water management practices [20,21]. Soil conservation is needed to reverse the process of land abandonment and improve agricultural production to ensure food security and sustainability. Therefore, there is a need to identify critical areas prone to erosion to provide the necessary information to establish soil conservation strategies, such as protected areas [22].

Researchers have developed different tools to estimate soil loss empirically, such as the Soil and Water Assessment Tool (SWAT), the Water Erosion Prediction Project (WEPP), the Universal Soil Loss Equation (USLE) and the Revised Universal Soil Loss Equation (RUSLE) [2,23,24]. Among these models, the USLE and its revised version (RUSLE) are the most widely used due to their simplicity, ease of use and ability to successfully integrate the various parameters of the ecosystem or natural areas [25,26]. Therefore, the RUSLE model will be used in this study.

RUSLE is an erosion prediction model, which estimates long-term annual soil loss with acceptable accuracy [27,28,29]. This model comprises the following five factors: rainfall erosivity (R factor), soil erodibility (K factor), topographic factor (LS factor), land use (C factor) and soil management and conservation practices (P factor) [30,31].

In addition to using empirical equations such as RUSLE in recent years, Geographic Information Systems (GIS) and Remote Sensing (RS) have become useful tools for natural resource management and disaster research. The use of these technologies, which facilitate the handling of many spatial data in a fast and efficient way [32], allows the generation of erosion risk models and cartographies through the analysis of a database to elaborate classifications, map algebra, etc. In this database, it is possible to integrate all the basic thematic parameters (R, K, LS, C and P), resulting in specific mappings, which are overlaid to finally establish synthetic erosive risk mapping [33]. For this reason, many researchers use GIS as the main approach to estimate soil erosion at all scales [34,35,36,37,38,39,40,41].

In this article, which was carried out in the Arribes del Duero natural park (Salamanca-Zamora), the objective is to determine the risk of soil erosion (not previously studied) by using the RUSLE model and, in addition to this, the use of GIS and RS, highlighting as a novelty with respect to other studies the use of satellite images of the highest current relevance and digital terrain models with high spatial resolution, allowing, in a quick way, the development of a cartography of erosion risk in a quick way to distinguish potential and actual erosion. In this natural park, the conservation of soil resources is a priority in order to promote economic activities that lead to population settlement. Finally, this cartography and the land use will determine the possible conservation practices or measures that could be implemented in the future to mitigate these losses.

2. Materials and Methods

2.1. Study Area

The study area (Figure 1) chosen for this work is the Arribes del Duero Natural Space, located to the west of the provinces of Salamanca and Zamora, on the border with Portugal. It is a protected area of 1061 km², made up of 38 municipalities and a population of about 17,000 inhabitants. The climate is characterised by mild winters and very hot and long summers in the valley areas, with an average annual temperature and rainfall of 17.1 °C and 500 mm, respectively, in contrast to the extreme continental climate that characterises the plains, with temperatures of 12.2 °C and rainfall of 750 mm [42]. Its landscape is characterised by an undulating peneplain (with a uniform height of 700–800 m) and the steep slopes formed by the canyons (with heights of 130 m) carved by the river system (Duero, Tormes, Uces, Huebra and Águeda rivers). In terms of vegetation, the “peneplain” is a rich mosaic, delimited by stone walls and pastureland, with species of the genus Quercus (holm oak, pyrenean oak, cork oak and gall oak), mixed with other tree species (ash trees) and scrub (scrubland and broom), pastures and dry crops (wheat, barley, rye and vines). For their part, olive and almond trees remain on the slopes and terraces, only displaced by myrtle, holm oak and juniper groves, where agricultural use has been abandoned [43,44]. Finally, it should also be noted that this is one of the areas with the greatest hydroelectric potential in the Iberian Peninsula.

2.2. Methodology

The methodology followed in this article combines field work with laboratory work. As a result, a series of water erosion risk maps of the soils of the Arribes del Duero Natural Park was obtained. The field work focused on obtaining representative samples of the different types of soils existing in the study area. The laboratory work considered of analysing the samples, taking care to establish the necessary parameters for calculating the different factors (granular-metric analysis, organic matter, structure, etc.) that determine the application of the RUSLE to the risk of water erosion. Finally, the data obtained in the field campaigns and the laboratory analyses were studied by applying graphic (Wischmeier nomogram, DTM generation, etc.) and empirical procedures (formulas for the calculation of parameters, RUSLE equation, etc.), in order to elaborate a database that has been implemented in a Geographic Information System (ArcGis 10.5) and to obtain different parametric and final erosion risk cartographies of the study area.

Therefore, the quantification of soil losses due to water erosion has been carried out by means of two mappings (Figure 2) [45].

1. Map of Potential Erosion or states of the terrain under hypothetical natural conditions: this constitutes the susceptibility of an area to erosion. In order to predict this risk, a detailed study of a series of factors or elements of the physical environment (mechanical resistance, rainfall, slopes, etc.) that condition erosion processes is performed. Knowing these variables, the potential erosion units can be inventoried and mapped by using erodibility indices (lithofacies and slopes) and erosivity indices (aggressiveness of rain).

2. Current Erosion Map: it considers current conditions and determines the degree of current soil loss in each area by considering the “present moment” and analysing the soil forming and protective factors, as well as its spatial distribution (types of crops and native vegetation and conservation practices).

For both mappings, the modified version of the RUSLE has been used to estimate the average annual soil loss under different conditions of use, climatic variation, relief and use of conservation practices. This model is expressed by Equation (1):

A = R × K × LS × C × P

(1)

where A is the soil loss per unit area in a given time in Tm/ha/year, R is the rainfall erosivity factor, K is the soil erodibility factor, LS is the topographic factor encompassing slope length (L) and slope steepness (S), C is the land use and management factor and P is the soil conservation practices factor.

2.3. Potential Erosion Risk Cartography

The Potential Erosion Map is obtained from the multiplication of the three factors R, K and LS. Once obtained, the reclassification of the erosive degrees is carried out in smaller intervals, using the criteria in

Table 1. Annual soil loss classification.

Types of Erosion	Loss of Soil (Tm/ha/year)
Very low erosion and tolerable soil loss	<5
Low erosion and tolerable soil losses	5.1–10
Mild erosion level	10.1–25
Moderate erosion level	25.1–50
Severe erosion level	50.1–100
Very severe erosion level	100.1–200
Extreme erosive level	>200

.

2.3.1. Rain Erosivity Factor (R)

This factor considers the average of the kinetic energy intensity values estimated from monthly and annual average rainfall [12]. This rainfall must have a continuous record of rainfall intensity variations and has been obtained from the database of the Geographic Information System for Agricultural Data (SIGA) [46].

To obtain this factor, firstly, rainfall data were collected from 35 stations in the study area, with data from more than 20 consecutive years. Then, given the dispersion of the stations, an interpolation was carried out using the weighted distance method (IDW) with ArcGis. In our study, we have worked with a cell size of 5 metres, thus obtaining greater precision. This operation was repeated each of the month of the year, thus obtaining the corresponding raster. The next step was the application of the Modified Fourier Index (IMF) [47], whose formula is (2):

$$IMF={\Sigma }_{i=1}^{12} \frac{p_i^2}{Pt}$$

(2)

where p_i is the precipitation of each month in mm and Pt is the mean annual precipitation in mm. Each of the previous raster corresponds to each of the 12 values of p_i, so that from them we can obtain the raster that represents the value of the mean annual precipitation, which will be the sum of all the monthly precipitation values. Finally, to calculate the R Factor, we have used Equation (3), which is the regression equation proposed by the I.C.O.N.A. for the region where our study area is located, which allowed us to calculate the value of R as a function of readily available precipitation variables, such as the total precipitation or maximum precipitation in a month [48]:

$$R=2.56\times IM{F}^{1.065}$$

(3)

2.3.2. Soil Erodibility Factor (K)

This is a sensitive parameter related to regional characteristics, soil structure and degree of weathering. It describes the susceptibility of the soil to detachment and transport of particles in quantity and flow rate for a specific predicted rainfall event. It is a quantitative value experimentally determined from soil texture, structure, organic matter content and permeability [10,12].

This factor is obtained using with the data obtained in the physico-chemical analyses carried out in the laboratory of the 38 soil profiles taken in the field work, necessary for the use of the corresponding Equation (4), whose values are shown in

Table 2. Erodibility values.

Soil Type	K Value
Gleyc luvisols	0.12
Chromic alisols	0.17
Chromic cambisols	0.19
Dystric gleysols	0.20
Eutric cambisols	0.22
Dystric cambisols	0.24
Eutric regosols	0.28
Dystric regosols	0.38
Dystric leptosols	0.39
Lithic leptosols	0.49

. The values obtained have been validated with the Wischmeier nomogram [12].

$$100K=\left[{10}^{-4}\times 2.71\times T^{1.14}\times \left(12-MO\right)\right]+4.2\times \left(E-2\right)+3.2\times \left(P-3\right)$$

(4)

where T is the texture parameter of the surface 15 cm, MO is the organic matter content (%), E is the structure and P is the permeability.

2.3.3. Topographic Factor (LS)

The topographic factor is determined by the length of the slope (L) and the slope (S). Thus, it establishes the need for an exhaustive knowledge of the spatial distribution, since soil erosion is intensified as a consequence of the concentration of runoff water towards the lower areas. As L increases, soil erosion increases, as well as increasing as a consequence of velocity and surface runoff [12].

To calculate this factor, the methodology followed by Zhanh et al., 2013 [49] was used, based on Equation (5), proposed by Moore and Burch (1986) [50]:

$$LS={\left(Flow accumulaion\times cell\frac{size}{22.14}\right)}^{0.14}\times {\left(\frac{sinslope}{0.0896}\right)}^{1.3}$$

(5)

where Flow Accumulation is the number of cells contributing to the flow in a given cell, cell size is the length of the size of one side of the cells and sin slope is the sine of the slope in radians.

First, the slope is calculated in degrees and then transformed to radians. Next, the calculation of the slope length is carried out, which is done on the basis of the flow accumulation raster, which represents the cells in which water accumulates when flowing from the cells with the highest altitude value. To obtain it, it is necessary to calculate the filling of sinkholes, the flow directions and the accumulation of flow. Once the slope and slope length are obtained, the LS Factor is calculated considering the above formula, which involves multiplying the flow accumulation by the cell size. In our study, the maximum flow accumulation value is 4796.821, which multiplied by the cell size (5 m) would result in a high maximum runoff length. Therefore, using these data would be erroneous, as it would overestimate the value of the slope length, resulting in exaggerated erosion values. To avoid this, it is necessary to establish a maximum length of 25 m, which is equivalent to 5 cells, thus obtaining a flow accumulation raster with a maximum value of 5. Finally, a reclassification of the flow accumulation is carried out (considering this maximum length), which will be used to calculate the LS Factor by means of the above formula.

2.4. Real Erosion Risk Cartography

This cartography is obtained from the Potential Erosion Map by adding two terms from the RUSLE, which are: the crop or vegetation factor (C-factor) and the Conservation Practices Factor (P-Factor).

2.4.1. Plant Cover Factor (C)

This factor analyses the influence of plant species, crop rotation and the degree of erosive susceptibility of the soil, which will influence its productivity. For its calculation, the management of plant masses and crops is considered, using the Forestry Map of Spain, scale 1:25,000, delimiting the study area. Finally, to obtain this value (

Table 3. Vegetation Cover Values.

Vegetation Cover Type	C Value
Mixed hardwood forests	0.003
Mediterranean scrub	0.04
Riparian forest	0.09
Poplar and banana plantations in production	0.09
Ash groves	0.09
Wild olive groves	0.18
Juniper groves	0.18
Cork oak groves	0.19
Meadows	0.19
Holm oak groves	0.19
Oak groves	0.19
Chestnut groves	0.22
Non wooded	0.24
Mixed conifers	0.42

), the values established for tree, shrub and mixed tree formations are considered, analysing the percentage of tree and shrub cover, type of herbaceous cover and thickness of plant debris as well as its extent [22]. For the herbaceous formations, the Wischmeier classification was considered [51]

2.4.2. Soil Conservation Practices Factor (P)

This factor analyses the existence of soil conservation practices in land use [27]. In this study, as in most soil loss estimates, it is not considered, given that we are interested in knowing the potential and real losses considering natural factors. On the other hand, human activities can increase this soil loss or reduce it through the implementation of specific actions: terraces and contour cultivation. As a result, this is a parameter to be subtracted from the potential erosion risk, so it is not considered, i.e., Factor P in our study area has a value of 1 [44].

3. Results and Discussions

The USLE factors (Figure 3), have been estimated with a 5 m × 5 m grid, being a very precise scale that allows us to provide the special distribution of the average annual soil erosion in the Arribes del Duero Natural Park with great detail and to resemble reality as much reality as possible, unlike other grids of lower precision.

3.1. RUSLE Factor Analysis

The R Factor has been obtained from the monthly and annual rainfall data of the different meteorological stations spread over the study area with normal data (20 years) by applying the IMF. Average annual rainfall ranged between 483.9 and 846.4 mm. The highest rainfall was concentrated between Saucelle and Aldeadávila. This factor shows values ranging between 156 and 291 (Figure 3A), being maximum between the above-mentioned localities.

The K factor has considered the soils studied by the authors in previous works, as well as their corresponding physico-chemical data [52]. With these data and the Wischmeier nomogram, values between 0.12 and 0.49 have been obtained (Figure 3B).

The topographic factor is influenced both by the length of the slope and by the gradient; therefore, geomorphological aspects are one of the main factors that determine the emission of sediments in a river basin. Our study area is characterised by the boxing in of the Duero River, i.e., a peneplain followed by a steep slope, the latter being more susceptible to erosion. Observing the results obtained (Figure 3C), the values range between 0 and 61. The areas with the steepest slope and the longest slope have the highest values.

Finally, the vegetation cover factor (Factor C) is a determining natural element in the protection of the soil against the erosive force of precipitation, because in addition to controlling the energy with which raindrops hit the surface of the soil, it slows down the speed of surface runoff. Additionally, stoniness must be considered because it acts as protection, reducing the inertia with which the drops fall, since the soil is covered by fragments of rock or gravel. In view of the results obtained (Figure 3D), C values between 0.003 and 0.42 are obtained.

3.2. Potential Erosion Risk Cartography

To make this mapping (Figure 4), the factors of the physical environment have been multiplied: R, K and LS, showing the susceptibility of the area to erosion, considering the existing conditions. In this way, the areas corresponding to the Duero River basin, as well as its main tributaries such as the Tormes, Águeda, Huebra and Uces, are where the greatest erosion is identified, with values of over 200 Tm/ha/year (reaching 5137.36 Tm/ha/year), being identified as an extreme erosive level.

This area coincides with high slopes that favour turbulent movements, giving rise to severe erosive forms, such as rills (Figure 5A) and gullies (Figure 5B), and, in addition, with poorly developed soils, such as Leptosols and Regosols, which have a higher sand content that act as erosive agents as it is dragged along by surface runoff. There are differences between these soils: the lithic and distric Leptosols are more susceptible than the Eutrophic regosols. The former because they have an A horizon with little thickness (<10 cm); the latter because they are desaturated. The clays are more susceptible to dispersion, compared to the Eutrophic regosols, which are saturated. On the other hand, the areas with a moderate and medium level of erosion, with values between 10.1 and 50 Tm/ha/year, correspond to peneplain areas in which medium development soils predominate, such as Gleysols and Distric and Eutric cambisols, which, unlike the previous ones, have a cambic. Among them, the Distric cambisols are the most susceptible as they are desaturated, compared to the Eutric cambisols and the Distric gleysols. Finally, the tolerable levels of erosion (up to 10 Tm/ha/year) are also found in peneplain areas, with the difference that the soils are more developed, such as Chromic cambisols, Gleyic luvisols and Chromic alisols. These soils are characterised by a higher content of dehydrated iron oxides, which act as a bridge between the clay and the organic matter, giving this complex greater stability to the soil. The other two types of soil are characterised by a Bt horizon of clay illuviation, i.e., with a higher percentage of clay, and also have a higher organic matter content with a more developed clay-humic complex, which protects the soil from possible erosion.

3.3. Real Erosion Risk Cartography

This map indicates the current erosion risk (Figure 6). It has been drawn up on the basis of the potential erosion risk considering the vegetation cover factor, which gives the soil certain protection, depending on its characteristics (height, density, percentage of cover and territorial extension). In order to classify the degrees of water erosion obtained, a reclassification is carried out considering the aforementioned criteria, expressed in Tm/ha/year and mm/year.

In view of the results, the areas suffering the greatest water erosion are those located in the river beds, especially in the case of the Duero River and its most abundant tributaries. These sectors are characterised by their steep slopes, increasing the speed of surface runoff as it descends, thus carrying away the materials most susceptible to erosion. In addition to this, it coincides with areas of vegetation with little protective power and with low percentages of cover, such as conifers and broadleaf trees, which further accentuates the erosive vulnerability of these soils. These areas show severe to extreme erosion, with values ranging from 50.1 to >200 Tm/ha/year and soil losses ranging from 4.18 to >16.70 mm/year. On the other hand, we find lower erosion values in the plain areas (<4.18 mm/year), as there are no or little slopes and the characteristic vegetation provides greater protection. It has a greater density and herbaceous cover and makes this area less vulnerable to erosion.

3.4. Erosion Validation and Mitigation: Land Uses Erosion

For the validation of erosion, the real erosion cartography and the land use map (SIOSE) have been considered, making it possible to determine the uses most susceptible to erosion in order to establish possible conservation practices to mitigate these losses. Thus, the most affected uses, ordered from highest to lowest, are as follows (Figure 7): coniferous and hardwood forests, crops, olive trees, non-citrus fruit trees and scrubland. In the case of coniferous forests, the needles tend to acidify the ground when they fall, preventing the development of herbaceous plants and leaving on the ground bare and exposed to runoff.

The hardwood forests characteristic of the study area does not have an appreciable undergrowth and, therefore, leave the soil exposed to erosion. Both olive trees and non-citrus fruit trees have a very low percentage of cover with hardly any protection for the soil. Lastly but to a lesser extent, scrubland, varies according to its density; if it is denser, it will provide greater protection to the soil.

Bearing this in mind and considering that the most affected areas have a steep slope, the practices that can be carried out are: reducing their length (and thus, the speed of runoff) and breaking up the slope. Some examples of such practices are contour cultivation (following contour lines), construction of “bancales” (Figure 8A) or terraces (Figure 8B), among others. The latter have already been carried out in the area, especially with vine crops. As they do not require highly developed soils they are adapted to this type of practice. Therefore, the use of this crop and this practice is very useful for reducing erosion in these areas where the soils are not very developed and the slope is steep.

4. Conclusions

Currently, there are numerous methodologies for calculating the losses caused by water erosion in soils, most of them based on contrasted formulas but difficult to update and implement quickly taking into account current geomatic resources (high-resolution orthophotos, metric and centimetric lidar models, use of UAVs, etc.). Thus, as a novelty in this article, USLE has been implemented by means of factorial or parametric analysis in GIS, making it possible to generate an automated RUSLE model so that, with new and future technologies, erosion values can be improved in great detail. In this way, taking into account the determining components of the territory, two maps are obtained, one of potential erosion and the other of actual erosion, from which soil losses due to water erosion in the study area are calculated and quantified.

In the study area, three zones are differentiated according to their degree of erosion. The first ones are zones with an extreme level of erosion, losses of more than 200 Tm/ha/year, high slopes, poorly developed soils and vegetation with little protective power, low density and cover, such as conifers and broadleaf trees. The second ones are areas with moderate and medium erosion levels, with values between 10.1 and 50 Tm/ha/year, corresponding to the peneplain, with medium developed soils and vegetation with a certain degree of protection. The last ones are areas with tolerable levels of erosion (up to 10 Tm/ha/year) also coincide with penillanura but with more developed soils. The vegetation provides greater protection, as it has a higher density and herbaceous cover and makes this area less vulnerable to erosion.

In turn, these maps make it possible to establish in a simple way the degrees of erosion established according to the FAO, expressed in Tm/ha/year and mm/year and, together with the land use map, they constitute a low-cost non-structural measure that helps to identify the areas where it is necessary and urgent to implement conservation practices that mitigate soil losses. With this in mind, these measures will focus on reducing the length and steepness of the slope such as contour cultivation, “bancales” or terracing.