A QGIS -plugin for gully erosion modeling

Gully erosion affects the landscape and human life in many ways, including the destruction of agricultural land and infrastructures, altering the hydraulic potential of soils, as well as water availability. Due to climate change, more areas are expected to be affected by gully erosion in the future, threatening especially low-income agricultural regions. In the past decades, quantitative methods have been proposed to simulate and predict gully erosion at different scales. However, gully erosion is still underrepresented in modern GIS-based modeling and simulation approaches. Therefore, this study aims to develop a QGIS plugin using Python to assess gully erosion dynamics. We explain the preparation of the input data, the modeling procedure based on Sidorchuk’s (Sidorchuk A (1999) Dynamic and static models of gully erosion. CATENA 37:401–414.) gully simulation model, and perform a detailed sensitivity analysis of model parameters. The plugin uses topographical data, soil characteristics and discharge information as gully model input. The plugin was tested on a gully network in KwaThunzi, KwaZulu-Natal, South Africa. The results and sensitivity analyses confirm Sidorchuck’s earlier observations that the critical runoff velocity is a main controlling parameter in gully erosion evolution, alongside with the slope stability threshold and the soil erodibility coefficient. The implemented QGIS plugin simplifies the gully model setup, the input parameter preparation as well as the post-processing and visualization of modelling results. The results are provided in different data formats to be visualized with different 3D visualization software tools. This enables a comprehensive gully assessment and the derivation of respective coping and mitigation strategies.


Introduction
As noted by de Vente et al. (2013), Vanmaercke et al. (2021) and Borrelli et al. (2022), gully erosion is generally considered to be less studied compared to other soil erosion processes such as splash, sheet, and rill erosion.Nonetheless, gully erosion is responsible for substantial soil loss and sediment production on catchment scale (Valentin et al. 2005, Castillo and Gómez 2016, Märker and Sidorchuk 2003).In general, gully erosion is related to the formation and subsequent expansion of rills and channels generated by concentrated water flow (Poesen et al. 2003).According to Poesen et al. (2003), gullies are defined by a critical cross-sectional area of at least one square foot.A channel of this size can no longer be reconstituted by normal tillage operations.Based on morphology, gullies can further be categorized as ephemeral gullies and permanent gullies.Ephemeral gullies are characterized by a short-lived nature, shallow depth (below 0.5m), and smaller amount of soil loss in respect to permanent gullies.Ephemeral gullies are normally are filled in by farm equipment annually (Castillo and Gómez 2016).Although gully erosion has also been observed in areas without human intervention or influence (Abdulfatai 2014), anthropogenic activities are the main cause of the formation of gullies, for instance, deforestation, cultivating soil in fallow lands, and associated change of the hydrological conditions in the rainfall-runoff system (Sidorchuk 1999).
Gully development poses a threat to agricultural lands as well as communities and transportation infrastructure (Jahantigh and Pessarakli 2011).Land damage by gully erosion is generally a permanent non recoverable damage (Kuhn et al. 2023).It has both, on-site as well as off-site effects.On-site effects include soil loss, changes in land cover, economic loss, damage of natural resources, decreasing water level, damage to construction sites, partition of properties, decreasing land values, and increasing land accessibility issues (Arabameri et al. 2019).While off-site effects include the initiation and spreading of badlands, flooding, the production and transportation of substantial volumes of sediments along the watershed to lowlands, as well as sedimentation and contamination in dams, reservoirs, and water bodies (Dercon et al. 2012).Once gullies are formed, erosion continues to remove soil through the drainage channel and the gullies will continue to grow either in headward or sidewall direction unless measures are taken to stabilize them (Anderson et al. 2021).
Gullies can occur in a wide range of climatic zones, surface characteristics, and in different land use and land cover conditions.Their morphology also varies extensively (Bennett and Wells 2019).Gully erosion has been reported in different countries around the globe (see Nyssen et al. 2014;Vanmaercke et al. 2021).In terms of evolution, gully erosion has a greater impact in arid and semi-arid regions of the world, resulting in rapid depletion of soil moisture and ground water (Jahantigh and Pessarakli 2011).Similarly, pasture and cropland are more susceptible to gully formation than forests and urban areas.Moreover, areas with less compact or soluble substrates like sedimentary rocks are more prone to gully erosion than areas with igneous or metamorphic rocks (Castillo and Gómez 2016).
Gully morphology was modeled and monitored using various techniques.Nadal-Romero et al. (2015) used Terrestrial Laser Scanning (TLS) in combination with Structure-from-Motion (SfM) photogrammetric techniques to model erosion in high erosion-prone areas.Gong et al. (2019) used Unmanned Aerial Vehicles (UAVs) to generate orthophotos and a Digital Elevation Model (DEM) to analyze the morphology of the gully and identify fragile areas of erosion.Image-based 3D-modeling using a terrestrial photo camera and PhotoScan software from Agisoft, was applied by Frankl et al. (2015), resulting in a high spatial as well as temporal resolution model.Rodríguez et al. (2022) conducted a surface runoff test on laboratory scale, to study erosion phenomena using SfM, which can be applied for real world gullies.All these methods require fieldwork and are time consuming, especially if there are multiple gullies to be modelled.
On the other hand, Geographical Information Systems (GIS) play a crucial role in modeling location-based data.Das (2017) highlights the unique quality of GIS to simulate soil erosion processes using multiple data layers such as satellite imagery, digital elevation models, plus vegetation and soil characteristics.El Jazouli et al. (2017) uses GIS to implement the Universal Soil Loss Equation (USLE) and spectral index methods Normalized Difference Vegetation Index (NDVI), Coloration Index (CI), and Form Index (FI)) to assess the soil erosion susceptibility and to predict soil erosion rates.
The morphological characteristics of gully systems help us to better understand the development of gullies in order to identify the direction of its evolution and the areas that are likely to be affected.In addition, temporal modeling of gully morphology facilitates recognizing the current evolutionary stage of the gully's lifetime according to Kosov et al. (1978), i.e., the dynamic or static stage.Furthermore, information on gully evolution is useful to develop coping and mitigation strategies to control gully erosion.Actions may comprise e.g. the implementation of retaining vegetation along gullies, divert water from erosion prone areas, and to plant and retain deep rooted perennial pastures.
In this study, we present a plugin tool for QGIS, based on Sidorchuk's (1999) simulation model and earlier works of Omran et al. (2022) to assess the temporal changes in gully morphology at watershed scale.
(meters above sea level) on a north facing slope near the village KwaThunzi.The gully system drains into the Upper Mkhomazi River.The area is covered with moist grassland of the Sub-Escarpment Ecoregion and belongs to the geomorphic province of the Eastern Coastal Hinterland of KwaZulu-Natal (Partridge et al. 2010).Belonging to this warm temperate climate zone with dry winters and warm summers (Koeppen-Geiger classification Cwb) (Conradie 2012), it is characterized by a mean annual temperature of 14.9° C and annual precipitation of 788 mm.Typical for the South African Summer Rainfall Zone, precipitation is concentrated in the humid summer months while the winter months from May to August are dry.The gully is incised into Late Pleistocene colluvial sediments on the foot slope of a dolerite dyke.The sediments are ~ 10 m thick and underlain by shale stones of the Beaufort Group.

