The effect of spatial resolution on identifying vulnerable areas for sheet and rill erosion using the RUSLE and MUSLE erosion model in ArcGIS: a case study (Puentes catchment, southeast Spain)

(1)

Bachelor thesis

The effect of spatial resolution on identifying

vulnerable areas for sheet and rill erosion using the

RUSLE and MUSLE erosion model in ArcGIS: a case

study (Puentes catchment, southeast Spain)

Maurits Valentijn Kruisheer

Student Number: 10608850

Future Planet Studies (major: Aardwetenschappen)

Supervised by: Dr. Erik Cammeraat

Universiteit van Amsterdam

July 25, 2016

(2)

2

Abstract

This research studies the effect of spatial resolution on the identification of areas vulnerable for sheet and rill erosion. This effect is investigated by comparing high-resolution DEM data (from DGPS measurements) with low-resolution DEM data (remote sensing 25-meter resolution) using two often applied erosion models, the Revised Universal Soil Loss Equation (RUSLE) and the Modified Universal Soil Loss Equation (MUSLE). Patterns of sheet and rill erosion are found to differ significantly

between the two studied resolutions and erosion models. Low-resolution data fails to incorporate distinctive topographical characteristics and is therefore incapable of correctly displaying erosion patterns, whereas high-resolution data computes discrete patterns of erosion and sediment yield. Areas with low vegetation and high topographic fluctuations are found to be most prone to sheet and rill erosion, where sediment yield is more concentrated in downstream areas and erosion is most profound in the hilly upstream areas. Erosion and sediment yield rates are found inaccurate due to the absence of a spatial and temporal resolution adjusted calibration factor, which

introduction is recommended in future research. However as erosion patterns are not affected by this, this research reveals that a frequently used method like remote sensing with a 25-meter resolution is inadequate for identifying erosion vulnerable areas in an area that is heavily subjected to geomorphic processes.

(3)

3

1. Introduction

Soil erosion and land degradation are widespread phenomena in the dryland zones of the Mediterranean. Long periods of drought alternated by short lasting high intensity rainfall events render the region excessively prone to water erosion. These short rainstorms are highly influential on the region’s geomorphic dynamics (Cammeraat, 2004; Baartman et al., 2011). In these areas, water erosion originates when rainfall exceeds the water infiltration capacity of the soil, causing a flow of water and particles in the topsoil (Lesschen, 2008). Further downhill, this flow can become stronger increasing its velocity and erosive force leading to the formation of rill and gully erosion. Erosion and runoff can have serious environmental and economic implications on the landscape, affecting the region through on-site and off-site effects. Loss of topsoil and soil fertility due to water erosion severely alters the productivity of agricultural areas in the catchment, while downstream areas are affected by accelerated reservoir sedimentation rates and enhanced occurrence and severity of flash floods (De Vente et al., 2013).

In order to mitigate these erosion effects, it is important to identify areas that are most vulnerable for erosion and to comprehend the accompanying rates of sediment loss. Erosion modelling is an often used method to do this. The accuracy of erosion models is depends on the quality of the input data and the spatial and temporal scales that are used, as they determine the erosion process that is the most dominant in the model (Lesschen, Schoorl & Cammeraat, 2009). Digital Elevation Models (DEMs) are key spatial input datasets in erosion modelling as they account for information on flow directions and accumulation, slope length and steepness among other topographical components (Moore et al., 1991 and Wu et al., 2008). Spatial resolution is an important factor for displaying this information and is defined in this research as the grid cell size in hydrological models.

Several studies have demonstrated the influence of DEM resolution on the outputs of hydrological models (Chaplot, 2005; Wang et al., 2015). In general, it is assumed that a higher spatial resolution improves the realism of the model’s predictive ability due to an increased homogeneity of

topographical characteristics within the same area (Farajalla and Vieux, 1995; Mulla, 2013). Coarse spatial resolution fail to incorporate small scale topographical features which can seriously affect modelling results (Vaze, Teng & Spencer 2010). Nowadays, most erosion modelling studies are performed in spatial resolutions obtained by remote sensing data with spatial resolutions mostly ranging between 5 to 100 meters (Savage et al., 2016). Although this enables hydrological modelling to be performed on far larger scales, the quality of model outcomes can become debatable (De Vente et al., 2013; Lin et al., 2013). On the other hand, other studies have also shown that the highest possible resolution does not necessarily generate the most accurate model outcomes (Baily et al., 2014; Gillin et al., 2015). Besides costs related to on-site and off-site effects of erosion and runoff, methods for increasing DEM accuracy also increase the cost of data acquisition, data processing and storage requirements (Robinson, 1994). Therefore, a critical assessment on the validity of erosion models at different spatial resolutions is needed to guarantee and improve the applicability of modelling techniques.

Most previously performed assessments on the influence of spatial resolution neglect two important aspects. First, there is a primary focus on detecting anomalies in erosion rates rather than

differences in erosion patterns. Misjudging erosion prone areas can have large implications on the effectiveness of local environmental management plans (De Vente et al., 2013). Second, comparison studies on spatial resolutions mainly are performed only on different remote sensing spatial

resolutions (Savage et al., 2016). Even high spatial resolutions of 1 meter largely fail to incorporate distinctive topographical characteristics such as gullies, incised terraces or ephemeral flow pathways

(5)

5 which can significantly alter the model outcomes (Vaze, Teng & Spencer, 2013). To address this problem, this research uses a high-resolution DEM entailing very small-scale topographical features constructed by Differential Global Positioning System (DGPS) altitude data. In this study, the high-resolution DEM is compared with a remote sensed 25 meter high-resolution DEM in modelling sheet and rill erosion using two often applied erosion models, the Revised Universal Soil Loss Equation (RUSLE) and the Modified Universal Soil Loss Equation (MUSLE) (Lal 2001; Bhattarai & Dutta, 2007). This comparison will reveal and quantify the effect of spatial resolution upon the identification of erosion vulnerable areas and is researched by performing a case study in southeast Spain. Accordingly, the main research question is stated as follows:

