SIMWE model application on susceptibility analysis to bank gully erosion in Alto Douro Wine Region agricultural terraces

Article history: Received 29 September 2016 Received in revised form 17 January 2017 Accepted 29 January 2017 Available online 3 February 2017 This paper focuses on the susceptibility evaluation to bank gullies on earthen embankments through the application of SIMWE (SIMulatedWater Erosion)model, using a high resolution digital elevationmodel (1meter spatial resolution). The results provided by themodel are comparedwith the hydrologic characteristics, soil texture and soil structure of the agricultural terraces. This approach demonstrates an association between the spatial distribution of erosive forms with high values of water depth and reduced water discharge that are consistent with the lower values of electrical resistivity. The areas with the highest percentage of erosive forms, related to sediment flux, transport capacity and sediment concentration susceptibility, assume medium values. These figures, combined with a low hydraulic conductivity and soil infiltration capacity, are consistent with the fine texture of soils, allowing increased runoff and the development of linear erosion. © 2017 Elsevier B.V. All rights reserved.


Introduction
The wine production in Alto Douro Wine Region -one of the world's oldest regulated and demarcated wine region -is based on a slope system organized in agricultural terraces once supported exclusively by dry stone walls, which has been undergoing the necessary changes for the introduction of technological innovations partially associated to the mechanization of vineyards work. In this sense, different forms of terrain framing have been implemented, namely the substitution of stone walls by earth embankments. This evolution raises a group of problems related to the hydric soil erosion and landscape preservation, since Alto Douro Wine Region is classified as UNESCO World Heritage Site since 2001.
Several works on hydrogeomorphic issues have been developed with the application of data acquisition techniques using detailed information of terraced areas.
To study the relationship between the hydrological processes and the stability of terraced areas, Camera et al. (2015) combined a simplified unsaturated-saturated hydrogeological model (STARWARS), which allow the determination of the changes in soil moisture and the groundwater distribution, being possibly to lead to the formation of local perched groundwater tables (PGTs), with a stability model (SARGLE). The authors consider the water table peaks are the hydrogeological occurrences that most affect the instability. So they tested them with daily and six-hour time steps. They admit that after 6 h part of the precipitation does not infiltrate leading to surface runoff, and they consider that a greater amount of 24 h average precipitation feeds the ground water table. The conclusions state that this hydrological model is a useful tool to determining the development of the water table. López Vicente et al. (2013) used the index of connectivity and a modified version of the revised Morgan, Morgan and Finney (RMMF) model to predict the hydrological connectivity and the rates of soil erosion. The index applied considered the characteristics of the drainage area and the flow path length used by a particle at the nearest sink. The authors conclude that the index represents the intense hydrological connectivity processes that occur in the irrigation channels and the walls of the agricultural terraces. Hussein et al. (2016) presents a study of terracing implementation based on RUSLE principles in a rain fed region. The results shows that terrace implantation is feasible in areas with a high slope degree and a seasonal precipitation of over 600 mm. Liu et al. (2013) analysed different scenarios, with terrain protection measures and different spatial terrace configurations simulated with the Water Erosion Prediction Project (WEPP) model, to determine soil erosion. The results showed that the soil loss increase from the top to bottom of slope with the higher precipitation intensity, generally in slopes without protective measures and no land consolidation. In different spatial configuration scenarios of high, medium and low grass coverage showed that protective measures are more efficient when located at the lowest position along slopes.
Advanced techniques have been used in terraced vineyards topography modelling, namely the Terrestrial Laser Scanner (TLS) and the Laser Scanner (LiDAR). Preti et al. (2013) used a centimetre resolution topography obtained from TLS survey in the analysis of terrace failure processes in vineyards. Tarolli et al. (2015) used Lidar elevation data to establish hydro-geomorphological analyses. The authors studied the erosion caused by agricultural terraces and roads through the application of the Relative Path Impact Index (RPII) as well the defined thresholds that allowed them to create scenarios. From this study they concluded that high resolution topography is very effective in the study of superficial erosion in terraced vineyards. Sofia et al. (2014) presents a study that demonstrates the potential of digital terrain models (DTM's) obtained from LIDAR data to determine spatial heterogeneity in terraced landscapes. Thus, the authors apply an algorithm, called Slope Local Length of Auto-Correlation (SLLAC), in two types of landscape with real and simulated DTM´s. (Sofia et al., 2016) The authors perform an analysis of automatically terracing mapping from DEM's obtained by Pleiades images and compare the results with DEM's obtained from LIDAR data. For the application of this methodology they applied two algorithms, the fast line segment detector (LSD) algorithm and the geomorphometric method based on surface curvature. From the conclusions they affirm that effective Pleiades DEMs are an alternative to the LIDAR DEM's in the monitoring of large agricultural areas.
The erosive forms of this area are not classical gullies as traditionally described in geomorphology studies, but are illustrative of erosional dissections referred to as 'bank gullies' (e.g. Poesen et al., 1996;Vandekerckhove et al., 2000). They usually are ephemeral gullies limited in length, depth and width, characteristics that are influenced by the morphology of agricultural terraces that do not exceeds 5.9 m of width and 63°of riser inclination (Fig. 1A).
Bank gullies are formed due to a retreat of headboard of a terrace in erodible hillslopes, upon the occurrence of a critical superficial flow that frequently coincide or follows high intensity rainfalls (Poesen et al., 1996;Meyer and Martínez-Casanovas, 1999;Vandekerckhove et al., 2000).
The changes of slope profile, after terracing ( Fig. 1B), builds a platform that can be divided in two sections: one excavated (point 2) and another filled (point 1). The filled section represents topsoil materials that have been (re)moved from above, representing friable materials derived from metasedimentary outcrop, and based at the board of the platform. This non-resistant material constitutes the point of abrupt increase of the slope angle that makes the transition between the platform and the earth embankment. This part of the terrace is under the influence of a high subsurface flow hydraulic gradient, combined with very week infiltration capacity, which provides surface runoff and promotes the linear erosion (Lesschen et al., 2008).
The bank gullies are located on the exposed area of the slopes, where the materials are desegregated and non-cohesive ( Fig. 1, Hte). They have small dimensions, with a length under 3 m and maximum depth about 0.40 m. For this reason, we have chosen to illustrate the inventory through a point at the head of each bank gully, on the terrace platform. The most acceptable geomorphological representation is by polygons or lines but, in this case, considering the small dimension and the great number of forms, becomes impracticable and don't bring additional information to modelling validation process. With a digital elevation model with 1 m of pixel, the representation projected in plan will always correspond to one point due to the morphology of agricultural terraces (i.e. Vandekerckhove et al., 2000;Kukemilks and Saks, 2013).
The soil erosion by surface runoff is a process difficult to model, given that it requires a detailed parameterization of all the concerning factors, which in turn depends on the scale analysis and the type of erosion forms to be modelled (Merrit et al., 2003). Revealing the extensive work on this issue, it is possible to identify more than 80 models that incorporate erosive processes (Lilly et al., 2009;Karydas et al., 2012).
The application of physically based models for erosion susceptibility evaluation in agricultural terraces, is fairly limited (e.g. Catalão and Pacheco, 2010;Figueiredo et al., 2013;Pacheco et al., 2014;Fernandes et al., 2014). In fact, "Gully erosion modelling has focused more on development of qualitative and empiric-statistical models than in the formulation of physically based models" (Meyer and Martínez-Casanovas, 1999, p. 320).
SIMWE, a physically based model designed by Mitas and Mitasova (1998), was developed to modelling linear erosion processes. Integrated in the software Geographic Resources Analysis Support System -Geographic Information System (Grass Gis), was applied by several researchers (e.g. Mitasova and Mitas, 1999;Koco, 2011;De Rosa, 2004;Fárek et al., 2014), allowing to obtain an association between the terrain's characteristics and the pattern balance between erosiondeposition.
The main purpose of this study is the evaluation of the susceptibility to linear erosion using the SIMWE model in a micro-watershed of 'Quinta de S. Luiz', an important wine estate of Demarcated Douro Region. Subsequently the results are compared with hydrological parameters obtained from field monitoring and are submitted to statistical validation with the field inventory, in order to verify if this model may be applied in a modified terraced area and a specific (vineyards) agricultural land use.
In this inventory we chose to maintain the term "bank gully" for forms that show a depth equal or above 30 cm (Morgan, 2005a,b). The metamorphic formations of Bateira, Ervedosa do Douro and Rio Pinhão forms the dominant lithology of the study area, mainly corresponding to phyllites and metagreywackes (Sousa and Sequeira, 1989) (Fig. 2C).
Originally the soils of this area are defined as leptosols associated to luvisols (Ribeiro, 2000), but according the World reference base for soil resources (WRB, 2014) and due to the terraces construction process and its long and intensive agricultural use, they are undoubtedly anthrosols. In this context, the observed texture reveals a high percentage of fine elements, variable between 32.1% and 64.9% of clay and silt. The average annual precipitation presents a value of 686.8 mm, according the 1983-2012 data series of the meteorological station of Adorigo.  However, the DDR is sporadically affected by intense but shortlasting precipitation episodes which are responsible for the initiation and development of erosion processes, as happened, for example, in the winter of 2000/2001 or in the beginning of the month of July in 2014, in Pinhão area. In fact, the climate of this area is "accentuated by its geographical position", revealing extreme values of average annual temperatures ranging from 2.2°C minimum to 29.6°C maximum (according to Sea and Atmosphere Portuguese Institute -IPMA), associated with reduced annual precipitation levels but presenting sporadic episodes of high intensity rainfall, of about 98.4 mm/day (Daveau, 1977(Daveau, , 1995. The presence of the Alvão-Marão and Montemuro mountains, to the west, as well the strong incision of the river system, are relevant geomorphological elements that form a barrier to the Atlantic wet air circulation, giving rise to a Mediterranean climate of continental influence.

Materials and methods
The SIMWE, developed by Mitas and Mitasova (1998, p. 505), is a physically based model that allows to simulate "(…) erosion, sediment transport, and deposition by overland flow, designed for complex terrain, soil, and cover conditions". Following the fundamental theoretical principles of the Water Erosion Prediction Project (WEPP) designed by Flanagan and Nearing (1995), its application uses the software GRASS GIS. It includes the calculation of two modules: the 'r.sim.water', that simulates the superficial water flow, and the 'r.sim.sediment' which calculates the sediment transport by the superficial drainage, according to the following equations (Julien et al., 1995): r (x, y)coordinates derived from DEM (Hann et al., 1994).  (DEM), was created from aerial photographs, acquired by a digital photogrammetric camera with a flight about 5000 m high. A stereopair of photos, with a pixel size of 50 cm, and longitudinal overlap of 60%, was oriented using 6 ground control points surveyed with a dual frequency GPS receiver, in differential, kinematic mode. Then, a stereomatching procedure was carried in order to obtain a dense point cloud, and finally, a regular grid DEM, with a grid spacing of 1 m. This DEM is in fact a digital surface model, since it represents elevation of any objects that exist on the ground, such as vegetation and buildings. In order to assess its accuracy, a set of points, also surveyed by GPS, on bare soil, in the Quinta de S. Luiz vineyards, were used as vertical check points. A root mean square error of 0.80 m in height was found, for a total of 121 check points. With this flight high and the variations of study area topography the accuracy of the model is reasonable. However to describe more accurately the erosive processes it is estimated that it would be necessary obtain a lower flight with a resolution that would allow a lower mean square error. All the photogrammetric processing was carried in Agisoft PhotoScan software.
This DEM works as input data for the calculation of the previously reported two modules of SIMWE, "r.sim.water" and "r.sim.sediment", as well as the "r.slope.aspect". This one uses the DEM to define the flow gradient vectors which will constitute the input data of the subsequent modules. The "r.sim.water" allows the simulation of the superficial water flow, and the "r.sim.sediment" module simulates the transportation of sediments by superficial drainage, relating the variation of sediment storage with the runoff along the slopes (Mitas and Mitasova, 1998).
Sequentially were estimated additional parameters that control the simulation model (Table 1). The rainfall was obtained from meteorological station of Adorigo, considering the 1983 to 2012 series, based on an event (6/2/2001) with 51 mm of daily precipitation on which was subtracted the average of the stabilized infiltration rate, (0.619), shown in Fig. 6. The other parameters were calculated through field survey and in situ tests.
Methodologically we measured these parameters at the top, middle and bottom of the watershed. The small local variations of the hydrological and textural of the materials along the watershed can be better analysed using this spatial monitoring.
Concerning the hydrological parameters, the saturated hydraulic conductivity (kfs) was measured through Guelph permeameter (a constant charge device), following the linear conductivity method and applying two water columns or pressure charges, at 5 and 10 cm (Reynolds and Elrick, 1986) (Fig. 4). The calculation of kfs was estimated from the following equation: where: Kfs is the saturated hydraulic conductivity expressed in cm/sec; Y is the value used for the constant associated to the linear method and it is expressed in cm 2 (2.16 cm 2 ); R1 is the value of the constant with stabilized water -5 cm (H1/ 60 s); R2 is the value of the constant with stabilized water -10 cm (H2/ 60 s).
The low permeability of soil, obtained in all tests, is consistent with its high percentage of clay and silt (Fig. 4).
To measure the electrical resistivity was used a resistivimeter (ABEM SAS 300C), disposing electrodes according Wenner array (Loke, 1999). Fig. 7. Results of SIMWE model. The electrical resistivity profiles are shown in the Fig. 5. Were also measured in different slope sections (top, middle and base), near the water lines, and the results derived from the following equation: where: ρ is the electrical resistivity; a is the distance between electrodes; R is the value registered by the equipment. According to Solberg et al. (2011), resistivity values that range 0-10 Ω * m correspond to fresh pore-water and stable clay structure. The values that range 10-100 Ω * m are related with fine-grained (silt) and leached materials. Values greater than 100 Ω * m area associated to dry crust clay deposits and coarse sediments. The infiltration capacity was evaluated using a double ring infiltrometer. The readings were realized with a 1-minute interval, until total infiltration was verified in the internal ring or until the stabilization of the level (Fig. 6), and the subsequent calculation follows the Libardi (1995) where: i represents the infiltration rate; l is the accumulated infiltration and t is the time.
In the upper area of the basin are registered fewer bank gullies due to the textural characteristics of the materials, since they present a large proportion of clay and silt. However, sub-surface water circulation process is fast due to the presence of voids in the soil, which can be suggested by measurements of hydraulic conductivity (cf. Fig. 4). Furthermore, resistivity profiles suggest a strong presence of water in the soil that promotes the occurrence of landslides rather than linear erosion caused by runoff, which was observed on the terrain. Instead, at the bottom of the basin was registered a large number of erosive forms, due to the fact that in this area infiltration capacity and hydraulic conductivity present lower values.
To validate the results of SIMWE model and terrain monitoring, we used the success-rate curve. The "success-rate curve", is based on the comparison between the prediction modelling and the erosion forms inventoried (Chung and Fabbri, 2003). Area Under the Curve -(AUC) evaluates the global prediction capacity of the model, ranging from 0 to 1 (Beguería, 2006).

Results
The results of modelling within the assigned conditions are presented in the maps of Fig. 7. These maps were represented into four classes (demonstrating reduced, medium, high and very high susceptibility), following an empirical representation focused in the lowest obtained results.
Generally, most of the bank gullies (91.2%) are concentrated in areas where the water depth assumes medium and high values (0.3-4.8 cm), and only 8.8% in very high class (N 4.8 cm), (Fig. 7a and Table 2).
The perspective is slightly different when it comes to the water discharge ( Fig. 7b and Table 2), with 100 % of the erosion forms located in the reduced to medium classes values (10-40 cm 3 /s). The category where the water discharge presents reduced values, there are 66% bank gullies.
Referring to sediment transport capacity ( Fig. 7c and Table 2), almost all erosive forms (98.5%) concentrates where materials mobilization reaches medium or high values (250.3 g/ms-9.4 kg/ms). In this sense the material texture, which presents a percentage of fines variable between 52.5-75.6%, can determine a larger particle detachment (cf. Fig. 3 and Table 2).
Regarding the sediment flux ( Fig. 7e and Table 2), the bank gullies (88.2%) once again overlap the medium or high classes (13.8-334.2 g/ms), contrasting with a percentage of 11.8% and 0% in the reduced and very high classes, respectively.
On the other hand, the classes of medium and high sediment concentration, (Fig. 6d and Table 2), register 70.6% of bank gullies. However, the very high class also registers a considerable number of erosive forms (28%).
Evaluating the balance between erosion-deposition (Fig. 7f), it is clear the concentration of forms in the sectors identified as susceptible to erosion (0.5-3.5 kg/m 2 s), registering 61.8% of gullies. However, there are still a considerable percentage of erosive forms identified by SIMWE in deposition areas (38.2%).
Therefore, analysing the erosion/deposition map, is possible conclude that the area under the curve (AUC) has a value of 0.619 of the susceptible areas to linear erosion. (Fig. 8).
Analysing the resistivity data, we see that Profiles 3, 5, 6, 7 and 8 have the higher values at the surface. For this reason they relate to the presence of superficially disaggregated material that results of the construction of the risers. The only difference between these four profiles is that the Profiles 7 and 8 are closer to the water line.
Profiles 1, 2 and 3, located at the lower end and medium area of watershed, have lower values of electrical resistivity at 2 m of depth, which indicates the presence of water. Nevertheless, there are no many erosive forms in this area. In follow-up the Profile 7, located in the upper section of the watershed, registered the lowest values of electrical resistivity clearly denouncing the presence of subsurface water. These values combined with the high infiltration capacity promote the occurrence of landslides, (present nearby terraces but not analysed in this text). That explains no superficial runoff and justifies the absence of erosive forms in this area.
Relatively to hydraulic conductivity, it is important refer that experiment only includes the 0.50 m of soil depth. So, it is clear that the points SL1 and SL3 present higher values of water circulation in the soil, although the experiments have been performed on a concave area. The points SL2, SL4 and SL7 illustrating low values in hydraulic conductivity. The SL5 and SL6 points located in the upper section of the watershed show medium values of hydraulic conductivity, considering that this area has the presence of water at 2 m of depth shown in the electrical resistivity profiles (Fig. 5). Concerning the infiltration capacity, the experiments reveal homogeny results, denoting a low infiltration capacity with values from 0.2 to 1.4 cm/min. It's interesting to note that the experiences that denoting higher values of infiltration rate are located in the upper sector of the watershed (points 5, 6 and 7).

Discussion
The SIMWE model has a medium adjustment to study area and presents a success rate AUC of 62%. This performance is due to the existence of different type of ground frame system. On the vertical vineyards there is no particles disaggregation derived from the terrace construction and the flow velocity is slower due to the vegetation. Due the construction process of agricultural terraces, the surface coarse material is desegregated and thus the infiltration flow is quick. Although, subsuperficial the hydraulic conductivity is very low, conditioned by the strong presence of clay and silt in the soil. Thus, it is possible to attend these results and compare them with the resistivity profiles. The presence of voids near the surface, from slightly cohesive material is consistent with high values of electrical resistivity, and a saturation on the underneath layer due to the preponderance of fine materials in the soil texture, explains the very low values of electrical resistivity.
Nevertheless, there are a great number of erosive forms in the lower area of the watershed. This clearly suggests the influence of the topography on the hydro-geomorphological dynamic modelling, namely where the upstream contributing areas to runoff are more significant (e.g. Prosdocimi et al., 2016). Therefore, the SIMWE model is sensitive to the influence of the morphology of agricultural terraces represented on the DEM. So, the erosive forms inventoried are related to medium and high classes of the water depth. The very high classes of water depth are not representative of soil erosion by gullies because are located in the stream line, with the highest contributing areas, is an evidence common to all topographic areas. Moreover, according to Morgan (2005a,b), the erosion is lower with the water depth increase, because the energy of rainfall is dissipated in the water and does not make changes in the surface soil, avoiding the rainfall splash impact on the soil surface. The medium and high classes of the sediment transport capacity and sediment flux parameters are coincident with most of erosive forms recorded. On the other hand, there is an elevated registry of erosive forms in high and very high classes of sediment concentration. The existence of disaggregate material at the surface and the morphology of agricultural terraces increase the resistance to water flow. On the other hand, with the decrease of the superficial flow velocity, the sediment transport capacity is decreased, like this it becomes evident the sediment concentration (Li and Abrahams, 1999). The existence of a large number of erosive forms in the higher classes of sediment concentration relates to the development of erosion in the sediment deposited during the previous sedimentation process that explains that the transport limited by erosion/deposition balance is favourable to the erosion process. The material is disaggregated and therefore the detachment of particles becomes more prominent. For this reason, the results of the transport limited by erosion/deposition show a considerable erosive forms percentage in deposition area. This result is due also to the fact that the erosive forms are very small and there isn't a pronounced distance between erosion and deposition areas. Relatively to water discharge, spatial variation seems directly associated to the general inclination of the slope, but there are factors that decrease the velocity. The main factors are related to the existence of vegetation or gravel on the earth embankment and on the platform. The vegetation slows the water flow, and gravel promotes the infiltration to the beneath soil layers through the macro pores, as evidenced in consistent studies e.g. Jomaa et al. (2012), Cerdà (2001), Figueiredo and Poesen (1998), Poesen and Lavee (1994). On the other hand, the morphology of the terrain, with the agricultural terraces implementation, acts as an attenuate element to the superficial flow that does not processed in a "free way".

Conclusion
The obtained results seems to indicate a moderate performance of the SIMWE model in the identification of the areas more susceptible to the bank gully erosion of the soils in terraced area with earthen embankments.
The erosion of this area is related mainly with linear flow processes, such bank gullies that reach a non-significant depth. On the other hand, there was a large number of bank gullies in the areas classified by the model with deposition processes. This is justified by the intrinsic morphology of the terrain. The deposition areas always coincide with the terrace platforms and the bank gullies affects the riser of the terrace but seldom they have evolved to the platform and therefor affected the deposition area.
There was a greater number of bank gullies in the downstream area of the watershed, where the agricultural terraces are less recent. Thus in the oldest terraces the soil materials are more compact and promotes runoff. In the most recent agricultural terraces the material is more disaggregated and enables the development of instability associated with subsurface flow, particularly landslides.
In the areas without continuous vegetation the flow proceeds superficially influenced by the weak infiltration capacity and hydraulic conductivity.
The validation process may be limited by the incompleteness of the inventory. The group of registered occurrences are represented by defect, since the erosion forms are frequently erased by the agricultural activity.
Since the AUC of the Success Rate Curve is about 0,62 we can state that the SIMWE model has a moderate to high adjustment to the bank gully erosion on this area.
That is because the area have a morphology conditioned by the length and inclination of the general slope and therefore by the contributing areas. Thus, the erosive forms on agricultural terraces are fully committed by these and other factors that come from anthropic intervention of this area, namely the surface disaggregation of material during the construction process and the presence or absence of vegetation.
For these reasons the SIMWE model is not robust enough when applied to agricultural terraces, this fact could be improved with a higher accuracy of digital terrain model.