Environmental data
Various environmental data are required for the gully model, including: Digital Terrain Models (DTMs) of different stages of the gully development, start and end points of the gully drainage lines (point vector data), a range of soil layers, depths and properties of soil layers, as well as discharge data at the gully's head.

Digital Terrain Model:
The plugin requires two DTMs, i.e., one DTM without gully features or an early stage of gully erosion and one older stage or post-erosion DTM.
In our study an old stage DTM of the study area was derived from drone images using Structure-from-Motion algorithm (Omran et al. 2022).The pre-erosion DTM was reconstructed filling the eroded parts in the preerosion DTM to recreate the initial surface.Both DTMs have a spatial resolution of 1 m.

Start point and end point(s) of drainage lines:
The starting point represents the outflow location of the gully watershed, whereas the end points characterize the upstream parts of the gully system.Generally, gully head cuts can be used to identify the end points.The end points are extracted from the DTM and saved as a point layer.The algorithm allows for multiple end points so that multiple side gullies can be modeled that drain into the main gully.3. Soil Layer: Gully erosion may incise into the existing substrate and/or soil layers.Each layer is represented by the following parameters: depth, critical velocity, stable slope angle, and the erodibility coefficient.The critical velocity is used to initiate erosion at the bottom of the gully channel.Hence, if the flow velocity of the run-off is lower than the critical velocity, erosion does not occur.The stable slope value is used to model the sidewall instability, i.e., the change in the top width of the gully caused by the change of the gully geometry from rectangular to trapezoidal shape.The erodibility