What are the effects of using different spatial resolutions on identifying vulnerable areas for sheet and rill erosion on a hillslope scale in southeast using the RUSLE and MUSLE erosion models?

First, the different rates and patterns between the high and low resolution model results are compared in order to address the influence of spatial resolution on the model outcomes. Second, differences between the RUSLE and the MUSLE method are observed to incorporate the influence of different models on modelling erosion at various spatial resolutions. Third, an analysis is made on all model geographical and topographical characteristics and outcomes from which an assessment is made on which areas are most likely to be affected by sheet and rill erosion processes. Hence, the subquestions of this research are stated as follows:

 How do sheet and rill erosion rates and patterns differ between high and low spatial

resolutions when using both the RUSLE and MUSLE erosion models?

 What differences can be observed in modelling sheet and rill erosion when comparing the

RUSLE and MUSLE erosion models?

 What geographical and topographical characteristics are most vulnerable for sheet and rill

erosion in this study area?

2. Methods

2. 1. Study area

The study area is situated in the Puentes catchment, located in the southeast of Spain in the province of Murcia (Figure 1). The study site comprises an abandoned terrace system and a part of its surrounding upstream area, consisting of semi-natural vegetation and fallow agricultural lands. The presence of steep slopes in combination with short lasting torrential rainfall events trigger soil and water erosion processes in this region (Solé Benet, 2006). A long history of anthropogenic influences in the form of agricultural intensification and deforestation further increased the region’s susceptibility to erosion as low vegetation cover and marginally developed soils became more abundant. These predominantly shallow, calcaric, silty textured soils are readily subjected to crusting, sealing and compaction influencing the soil infiltration capacity and runoff dynamics (De Wit, 2001).

(6)

6

Figure 1: Study area location

2.2. DEM construction

Two DEMs with different spatial resolutions are constructed; one DEM derived from a high

resolution altitude survey using a DGPS device and one DEM derived from 25m LiDAR data. A DGPS-device from TOPCON is used to manually measure altitude data at centimetres precision as multiple satellites accurately depict the coordinates of each point measurement. Subsequently, the position from the Rover relative to the Base Station is used to calculate the Z-coordinate at each

measurement. Previous altitude measurements of the terrace system by Scholten (2015) were combined with newly gathered data from the surrounding upstream areas. The density of point measurements increases along with the degree of sudden varying topography in the landscape in order to cover a high extent of topographical features such as gullies, rills and ephemeral flow paths. This data is then interpolated using the “Topo-To-Raster” tool in ArcGIS 10.1 to construct a DEM. This tool is designed to create hydrologically correct DEMs able to follow abrupt changes in terrain such as cliffs and ridges, which are readily present in the fieldwork area. The chosen cell size of the high-resolution DEM is 0,1 meter. The second DEM is obtained from a 25m resolution aerial photo series taken in 2009 (Instituto Geográfico Nacional, 2009) In order to allow the DEM to function properly, sinks are filled using the “Fill” tool in ArcGIS.

2.3. Erosion modelling

The Universal Soil Loss Equation (USLE) and its revised and modified version (RUSLE and MUSLE) are most extensively used in erosion modelling studies due to their relatively simple applicability and low input data requirements (Lal 2001; Bhattarai & Dutta, 2007). The RUSLE and MUSLE evolved from the USLE and were developed to significantly improve the functionality and universality of the erosion model (Wischmeier, 1976). Both models predict distinctively different processes: The RUSLE is developed to calculate the annual soil loss per hectare per year (erosion), while the MUSLE is used to calculate the sediment yield after a single rainfall event. Sediment yield is defined as the amount

(7)

7 of eroded soil that is delivered to the stream network and is depicted by calculating the amount of erosion from slopes, channels and mass wasting minus the sediments that are deposited on the way before reaching the watershed flow outlet. Thus, sediment yield incorporates both erosion and sediment deposition, which is why it cannot be used interchangeably with erosion.

2.3.1. RUSLE

The Revised Universal Soil Loss Equation (RUSLE) is a simple equation calculating annual soil loss in t ha-1_yr-1_{(A) based on five factors: a rainfall erosivity factor (R), a soil erodibility factor (E), a cover} crop management factor (C), a slope length and steepness factor (LS) and a support practice factor (P) (Figure 2a-e).

The rainfall erosivity factor is derived from a research by Sanchez & Senent-Aparicio (2016) which used weather data from a station located within 2 kilometres from the study site. Using data on the kinetic energy from storms, duration and the maximum 30-min rainfall intensity averaged over twenty years (1993-2013), the average rainfall erosivity factor is calculated as 657.629 MJ mm ha-1_h -1_yr-1_{. Although Cammeraat (2004) has pointed out that local weather differences can have a}