Methodology
Two stages of gully development have been observed; a dynamic stage and a static stage (Sidorchuk 1999).The dynamic stage of the gully is short and generally covers only about 5% of its total lifetime (Fig. 2).However, the main gully erosion dynamics occur in this dynamic phase such as upstream headcutting, channel bed incision, and widening of channel cross-profiles (Campo-Bescós et al. 2013).Almost 90% of gully length, 80% of its depth, and 60% of the gully area is already developed in this initial phase (see Kosov et al. 1978).Because of the rapid and intense processes of this stage, mitigation and intervention measures should concentrate in this dynamic initial stage because it is more difficult and cost intensive to remediate already developed gullies.In contrast, in the static stage, the gully reaches its maximum volume and becomes more or less stable.This means the flow velocity is less than the threshold value for erosion initiation but is greater than critical velocity of wash load sedimentation.There are several well-defined parameters used to model gully erosion such as: i) the properties of the respective soils, ii) runoff or rainfall data, iii) catchment area, as well as iv) topography and morphology of the gully itself.These parameters are generally used in the simulation models to generate the morphology of the gully.Sidorchuk (1999) proposed a simulation model based on the equation for mass conservation (Eq. 1) and deformation (Eq.2) to model the temporal changes in the longitudinal morphology of the gully.
Here Q s = QC is sediment discharge (m 3 /s), Q = water discharge (m 3 /s); C = mean volumetric sediment concen- tration; X = longitudinal coordinate (m);C w = sediment concentration of the lateral input from the contributing catchment; t = time (s); q w = specific lateral discharge (m 2 /s);W = flow width (m); M 0 = detachment rate of the soil particles at the gully bottom (m/s);D = flow depth (m); M b = detachment rate of the soil from the channel banks (m/s); Z =gully bottom elevations (m); V f = fall velocity of sediment particles in turbulent flow (m/s), = soil porosity.
Both equations are used to calculate the gully's incision rate (for more detailed information see: Sidorchuk (1999Sidorchuk ( , 2021))).The gully's top and bottom width is represented by a trapezoidal shape.As stated by Sidorchuk (1999), during the initial stage of gully formation, flowing water erodes a rectangular channel in the top soil or gully bottom.Afterwards, the vertical sidewalls of the gully may become unstable and subject to mass movements and further erosion.Consequently, the rectangular cross-section transforms into a trapezoidal shape (Fig. 3).
To numerically solve the differential equations shown above we used the algorithm proposed by Sidorchuk (1999).To prepare the input data and visualize the results a QGIS plugin was developed and implemented with Python, allowing the assessment of the gully's geometry and evolution dynamics by wrapping the Sidorchuk model.The plugin is developed using the "Plugin Builder" in QGIS.This builder provides a Python framework to integrate the plugin into QGIS.
The plugin considers three groups of parameters which directly impact the gully erosion: i) topographical parameters including elevation points, stream network, and flow accumulation of gully sites, ii) soil parameters including critical velocity, soil erodibility, and slope stability, and iii) hydrological parameters including the discharge data at the start point of the gully.The simulation model of Sidorchuk (1999) was integrated into the plugin to model the gully morphology; additionally, three different algorithms were (1) Fig. 2 Gully morphology during its lifetime.1 = length; 2 = depth; 3 = area; 4 = volume (Sidorchuk 1999) integrated to predict the gully flow width, representing specific discharges for different environmental conditions.After an actual simulation run (including data preparation, integrating Sidorchuk's (1999) simulation model of gully erosion, transforming the local coordinates) a sensitivity analysis is performed.The first step regarding the gully model input preparation is to provide all the necessary topographical data such as DTM, start-end point vector file, and depths of soil layers.The second step concerns the integration of the soil characteristics and discharge data.Subsequently, we applied the model and the results were post-processed to be visualized in 2D and 3D (Fig. 4).

Data preparation
The data preparation procedure involves the preparation of input files for each gully channel to be simulated.The input file contains details of each gully point defined by a longitudinal step size along the gully's flowline.The main input is based on the DTM which characterizes the situation before gullying starts.A second DTM showing the evolved gully or a later stage in gully evolution that is subsequently used for the validation of the simulation.From the initial DTM, the location of the start and end points and the soil layer depths are extracted.For this we use the drainage network, derived with a flow accumulation algorithm (D8) from the preprocessed DTM.Our hypothesis is that the gully will potentially evolve along the existing drainage paths or drainage network that is derived from the initial topography.
In case of multiple end points, only the longest drainage segment, having the maximum distance to the gully outflow (start point), is used to represent the main stream of the future gully and all other points are taken as tributaries.In this case, the start point of the smaller gullies or tributaries are the points intersecting the major gully segment (Fig. 5).Subsequently, the drainage network based on a D8 flow accumulation algorithm is used to derive 3D world coordinates, i.e., Easting, Northing and height along the main drainage lines.Theses coordinates are used in the post-processing to embed the simulated gully into the initial DTM.In addition, the distance of gully points (along the gully channel) from the starting point are calculated, used for the numerical solution of the differential equations.The data preparation steps are shown in detail in Fig. 6.

Implementation of Sidorchuk's (1999) gully simulation model
The gully erosion simulation algorithm by Sidorchuk (1999) is integrated into the plugin using the data from the data preparation step as well as additional parameter files as input.The generated output contains the gully geometry for each gully channel (Fig. 7A, on the left) including: i) the total volume change for each year and for each point along the gully channel, ii) the initial surface elevation, iii) the gully depth after erosion, iv) the change in height, v) gully bottom width, vi) gully top width, and vii) the area of gully erosion (Fig. 7B).
The simulation model takes the flow accumulation one-by-one for each sub-gully, and processes it using the provided parameters as well as discharge values (Fig. 8A and B).The parameters include the critical velocity (m/s), stable slope (m/m), and erodibility coefficients (cm 3 /N-s).The critical velocity is used to estimate the minimum flow velocity required to initiate an erosion process (Shidlovskaya et al. 2016).According to the particle diameter some soil textures, for instance boulders, have a critical velocity of up to 4.0 m/s, while fine and medium sand have a critical velocity of 0.32 to 0.57 m/s.Major grain size classes and their critical velocities are shown in Table 1 (Sidorchuk 1999).The stable slope coefficient is used to identify the stability of the sidewall for a specific slope.However, slope stability also depends on the degree of consolidation (Gong et al. 2019).
The erodibility coefficient is the difference between the shear stress created by the flow and the soil critical shear stress (shear stress required for erosion initiation) and is used to estimate the amount of erosion (Mazurek 2010).Moreover, discharge data for the starting point representing the outlet point of the major gully segment, is required.This discharge value is split according to the flow accumulation values at the starting point of each gully side channel.
To model the flow width (W) as a function of the discharge (Q), different empirical formulae have been proposed.In the implemented plugin, the user may select between three different functions: (Here, W = Gully width (m); Q = Discharge value (m 3 /s)).
1. General Modeling (Sidorchuk 1999): Used as a default approach, as a generalized algorithm for the prediction of gully's width.
2. Revised Modeling (Sidorchuk (1999), adapted for permafrost areas): (3) W = 3.0Q 0.4 Fig. 5 Gully Modeling Structure with initial starting point, longest drainage line defined as main gully, as well as the distance from the starting point to the end points Fig. 6 Data Preparation including the spatial data input as well as the calculation procedures and output files Used to model areas where there is less erodible underground soil layer.This tends to result in gully systems that areas wider than deeper.
3. Cropland Area Modeling (Nachtergaele et al. 2002) Used to model cropland areas or particularly seedbed conditions, where soil is loose and smoothed.