significant impact on the erosion rates of high intensity rainfall event, a uniform R-factor is used due to the small size of the study area. The K-factor represents both susceptibility of soil to erosion and the rate of runoff, and is measured on the terrace system by Scholten (2016) using data on soil texture, soil structure and organic matter content (K-factor = 0,067 tonnes ha h ha-1_MJ-1_mm-1_{). The} K-factor of the remaining areas was determined by predicting soil texture from conducting a simple field survey, which was used as an indicator for K-factors throughout the area (Wischmeier &

Mannering, 1969). The C-factor reflects the effect of cropping and management practices on erosion rates and is derived from Panagos (2015). The LS factor reflects the effect of topography on soil erosion; a product of the slope length and slope steepness. It can be calculated in ArcGIS from the flow accumulation, cell size and the slope (Moore and Burch (1986) (Figure 3):

The flow accumulation is measured using the Flow Accumulation Tool in Spatial Analyst. The slope is measured in degrees using the Slope Tool and converted to radians in order to function properly in ArcGIS. The P-factor is the impact of support practices on the average annual erosion rate. This is kept at 1 as field specific data on management practices was not available (Zhang et al., 2009).

(8)

8

Figure 2a-e: The RUSLE input maps: Rainfall erosivity factor (a), soil erodibility factor (b), cover crop management factor (c), slope length and steepness factor (d) and the support practice factor (e).

(9)

9 2.3.2. MUSLE

The Modified Universal Soil Loss Equation (MUSLE) is similar to the RUSLE; however the rainfall erosivity factor is replaced by runoff volume (Q) and peak flow rate (qp) factors. In this research, the general MUSLE equation is used as it was developed by Williams (1977):

Where Y is the sediment yield added to the stream network in metric tonnes, Q is expressed in m3 and qp in m3_s-1_{. All other factors are similar to those used in the RUSLE.}

The runoff volume is calculated using the curve number method: Mass rainfall is converted to mass runoff by using a runoff curve number (CN) based on plant cover, interception, surface storage and amount of impervious areas (NRCS, 1986). The following equation is used to calculate runoff:

Where Qd is the runoff depth in mm, R is the event rainfall in mm and S is the potential retention parameter. The latter is based on the hydrological soil group, land cover and hydrological condition of the soil according to regulations of the NRCS (1986). Subsequently, the curve numbers are adjusted according to soil moisture conditions; which where dry throughout the study area. The potential retention parameter is calculated by:

Thus, a low curve number results in low runoff volume rates (Figure 4a-b). Lastly, the runoff depth is multiplied by the cell area (m2_{) to derive the runoff volume (Q).}

The peak flow rate is calculated by the following equation:

Where qu is the unit peak discharge in m3_s-1_km-2_mm-1_{, A is the drainage area in km}2_{, Qd is the runoff} depth in mm and Fp is the percentage of pond and swamp area over the watershed area (value of 1 used as a default). The unit peak discharge is dependent on the time of concentration (Tc); the time that is needed from each point to flow to the watershed outlet (Figure 4c). Time of concentration is the product of the travel time of three different flows: sheet flow, shallow concentrated flow and channel flow. Decades of research have led to a clear agreement of when these flows alternate each other (NRCS, 1986). First, sheet flow persists for a flow path up until 300 ft. (i.e. 91,4 meters). After this, the soil becomes saturated and flow continues as concentrated shallow flow. Following shallow concentrated flow, flow may continue in natural channels. In this research, channel flow is assumed where 5% of the maximum flow accumulation occurs. First, map algebra is used to make three binary maps indicating where each of these three flows occurs. Second, these maps are multiplied by each unique mathematical equation:

(10)

10 Where n is Manning’s roughness value, based on the land cover according to Chow (1959). L is the flow length in meters, P2 is the 24-h rainfall in mm with a 2-year return period, which is set at 48 mm based on maximum daily precipitation data from the Puentes reservoir between 1998 and 2010 (Dykes, Mulligan & Wainwright, 2015). S is the slope in m/m and R is the hydraulic radius of

channels, which is assumed 0,4 ft. (i.e. 0,12192 m). S is calculated by converting the slope in degrees to radius and then using the tangent:

The sum of the concentration times from all three flows gives the total time of concentration (Tc):

Subsequently, the unit peak discharge is derived from Tc and coefficients from the Urban Hydrology for Small Watersheds manual (NRCS, 1986). Dividing the event rainfall R by 0.2S accordingly, gives the coefficients: C0 = 1.83842, C1 = -0.25543, C2 = -0,02597. Lastly, these coefficients are used to calculate the unit peak discharge:

The peak flow is then used to calculate the total sediment yield after a rainfall event of 48 mm (Figure 4d).

2.4. Interpreting results

As high intensity rainstorms are found to be essential drivers of erosion through overland flow, erosion patterns of sheet and rill erosion should largely resemble each other in both models (Duvert et al., 2012). However, the impact of spatial resolution on incorporating topographical features in the erosion models can cause significant changes in displaying erosion vulnerable areas. Besides differences likely to be observable between spatial resolutions of the same erosion model,

differences between different models of the same resolution can be enlarged as a consequence of spatial resolution. Critically analysing the differences between both models and both resolutions at various scales is done in order to include major and minor anomalies. This includes testing

correlation between model outcomes and each depending factor. Lastly, an assessment is made on which areas are most vulnerable for erosion and sediment yield by combining the model outputs with existing literature on erosion models and geomorphic dynamics.

(11)

11

Figure 4a-d: Some of the flow characteristics needed for the calculation of the MUSLE: runoff curve number (a) and runoff volume (b) are determined by the rainfall event, land cover, soil hydrological group and antecedent moisture condition. The time of concentration (c) and peak flow rate are mainly a product of the type of flow distribution throughout the study area, the soil infiltration capacity and the rainfall event.

(12)

12

3. Results

3.1. Spatial Resolution

3.1.1. RUSLE

Figure 5a-d shows a clear visualization of the impact of spatial resolution on the RUSLE and MUSLE model outcomes. The most striking feature of this visualization are the differences in erosion and sediment yield patterns observed between different resolutions in both the RUSLE and MUSLE models. Distinct areas lower in the terrace system show to be most subjected to erosion in the low spatial resolution models while being almost completely ignored as erosion vulnerable areas in the high resolution models. This however is partially a matter of interpretation. While the highest erosion and also SY rates can indeed be found in this area, the large majority of this area has very low rates. Aggregating erosion and SY of the high resolution data on larger areas reveal a different pattern of areas vulnerable to geomorphic processes (Figure 6). Nevertheless, erosion vulnerable areas in the high-resolution model clearly coincide with topography; steep slopes result in higher LS-factor values, causing erosion rates to increase accordingly. This can be observed on the steep slopes in the upstream areas of the catchment which have higher erosion rates than the flatter areas downstream, but also on the terrace system. Here, topographical features such as gullies and rills cause abrupt and substantial changes in steepness enhancing the LS-factor values and determining erosion patterns (Figure 7a-c).

In the low-resolution erosion model, erosion patterns clearly differentiate in several ways. First, the erosion rates are more or less evenly dispersed throughout the study area. Second, topographical features such as gullies and rills are omitted from the model. Third, areas most prone to erosion do not coincide with those from the high-resolution model. However as with the high-resolution erosion map, it can be concluded that the erosion rates coincide with the LS-factor map (Figure 3 and Figure 5b). Moreover, the erosion rates calculated in the high-resolution maps are vastly

different than those in the low-resolution maps. Despite the vastly smaller range of steepness, slope length (Figure 8) and erosion rates in the low-resolution erosion model, the average erosion rate is higher than in the high-resolution map. This is a vast difference of a factor 3.36; 0.170 t ha-1 _yr-1_in the high-resolution map against 0.571 t ha-1_yr-1_{in the low-resolution map. Noticeable is that while} the maximum values in the high-resolution models are approximately 30 times higher than the low-resolution model outputs, the average erosion rates are significantly lower.

3.1.2. MUSLE

When analysing the outcomes of the MUSLE models at different resolutions, similarities can be observed with the RUSLE models. Again, topographical features such as gullies and rills that are causing the highest sediment yield rates in the high-resolution model are omitted in the

low-resolution model. However, several differences can be observed as well. First, the lower areas show higher sediment yield rates in both models than upstream areas. This is contradictory to the RUSLE maps where erosion rates are evenly dispersed or even somewhat concentrated in the upper areas. Second, both spatial resolutions of the MUSLE appear to show more resemblance with one another compared to the RUSLE models. Third, the rates in both sediment yield models are vastly different compared to the erosion models. Sediment yield as calculated in the high-resolution map depicts averagely 1.766 t ha-1_{across the study area while the low-resolution model shows an average} sediment yield of 0,060 t ha-1_{respectively. However, the crucial difference here is that in the MUSLE} models the low-resolution map is the model yielding vastly higher rates instead of the

(13)

13

Figure 5a-d: Erosion and sediment yield at low and high spatial resolutions. Erosion at high spatial resolution (a), erosion at low spatial resolution (b), sediment yield at high-resolution (c), sediment yield at low-resolution (d).

(14)

14

Figure 6: Aggregation maps of the erosion (left) and sediment yield (right) high-resolution calculations. Each cell equals 1 m2_{and therefore is the mean value of 100 cells as used in Figure 5}_a_{and Figure 5}_c_.

Figure 7a-c: The decisiveness of topography on creating erosion patterns. Slope (a), LS-factor (b) and erosion rate (c) for the terrace system. Red indicates high values, dark green indicates low values. Terraces are slightly turned eastwards for clearer projection.

(15)

15

Figure 8: The distribution of slopes on the high-resolution (left) and low-resolution DEMs. The y-axis represents the percentage of cells occupying each slope value. Each bar represents an integer slope angle in degrees.

3.2. RUSLE versus MUSLE

Comparing the RUSLE with the MUSLE model outcomes demands for careful interpretations. As the RUSLE computes annual erosion rates while the MUSLE calculates sediment yield after a rainfall event, both models depict different features at different time scales. Scaling both high-resolution maps relatively to their maximum output as shown in Figure 5ac (break values at 0.05%, 0.1%, 0.5%, 1%, 10% and 100%) to some extent enables making a comparison. This clearly reveals the larger relative differences in sediment yield in the MUSLE model compared to erosion rates in the RUSLE model; far more areas are classified in the middle classes in the RUSLE erosion model. The same is true for the low-resolution models. While the rates of erosion and sediment yield in the high-resolution models are different at some parts of the study site, the patterns are almost completely identical. This can clearly be observed when looking at zoomed-in displays of both models (Figure 9abef).

However, the aggregated maps of the high-resolution RUSLE and MUSLE models (Figure 6) reveal significant differences in erosion and SY throughout the study area. Here, it suddenly becomes clear that the areas that are most vulnerable for erosion and SY are distinctly different from one another. Erosion is most profound on the incised terraces in the upper area of the study area, while the SY maximum is located on the uppermost terrace of the terrace system. Moreover, a difference map as created in Figure 9g carefully displays the differences in both models relative to their maximum value. From the two gullies that can be observed, the right gully is red meaning that the RUSLE erosion rate relatively to its maximum here is higher than the relative SY rate of the MUSLE while the opposite is true for the left gully, therefore displayed blue. Figure 10 indicates peak flow rate and accumulation for the same area; peak flow is highest on the left gully while flow accumulation is highest on the right gully. While the patterns appear to be completely identical in this area, the relative rates show minor but significant differences. When looking at the entire terrace system in Figure 11, the predominantly blue colour again clearly indicates the larger relative differences in the MUSLE compared to the RUSLE model. Despite the differences in rates observed between both high-resolution models, the geographical patterns of these fluxes are still identical. This is not true for the low-resolution maps, where large differences can be observed both in patterns as in rates between different spatial resolution in both the RUSLE and MUSLE models.

Lastly, a comparison is made between all four models based on the impact of each individual factor on the outcomes of the model using the Pearson’s correlation coefficient (Table 1a-d). From this, several things can be concluded. First, both erosion maps are both only significantly correlated with

(16)

16 the LS factor showing the high influence of this factor on the model outcomes, most extensively on the high-resolution map. Second, both high-resolution maps are highly correlated with the LS factor map (r = 0,95 and 0,80) while both low-resolution maps clearly yield lower values. Third, both SY maps are less correlated with the LS factor map than their RUSLE counterparts. Lastly, the

correlation with runoff volume (Q) and peak flow rate (qp), but also K factor is significantly higher in the low-resolution SY model than in the high-resolution SY model. Although correlation must by no means be confused with causation, the beneath graphs can cautiously be interpreted as somewhat showing the relation between each component on the model outcomes.

(17)

17

Figure 9a-g:Patterns of erosion and sediment yield on the terrace system at different zoom-in scales. Erosion (a) and sediment yield (b) on the terrace system, sediment yield at closer scales (c,d) and sediment yield (e) and erosion (f) at a high close up. The difference between the RUSLE and MUSLE relative to their maximum (g) is calculated by subtracting sediment yield from erosion.

(18)

18

Figure 10: Peak flow rate (left) and flow accumulation (right) at the same area as Figure 9efg. On the left image, red indicates high values and green low values, while high values on the right image are displayed white.

Figure 11: Difference RUSLE and MUSLE for the entire terrace system. Again, the difference is calculated by subtracting MUSLE from the RUSLE in percentages relative to their maximum.

(19)

19

4. Discussion

The primary focus of this research is to assess the effect of spatial resolution upon identifying sheet and rill erosion vulnerable areas on a hillslope scale. This will be answered by thoroughly explaining each individual research subquestion, followed by a clear disquisition on the strengths and

weaknesses of the used method along with a summation of possible research improvements and future recommendations. Lastly, the interpretation and the degree of extrapolation of the research findings will be briefly discussed. The first two research subquestions are based on the behaviour of the models and therefore require a preponderant technical elucidation of the processes driving the model operation. The last research question requires a clear and brief overview of the similarities and differences in all model patterns and rates with the existing literature on geomorphic processes, from which a critical assessment is made on which areas are likely to be most vulnerable for erosion and sediment yield.

4.1. Spatial resolution

Spatial resolution defines the amount of detail within a model and therefore the incorporation of topographical features such as gullies, terraces and rills (Theobald, 1989). Low spatial resolutions fail to incorporate these features as a result of averaging: only one value is used for a relatively large area impeding the incorporation of any potential fluctuations within that area (Thompson, Bell & Butler, 2001). The transition of information over the same area therefore becomes less gradual creating a discretized model as can be seen in the LS factor maps in figure 3. Averaging can clearly be depicted when observing the model rates. Albeit the vastly smaller range of topographic features in the low-resolution erosion model (i.e. slope), total erosion rates are approximately three times higher in the high-resolution model. Schoorl, Sonneveld & Veldkamp (2000) state that the large surface area of each low-resolution cell generally results in an overestimation of erosion rates and is similar to many coarse resolution studies. They explain that lower spatial resolutions are always accompanied by a larger artificial mathematical overestimation, which can be corrected for by introducing a calibration factor.

This mathematical overestimation can also explain the large difference in average sediment yield observed between both resolutions in the MUSLE models. As the MUSLE model incorporates more spatial resolution-sensitive parameters than the RUSLE, more calibration steps are needed to obtain correct model outcomes (Beven, 1997). The replacement of the constant rainfall erosivity factor by two continuous factors runoff and peak flow rate increases the amount of spatial

resolution-sensitive factors and therefore enhances the error margin of rates between both resolution models. This again stresses the importance of introducing a calibration factor for obtaining accurate rates.

4.2. RUSLE and MUSLE comparison

The main difference between the RUSLE and the MUSLE model structures is the replacement of the rainfall erosivity factor by a runoff and peak flow rate factor. First, this explains the lower rates of sediment yield compared to erosion rates in the upstream areas of the study site. Shallow

concentrated flow and channel flow have the highest peak flow rates and can be found in areas with flow lengths upwards of 91,4 meters. As these flows generally occur in more downstream areas, sediment yield rates increase here accordingly (Figure 4d).

(20)

20 Second, the addition of two continuous variables is likely to cause the reduced correlation between sediment yield and the LS-factor in comparison to both RUSLE models (Table 1). Simultaneously, the reliance of the replaced factors is increased, especially on the peak flow rate. This can clearly be seen when looking at a detailed RUSLE-MUSLE difference map as can be seen in Figure 9g along with an accumulation and peak flow rate map of the same area (Figure 10). The flow accumulation map indicates higher values on the right gully which lead to a higher LS-factor. As the LS-factor is the only controlling factor on this small scale in the RUSLE model, the right gully yields higher erosion rates than the left. However in the MUSLE model, sediment yield is also determined by the peak flow rate; explaining higher runoff volumes on the left gully. This leads to overall higher relative sediment yield rates on the left gully in the MUSLE model compared to the RUSLE model, while showing relative lower rates in the right gully. As a result, the difference map displays the left gully blue and the right gully red.

An often used method to calculate SY is by calculating erosion and subsequently applying delivery ratios (Ferro & Minacapilli, 1995). Because erosion is partly dependent on the sediment being discharged with runoff and the sediment concentration of this flow, the inclusion of these variables in the MUSLE is known to generally yield more accurate results than applying delivery ratios (Kinnell, 2005). However, delivery ratios are not used in this research as they require an analysis of field data not measured in this study (Fernandez et al., 2003). Therefore, a more accurate comparison than executed in Figure 11 is prohibited.

4.3. Vulnerable areas for sheet and rill erosion

At first, this research shows that low-resolution model calculations are highly likely to display erosion and SY patterns less accurately than high-resolution models. Therefore, high-resolution models can best be used to allocate areas vulnerable for sheet and rill erosion in this research. The results reveal strong alignment with topography, with highest rates of erosion and sediment yield observed on steep areas, sudden topographic fluctuations such as gullies, incised terraces and flow accumulation pathways. However as previously stated, the allocation of vulnerable areas is an ambiguous process as it requires a consideration between intensity and degree of dispersion of both erosion and sediment yield. The aggregated erosion and SY maps of Figure 6 perhaps are the best indicators of vulnerable areas, because cell sizes of 10 centimetres as used in the high-resolution models are likely to be too precise for executing mitigation strategies. In these models, three areas can be allocated as most vulnerable areas in both models: the small incised terraces in the upstream area, the steep ridge on the middle left of the terrace and the incised uppermost terrace of the terrace system downstream. These areas are characterized as zones with high topographic fluctuations and low vegetation cover and in line with sheet and rill erosion vulnerable areas found in other research (de Vente, Poesen & Verstraeten, 2005). Subsequently, the research result of finding more profound sediment yields rates in downstream areas compared to erosion is also in line with comparable studies (Birkinshaw & Bathurst, 2006).

4.4. Research flaws, improvements and future recommendations

The main flaw of this research is the lack of actual data on the erosion and sediment yield rates of this area needed to verify the research results. However as this has not been the objective of this research, this might better be regarded as a recommendation for future research. While erosion is normally found to be difficult to measure manually, SY can relatively easy be calculated by

(21)

21 & Fletcher, 1999). These measurements will need to be corrected for as the study area represents only a part of the Puentes catchment.