Output postprocessing
The output of the simulation together with the coordinates generated in the data preparation step, are used to derive three GISlayers representing the gully's bottom along with the gully's left and right walls.Every layer contains the gully's 3D morphology for one year, which is used to visualize the temporal changes in the gully's morphology.In addition, a DTM raster is created for each year.Every DTM raster contains all gully channels of the same year, to visualize the overall erosion in the study area.

Validation, sensitivity and uncertainty analysis
We validate the model comparing the simulated morphology of the gully with the drone based DTM from 2018.Therefore, (4) W = 3.17Q 0.368   (5) W = 2.51Q 0.412   we delineate and compare the main morphometric characteristics of the gully.Moreover, we conduct a sensitivity and uncertainty analysis.A sensitivity analysis is a common practice used in various disciplines, e.g., social science, data science, machine learning, and environmental modeling (Saltelli et al. 2021).It is a technique to determine the reliability of experimental results by examining how the results are affected through changes in methods, models, and values of unmeasured parameters (Souza et al. 2016).Sensitivity analysis is on its way to becoming an essential part of mathematical modeling (Razavi et al. 2021).In this research, sensitivity analysis is used to determine how various parameters, used in the plugin, are affecting the generated morphology of the gully.The one-factor-at-a-time (OAT) approach is carried out to individually check the sensitivity of each parameter, while keeping other parameters fixed.The resulting morphology of the gully is analyzed to evaluate the most sensitive parameters of the gully erosion modeling application.

Gully modelling plugin
The implemented QGIS plugin can be used to simulate gully erosion for different time steps, e.g., daily, monthly or yearly time steps (Fig. 9).The plugin was developed and tested on For testing the plugin, a study area in KwaZulu-Natal, South Africa was used.The topographical parameters (i.e., DTMs, start and end points), the hydrological parameters (i.e., discharge data) and the soil characteristics for two soil layers were used to simulate the gully erosion dynamics.The soil layer shows a depth of 1 m, whereas the substrates reach 5 m depth.The characteristics of both soil layers are shown in Table 2.

Plugin output
The simulated gully morphology and its temporal changes can be visualized in two dimensions (2D) (Fig. 10A, on the left) and in three dimensions (3D).The DTM generated is shown in Fig. 10B.In order to visualize the DTMs, the raster DTM was converted to a Triangulated Irregular Network (TIN).Thus, the steep edges of the gully side walls were preserved.This conversion is part of the plugin.
The gully morphology evolution generated by the plugin (Fig. 11) is corresponding to the general gully evolution as postulated by Kosov et al. (1978), shown in Fig. 1.It can be observed that in the early stage of gully's lifetime the simulated erosion is higher; 78% of gully length, 53% of gully depth, plus 20% of both gully area and volume are formed during the initial 10% of its lifetime.However, the simulated morphological changes depend on the soil characteristics and the discharge data.Thus, it may vary with different parameters.Because of the change in depth and width of the gully at every point of the gully channel, the

Gully validation
For validation, the simulated gully morphology was compared to the morphology of the measured 2018 DTM.The hypothesis is, that according to Kosov (Fig. 2), gully length and depth develop quickly reaching close to final values, whereas the volume and area evolve much slower.Hence a comparison of length and depth allow a rough validation of the model, whereas the other two allow for an estimation of the gully evolutionary stage.Moreover, we compare crosssection profiles of the final DTM generated for a fifty-sixyear run and of the DTM obtained from the drone in 2018 (Fig. 12).It can be observed that the morphology of the generated gully matches with the existing gully at profiles 1 and 2 (Fig. 12).The statistical correlation has been calculated between the profiles at different locations in both DTMs: the modelled one (Plugin) and the DTM derived from dronebased photo-analysis.The result showed that there are high correlation values between profiles 1 and 2 at both DTMs ranging between 0.78 to 0.92 as seen in Fig. 12.Additionally, the length of the gully has been compared between both DTMs.The length values of the gully from the starting point to the end point reached 355 m for the 2018 drone based DTM, while the value for the modelled DTM (Plugin) reached 334 m.The area of the gully has been evaluated in both DTMs.Generally, the morphometric parameters such as cross section, length and depth correlate well between the modelled and real gully topography, as shown in Fig. 12.
Using the plugin, we performed a sensitivity analysis using different soil parameter values to evaluate their impact on the model results and thus to assess the overall uncertainty of the model results.The effect of on the different models for the relationship of gully width and discharge have been tested (Fig. 13A, B).In general, we observe that the erosion rates using the cropland model is higher than for the other two models (Fig. 13A, B).Both figures show that over time the changes of the top and bottom width using the cropland model is higher than for the revised model, which was already confirmed by Sidorchuk and Sidorchuk (1998).They stated that the rates of erosion will depend mainly on the temperature of the soil and the rate of runoff, both will increase with cropland soils rather than with permafrost soils.

Sensitivity and uncertainty analysis of gully plugin output
The sensitivity analysis is carried out to check the individual effect of soil parameters such as critical velocity, slope stability and erodibility coefficient on the morphological gully evolution.For this analysis only one soil layer is considered.The model is run for different settings, sixty-four model runs in total.In each run, one soil parameter is changed and the other parameters are kept constant.

Critical velocity parameter
The critical velocity is selected for different soil types (0.32 for Fine Sand, 0.70 for Soft Sandy Loam, 0.75 for Soft Loam, 1.00 for Hard Sandy Loam, 1.15 for Hard Loam, and 1.4 for Hard Clay: values in m/s, see Table 1).The generated gully bottom morphology during the fifty-seven years of gully formation is shown in Fig. 14.Generally, soils with higher critical velocities show a lower potential to erode, even in areas with higher flow accumulation (i.e., downstream areas).Instead, low critical flow velocities are coming along with notably high erosion potential, especially in areas with high flow accumulation (Fig. 14).The figure depicts the changes in gully volume over time for different soil types.Fine Sand and Soft Sandy Loam experience higher erosion compared to Soft Loam, whereas Hard Sandy Loam, Hard Loam, and Hard Clay are characterized by smaller changes.