Overall, it is important to acknowledge that the empirical approach of calculating erosion and sediment yield used in the RUSLE and MUSLE soil loss equations have the major benefit to be relatively easy applied in erosion studies while having the downside to simplify important feedback mechanisms decreasing the accuracy of the model outcomes (Yoder & Lown, 1995). Plant cover and land use are considered very important erosion parameters due to their influence on the infiltration rate, crust formation, organic matter content and aggregate stability (Garcia-Ruiz, 2010; Solé-Benet et al., 2010). However, in the RUSLE and MUSLE model structures these feedback mechanisms are not incorporated. This is likely to be the reason that plant cover and land use are clearly being inferior to topography in all model outcomes. Another likely reason for this is the vastly reduced accuracy from which these features are measured and mapped in comparison to topography. Lastly, Yoder and Lown (1995) report that it is important to acknowledge that some values used for the non-continuous factors (K-factor, C-factor) in both model equations in this research are based on statistical relationships mostly derived from studies performed in the United States. These values are likely to be different from European study sites due to differences in human intervention,

topography and climate (Hammad et al. (2004). Incorporating data from Klaver (2016) on the dispersion of erodibility throughout the terrace would significantly improve the research results, provided that proper units were used.

This research also shows the importance of introducing a calibration factor attuned on the number of spatial and temporal resolution-sensitive parameters. Another option would be the introduction of a fractal method: a method that compensates the roughness originating from the non-gradual transition between cells in coarser spatial resolutions (Naipal et al., 2015). These factors are not incorporated in this research as the emphasis lies on the detection of sediment yield and erosion vulnerable areas instead of rates, which would not change after introducing such a factor. However, it would significantly increase the research accuracy.

Another improvement of this research would be the adaption of the flow routing algorithm. Now, flow accumulation is calculated by a D8 single-flow algorithm by O’Callaghan & Mark (1984) as this is the only available flow-routing algorithm in ArcGIS 10.1 (Schäuble, Marinoni & Hinderer (2008). In this algorithm, all matter from the source cell is transferred to a single cell downslope whereas in reality matter from the source cell can be divided over multiple downstream “cells” (Desmet & Govers, 1996). A multiple flow algorithm is capable of simulating this and would be a model

improvement as it is found to give more realistic representations compared to single-flow algorithms (Butt & Maragos, 1998; Wilson et al., 2007).

4.5. Interpretation of research outcomes

When interpreting the results from this research, it must be realized that sheet and rill erosion are not the only processes that determine the geomorphic dynamics of this region. As Plata Bednar et al., (1997) conclude, it is most likely that gully and river channel erosion accounts for the major part of erosion in the Puentes catchment. Furthermore, the influence of spatial and temporal scale on hydrological modelling outcomes hinders the conclusiveness of this research. Small scale models have the major benefit to incorporate detailed information about the soil conditions and their spatial pattern distribution while the connectivity between morphological components from

different parts of the catchment is neglected (Lesschen, Schoorl & Cammeraat, 2009). The reverse is true for models on catchment scales. The omission of the remaining catchment area from this study

(22)

22 is an example of this. Finally, since erosion and SY rates and patterns of this research are linked to regional specific characteristics, the results must be interpreted accordingly. Still, the research findings can be used as an indicator of the geomorphic dynamics of this area on similar spatial and temporal scales and as a procedure of identifying areas vulnerable for sheet and rill erosion.

5. Conclusion

The objective of this research is to assess the influence of spatial resolution on identifying vulnerable areas for sheet and rill erosion. By focusing on erosion patterns rather than rates while using high-accuracy DGPS data, this research addresses the effect of spatial resolution on erosion modelling using two under researched methods. DGPS data is compared with remote sensing data (25-meter resolution) using the popular RUSLE and MUSLE erosion models in order to address the main research question: “What are the effects of using different spatial resolutions on identifying

vulnerable areas for sheet and rill erosion on a hillslope scale in southeast using the RUSLE and MUSLE erosion models?” As misjudging erosion prone areas can nullify the impact of erosion

mitigation strategies, a thorough investigation is evident.

In this research, spatial resolution greatly influences sheet and rill erosion patterns in both the RUSLE and the MUSLE erosion models. Both low-resolution models are unable of correctly displaying erosion and sediment yield patterns as they fail to incorporate distinctive topographical

characteristics. On the other hand, high-resolution models produce discrete patterns of erosion and sediment yield, which are highly correlated with the topography of the area. Furthermore, sheet and rill erosion patterns also vary significantly between the RUSLE and the MUSLE erosion models. Sediment yield is most abundant in the downstream areas as it relies on flow accumulation and routing patterns, whereas erosion is most profound in the upstream areas of the study site where highest topographic fluctuations can be found. Overall, the areas found to be most vulnerable to sheet and rill erosion in both erosion models are all characterized by high topographic fluctuations, either in the form of steep slopes or gully-incised terraces, together with low vegetation cover. A major flaw of this research is that the rates of all models are unusable for interpretation and cannot be compared to another due to the lack proper calibration. However as this does not affect the pattern display of erosion vulnerable areas, the significance of this research remains unchanged. Still, minor improvements could significantly improve the erosion modelling accuracy in all models. Introducing a calibration factor attuned on the number of spatial and temporal resolution-sensitive parameters and replacing the single-flow algorithm by a multiple-flow algorithm are a good example of this. Most importantly, it remains highly recommended that the results from this research are verified by actual data measurements on erosion and sediment yield in the future.

Overall, the reduced modelling accuracy and the allocation of erosion vulnerable areas as observed in this research correspond with the general scientific consensus on spatial and erosion modelling. The difference in erosion patterns between both spatial resolutions clearly displays the importance of incorporating small topographical fluctuations, but also stresses the large modelling inaccuracy associated with remote sensing data observed in this research. As increased spatial resolution also comes with increased costs of data collection and processing, a careful consideration between financial costs and modelling accuracy will determine the appropriate spatial resolution in future erosion modelling studies. However, this research reveals that a frequently used method like remote sensing with a 25-meter resolution is inadequate for the identifying of erosion vulnerable areas in an area that is heavily subjected to geomorphic processes.

(23)

23

References

Aksoy, H., & Kavvas, M. L. (2005). A review of hillslope and watershed scale erosion and sediment transport models. Catena, 64(2), 247-271.

Baartman, J. E., Veldkamp, A., Schoorl, J. M., Wallinga, J., & Cammeraat, L. H. (2011). Unravelling late Pleistocene and Holocene landscape dynamics: the upper Guadalentín basin, SE Spain. Geomorphology, 125(1), 172-185.

Bailey, S. W., Brousseau, P. A., McGuire, K. J., & Ross, D. S. (2014). Influence of landscape position and transient water table on soil development and carbon distribution in a steep, headwater catchment. Geoderma, 226, 279-289.

Beven, K. (1997). TOPMODEL: a critique. Hydrological processes, 11(9), 1069-1085.

Bhattarai, R., & Dutta, D. (2007). Estimation of soil erosion and sediment yield using GIS at catchment scale. Water Resources Management, 21(10), 1635-1647.

Birkinshaw, S. J., & Bathurst, J. C. (2006). Model study of the relationship between sediment yield and river basin area. Earth Surface Processes and Landforms, 31(6), 750-761.

Butt, M. A., & Maragos, P. (1998). Optimum design of chamfer distance transforms. IEEE

Transactions on Image Processing, 7(10), 1477-1484.

Calhoun, R. S., & Fletcher, C. H. (1999). Measured and predicted sediment yield from a subtropical, heavy rainfall, steep-sided river basin: Hanalei, Kauai, Hawaiian Islands. Geomorphology, 30(3), 213-226.

Cammeraat, E. L. (2004). Scale dependent thresholds in hydrological and erosion response of a semi-arid catchment in southeast Spain. Agriculture, ecosystems & environment, 104(2), 317-332.

Chaplot, V. (2005). Impact of DEM mesh size and soil map scale on SWAT runoff, sediment, and NO 3–N loads predictions. Journal of hydrology, 312(1), 207-222.

Chow, V. (1959). Open channel hydraulics.

Desmet, P. J. J., & Govers, G. (1996). A GIS procedure for automatically calculating the USLE LS factor on topographically complex landscape units. Journal of soil and water conservation, 51(5), 427-433. De Vente, J., Poesen, J., & Verstraeten, G. (2005). The application of semi-quantitative methods and reservoir sedimentation rates for the prediction of basin sediment yield in Spain. Journal of

Hydrology, 305(1), 63-86.

De Vente, J., Poesen, J., Verstraeten, G., Govers, G., Vanmaercke, M., Van Rompaey, A., ... & Boix-Fayos, C. (2013). Predicting soil erosion and sediment yield at regional scales: Where do we stand?. Earth-Science Reviews, 127, 16-29.

de Wit, A. M. W. (2001). Runoff controlling factors in various sized catchments in a semi-arid Mediterranean environment in Spain.

Farajalla, N. S., & Vieux, B. E. (1995). Capturing the essential spatial variability in distributed hydrological modelling: infiltration parameters. Hydrological processes, 9(1), 55-68.

Fernandez, C., Wu, J. Q., McCool, D. K., & Stöckle, C. O. (2003). Estimating water erosion and sediment yield with GIS, RUSLE, and SEDD. Journal of Soil and Water Conservation, 58(3), 128-136.

(24)

24 Ferro, V., & Minacapilli, M. (1995). Sediment delivery processes at basin scale. Hydrological Sciences

Journal, 40(6), 703-717.

García-Ruiz, J. M. (2010). The effects of land uses on soil erosion in Spain: a review. Catena, 81(1), 1-11.

Gillin, C. P., Bailey, S. W., McGuire, K. J., & Prisleyt, S. P. (2015). Evaluation of LiDAR-derived DEMs through terrain analysis and field comparison. Photogrammetric Engineering & Remote

Sensing, 81(5), 387-396.

Hammad, A. A., Lundekvam, H., & Børresen, T. (2004). Adaptation of RUSLE in the Eastern Part of the Mediterranean Region. Environmental Management, 34(6), 829-841.

Instituto Geográfico Nacional (2009) Modelo Digital del Terreno – MDT25. Centro de Descargas del

Centro Nacional de Información Geográfica. Derived from:

http://centrodedescargas.cnig.es/CentroDescargas/index.jsp

Kinnell, P. I. A. (2005). Why the universal soil loss equation and the revised version of it do not predict event erosion well. Hydrological Processes, 19(3), 851-854.

Lesschen, J. P. (2008). Multi-scale interactions between soil, vegetation and erosion in the context of agricultural land abandonment in a semi-arid environment.

Lesschen, J. P., Schoorl, J. M., & Cammeraat, L. H. (2009). Modelling runoff and erosion for a

semi-arid catchment using a multi-scale approach based on hydrological

connectivity. Geomorphology, 109(3), 174-183.

Lin, S., Jing, C., Coles, N. A., Chaplot, V., Moore, N. J., & Wu, J. (2013). Evaluating DEM source and resolution uncertainties in the Soil and Water Assessment Tool. Stochastic environmental research

and risk assessment,27(1), 209-221.

Moore, I. D., & Burch, G. J. (1986). Physical basis of the length-slope factor in the Universal Soil Loss Equation. Soil Science Society of America Journal,50(5), 1294-1298.

Moore, I. D., Grayson, R. B., & Ladson, A. R. (1991). Digital terrain modelling: a review of hydrological, geomorphological, and biological applications. Hydrological processes, 5(1), 3-30. Mulla, D. J. (2013). Twenty five years of remote sensing in precision agriculture: Key advances and remaining knowledge gaps. Biosystems Engineering, 114(4), 358-371.

Naipal, V., Reick, C. H., Pongratz, J., & Van Oost, K. (2015). Improving the global applicability of the RUSLE model-adjustment of the topographical and rainfall erosivity factors. Geoscientific Model

Development, 8, 2893-2913.

Natural Resources Conservation Services (NRCS), 1986. Urban hydrology for small 300 watersheds.

Technical Release 55, 210-VI-TR-55, 728.

O'Callaghan, J. F., & Mark, D. M. (1984). The extraction of drainage networks from digital elevation data. Computer vision, graphics, and image processing, 28(3), 323-344.

Panagos, P., Borrelli, P., Meusburger, K., Alewell, C., Lugato, E., & Montanarella, L. (2015). Estimating the soil erosion cover-management factor at the European scale. Land Use Policy, 48, 38-50.

(25)

25 Pérez-Sànchez, J., & Senent-Aparicio, J. (2016). Estimating Rainfall Erositivity in Semiarid Regions. Comparison of Expressions and Parameters Using Data from the Guadalentín Basin (SE Spain). Soil

and Water Research,11(2), 75-82.

Robinson, G. J. (1994). The accuracy of digital elevation models derived from digitised contour data. The Photogrammetric Record, 14(83), 805-814.

Savage, J. T. S., Bates, P., Freer, J., Neal, J., & Aronica, G. (2016). When does spatial resolution become spurious in probabilistic flood inundation predictions?. Hydrological Processes.

Schäuble, H., Marinoni, O., & Hinderer, M. (2008). A GIS-based method to calculate flow accumulation by considering dams and their specific operation time. Computers &

Geosciences, 34(6), 635-646.

Schoorl, J. M., Sonneveld, M. P. W., & Veldkamp, A. (2000). Three-dimensional landscape process modelling: the effect of DEM resolution. Earth Surface Processes and Landforms, 25(9), 1025-1034. Solé-Benet, A., Lázaro, R., Domingo, F., Cantón, Y., & Puigdefábregas, J. (2010). Why most agricultural terraces in steep slopes in semiarid SE Spain remain well preserved since their abandonment 50 years go?.

Theobald, D. M. (1989). Accuracy and bias issues in surface representation. Accuracy of spatial

databases, 99-106.

Thompson, J. A., Bell, J. C., & Butler, C. A. (2001). Digital elevation model resolution: effects on terrain attribute calculation and quantitative soil-landscape modelling. Geoderma, 100(1), 67-89. Vaze, J., Teng, J., & Spencer, G. (2010). Impact of DEM accuracy and resolution on topographic indices. Environmental Modelling & Software, 25(10), 1086-1098.

Wang, H., Wu, Z., & Hu, C. (2015). A comprehensive study of the effect of input data on hydrology and non-point source pollution modelling. Water Resources Management, 29(5), 1505-1521.

Williams, J. R., & Berndt, H. D. (1977). Sediment yield prediction based on watershed hydrology. Transactions of the ASAE, 20(6), 1100-1104.

Wilson, J. P., Lam, C. S., & Deng, Y. (2007). Comparison of the performance of flow‐routing algorithms used in GIS‐based hydrologic analysis. Hydrological processes, 21(8), 1026-1044.

Wischmeier, W. H. (1976). Use and misuse of the universal soil loss equation. Journal of soil and

water conservation.

Wu, S., Li, J., & Huang, G. H. (2008). A study on DEM-derived primary topographic attributes for hydrologic applications: sensitivity to elevation data resolution. Applied Geography, 28(3), 210-223. Yoder, D., & Lown, J. (1995). The future of RUSLE: inside the new revised universal soil loss equation. Journal of Soil and Water Conservation, 50(5), 484-489.

(26)

26

Appendices

Table 1: Pearson correlation coefficients

a.

b.

c.

Layer C factor K factor LS factor (low res) Q (low res) qp (low res)

Sediment yield -0.04169 0.27951 0.47972 0.32669 0.52955 d.

Layer C factor K factor LS factor (high res) Q (high res) qp (high res)

Sediment yield 0.04273 0.04610 0.80234 0.07960 0.13151

Table 1a-d: Pearson correlation coefficients indicating the correlation between each factor and the outcome of all four model: erosion at low-resolution (a), erosion at high-resolution (b), SY at low-resolution (c) and SY at high-resolution (d).

Layer C factor K factor LS factor (low res) Erosion -0.01214 0.06928 0.62102

Layer C factor K factor LS factor (high res) Erosion 0.04013 -0.03505 0.94745

(27)

27

Topography and hydrology

(28)

28

(29)

29

(30)

30

(31)

31

(32)

32

(33)

33

(34)

34

(35)

35

(36)

36

(37)

37

Appendix 1k: Data measuring points

(38)

38

RUSLE factors

(39)

39

(40)

40

(41)

41

(42)

42

(43)

43

(44)

44

MUSLE specific factors and additional data

(45)

45

(46)

46

(47)

47

(48)

48

(49)

49

(50)

50

(51)

51

(52)

52

(53)

53

The effect of spatial resolution on identifying vulnerable areas for sheet and rill erosion using the RUSLE and MUSLE erosion model in ArcGIS: a case study (Puentes catchment, southeast Spain)

Bachelor thesis