Slope stability parameter
The slope stability parameter was chosen with different slope gradient values: 0.5, 0.7, 1.0 and 1.3 (m/m).We observed that the top width of the gully gradually changes with decreasing slope stability values, while there is not much difference in the bottom width of the gully (Fig. 15).Thus, there is minor change of the gully erosion volume.It can be concluded that the slope stability value is less sensitive and therefore does not significantly affect gully erosion.

Erodibility coefficient parameter
Similarly, we tested the effects of different erodibility coefficient values (0.001, 0.002, 0.003 and 0.004 cm 3 /N-s).It could be observed that the change of the erodibility coefficient value does not affect the change of gully erosion volume significantly, whereas the shape of the gully changes slightly for the different erodibility coefficient values (Fig. 16).The overall volume change over time of gully morphology is illustrated in Fig. 16.Thus, it can be stated that the gully erosion morphology is less sensitive regarding the erodibility coefficient.
The sensitivity analysis identifies as most sensitive parameter for gully erosion the critical velocity of the soils and substrates.Figure 17 shows the gully volume changes for different soil parameters.The volume of the gully increases with low critical velocity values, for example Fine Sand or Soft Sandy Loam soils.The result of the sensitivity analysis matches with the main results of Sidorchuk (1999), who stated that besides the discharge factor the erosion process is mainly controlled by the critical shear velocity of the soil.

Conclusion
Gully erosion is affecting more and more different areas in the world, particular those affected by climatic changes.The developed plugin is based on open-source Python libraries and can be used to model gully erosion dynamics through temporal changes in the gully morphology.Therefore, we use easily The plugin also provides options to select specific land use and climatic settings affecting gully formation (e.g., cropland or permafrost).We show with our case study that the network of existing gullies in a watershed can be modeled relative to the sub-gully end points provided as input to the plugin.The latter can be derived from any DTM with sufficient resolution.
The plugin was tested in KwaZulu-Natal South Africa.The area is highly susceptible to gully erosion due to highly erodible colluvial and alluvial valley fillings of the Masotcheni formation (see Bernini et al. 2021).The soils and/or substrates are characterized by different soil formation periods.Thus, gully erosion is destroying valuable fertile land suitable for agricultural production.The quality of the generated gully models depends on the accuracy of the input parameters: e.g., soil characteristics, discharge data, and DTMs.The resulting morphology can be analyzed to predict and estimate the amount of soil loss per year.Generally, we have shown that the model yields valuable results in terms of the morphometric gully parameters depth and length compared to the existing topography in 2018 that we used for validation.Moreover, we conducted a sensitivity analysis of the main input parameters of the plugin effecting gully morphology changes.Discharge data mainly controls the erosion and the volume increase of the gully followed by the critical velocity which is used to estimate the flow velocity required for erosion initiation.Slope stability controls the instability of gully walls and affects mainly the widening of the gully but is less sensitive.Moreover, our results show that also the erodibility coefficient controlling the overall erosion in a gully system seems to be less sensitive.The outputs generated using the plugin are 3D vector as well as raster data.This enables users to examine gully erosion and morphology within a GIS or other 3D visualization software.We have shown that the gully erosion plugin represents a valuable tool to prepare, run, and visualize the morphological dynamics of gully evolution throughout yearly timesteps.Due to the use of Open-Source software, our approach is especially suitable to assess gully erosion and related degradation processes in developing countries.Our plugin is a free alternative to costly commercial software and secondly, because our Graphical User Interface and streamlined workflow simplifies the access and opens the opportunities of gully modeling to a wider range of users.This is a step forward to more sustainable agricultural land use and agricultural landscape planning allowing scientists, planners, and administrative stakeholders to prevent, mitigate, or cope with negative impacts of gully erosion.

Fig. 7 AFig. 8 A
Fig. 7 A Simulation model integration (on the left).B Simulated output parameters with bottom and top widths.The number of profiles depends on the spatial resolution of the Digital Terrain Model (on the right)

Fig. 9
Fig. 9 Graphical User Interface (GUI) of the Gully Erosion Plugin with explanations of the fields to be filled in

Fig. 12
Fig. 12 Cross Section Profiles for the generated DEMs

Fig. 13 A
Fig. 13 A Gully top width evolution over 57 years for different environmental settings (on the left).B Gully bottom width evolution over 57 years for different environmental settings (on the right)

Table 2
Soil characteristics at gully bottom