A GIS-based simulation and visualization tool for the assessment of gully erosion processes

ABSTRACT

Gully erosion is known to affect soil productivity, and to limit land use in many parts of the world. To assess the susceptibility of soils to gully erosion, the triggering factors must be analyzed. A simulation of gully development over time was implemented in Python – GIS modules based on the topographic characteristics of the study area, soil properties, and measured runoff. The results show that the developed modules allow a successful data preparation, process simulation and spatial visualization of the obtained results in 2D and 3D as well as an analysis of time series of gully evolution .

KEYWORDS:

• Soil erosion
• Gully simulation model
• Python application
• GIS
• visualisation

1 Introduction

Soil erosion is recognised as a global environmental and agricultural issue threatening especially semi-arid areas such as in Australia, Russia, South Africa, Chile, Iran and the Mediterranean Regions of Europe. It has been estimated that about 75 billion tons of soil per year is eroded by wind and water-related processes often impacting agricultural land (Borelli et al., 2017). Around 70% of South Africa are affected by soil degradation (Hoffman et al. 1999, Garland et al. 2000) and 50% of the country’s area has a moderate-to-severe erosion potential exceeding 12 t/ha/a (Le Roux et al. 2008). Due to mountainous topography, erodible soils and cyclonic precipitation patterns, the Provinces of KwaZulu-Natal, Eastern Cape and Mpumalanga are most susceptible to these processes (Le Roux et al. 2008). According to the national soil erosion hazard assessment (Msadala et al. 2012), our study area lies in a region of moderate-to-high erosion susceptibility.

Soil erosion by water, such as gullying, is acknowledged as a global main driver of land degradation (Park et al. 2009). Gully erosion is typically defined as a deep channel that has been eroded by concentrated water flow, removing surface soils and materials (Torri et al. 2018). Gullies are distinguished from rills based on a critical cross-sectional area that corresponds to the size of the channel that can no longer be erased by normal tillage operations (Baade 1994, Bull and Kirkby 1997). Gully erosion is related to a variety of on-site and off-site damages, particularly impacting agricultural areas. On-site damages include the direct loss of arable land resulting in significant soil losses, reduced soil fertility, loss of crops and vegetation cover, as well as damages to infrastructures such as roads, power lines and communication networks (Ionita 2006, Haregeweyn et al. 2008a, Hayas et al. 2017). As an outcome of their high erosion rates, gullies enhance the connectivity of slope systems and lead to an increased sediment delivery. Estimates indicate that gullies are accountable for up to 80% of a catchment’s mean sediment yield (e.g. Flügel et al. 2003, Maerker 2001; Vanmaercke et al. 2012). In particular, the produced sediments washed into the drainage network are related to off-site damages including the reduction of water quality within the adjacent river systems. Moreover, the gully-related sediments can be a transport medium for chemical, physical and biological pollutants. This strongly impairs freshwater ecosystems and their ecosystem services (e.g. Kroon et al. 2016) and leads to increased costs in physical, chemical and biological water treatment in order to provide drinking water as well as water for industrial (e.g. cooling, hydropower) and agricultural use (irrigation) (e.g. Zabihi et al. 2018).

Even though gully erosion has a major quantitative impact on a catchment’s sediment yield (e.g. Marker et al. 2003), it is rarely taken into consideration by erosion models (Bull and Kirkby 1997, Poesen et al. 1998, Vanmaercke et al. 2021). One main reason is the high spatial and temporal heterogeneity of gully erosion process dynamics (Sidorchuk et al. 2003, Sidorchuk 2021). For this reason, the mapping of gully erosion is relevant for the identification of hazards and, if necessary, initiation of soil conservation measures. Therefore, many studies discuss the occurrence, modelling and prediction of gully erosion events (see Vanmaercke et al. 2021). The prediction and modelling of soil erosion is crucial for decision-making for soil conservation and landscape planning (Chaplot et al. 2005). Due to the knowledge gain provided by simulation models, many models have been applied to map, simulate and predict gully erosion processes as well as related geomorphic features.

In terms of a quantitative assessment of gully erosion, different modelling approaches have been proposed. In general, the potential location of a gully can be identified with probabilistic models (e.g. Maerker et al. 2011, Angileri et al. 2016; Meliho et al., 2018; Azrah et al. 2019, Nhu et al. 2020, Lei et al. 2020, Cama et al. 2020, Maerker et al. 2020). The idea of these models is to produce a map-based classification of gully erosion susceptibility at local or regional scale depending on the most important predictor variables triggering gully erosion in the study areas. Empirical models offer an alternative approach by predicting retreat rates of a gully headcut erosion (e.g. Vandekerckhove et al. 2001; Marzolff et al. 2011, Poesen et al. 2011, Frankl et al. 2012, Li et al. 2015). However, only a few studies have focused on gully widening and depth erosion (Archibold et al. 2003). Moreover, empirical models are limited in terms of their predictive capacity. Consequently, there is a strong need for process-oriented models that are able to assess most of the gully-related processes, such as headcut retreat, gully-widening and gully-deepening.

Moreover, the process-oriented models should be simultaneously also sufficiently flexible to handle different spatio-temporal scales. This is even more essential for the assessment of both expected future climate and land use changes. Only a few processes-based models such as DIMGUL/STABGUL (Sidorchuk 1999, 2021, Sidorchuk et al. 2003) are capable of simulating gully evolution dynamics. DIMGUL/STABGUL models allow the simulation of gully morphology considering the causes for their development. Hence, a stable gully morphology can be predicted (Zorina 1979, Mirtskhulava 1988), as well as the conditions supporting the incision of ephemeral gullies (Poesen and Govers 1990), gully head growth rate (US Soil Conservation Service 1966, Trofimov and Moskovikin 1983) and the transformation of the length profile (Sidorchuk 1996).

The process of gully formation needs precipitation events that lead to an infiltration surplus or saturation surplus and produces a water discharge that exceeds the runoff threshold required to remove and, respectively, transport sediments. Additionally, topographic circumstances including slope, exposition and catchment area control the amount and volume of water runoff. Depth and morphology of the gully cross-section are dictated by soil erodibility as well as the lithostratigraphic characteristics of the area. Kosov et al. (1978) explained the mechanisms of gully development dynamically and statically pointing out that there are two stages of gully development leading to two types of gully erosion models: 1) the dynamic model predicting rapid changes of gully morphology in the first stage of gully development covering only ca. 5% of the entire gully lifetime, and 2) the static model that calculates the final morphometric parameters of stable gullies covering the rest of the gully lifetime (see also Sidorchuk 1999). Hence, the two models are used to predict both the development of gullies and their final morphologies.

Traditionally, the modelling of gully and channel erosion in GIS environments has raised issues related to computing capacity and efficiency, specific input data requests and inflexible programming interfaces. Generally, GIS is acknowledged as a versatile tool for the analysis and visualisation of spatial model results, incorporates sophisticated statistical analytical tools and offers computational solutions for model calibration and validation (Poeter et al. 2005, Mitasova et al. 2013). Recently, GIS has become increasingly relevant for spatial environmental modelling and specifically for gully modelling, as it offers Application Programming Interfaces (API) and supports script languages such as Python, and libraries for raster algebra. A model that has been implemented in a GIS environment can be used as a tool via the user interface and the input and output are stored in standardised geodatasets. Therefore, there is no need for transformation between different, model-specific formats. Further advantages of such module development include customisable development tools, integration of built-in raster calculators (Wesseling et al. 1996) and visual programming languages.

The objective of this study is to simulate processes of gully incision and thus, to estimate the gully evolution over longer time periods using a Python-based tool integrated in a GIS-environment. This simulation is based on the gully models proposed by Sidorchuk and Sidorchuk (1998) and Sidorchuk (1999, 2003) and was originally provided in the form of a Fortran code. Sidorchuk’s models allow the simulation of the gully evolution considering different lifetime stages (see above). The main aim of this work is to implement a gully simulation model with the respective pre- and postprocessing procedures in a GIS-environment and provide a 2D/3D visualisation tool to illustrate the obtained results.

2 The DIMGUL theoretical framework

Gully formation originates under certain conditions when the force of the water flow exceeds a certain threshold value, and the shallow water flow begins to concentrate along flow lines. Particularly, these thresholds are reached at knickpoints or at depressions along the slope. Several flow lines may merge and increase the probability for gully development, eventually merging the knickpoints and forming a single incised stream (Haan et al. 1994). The gully form then progresses via headwall migration and the channel widens.

As described by Sidorchuk (1999), the first phase of gully initiation is intense and characterised by the instability of morphological characteristics like volume, area, length, width and depth. Instead, the last stage describing the stable gully evolution is associated with sediment mobilisation and deposition at the gully bottom. Moreover, lateral widening and mass movements of the gully side walls are observed. Empirical measurements by Kosov et al. (1978) demonstrate that the first stage evolves rapidly in the first 5% of the gully’s lifetime, but accounts for 60% of its area and 35% of its volume (Figure 1). In the remaining lifetime gully area and volume are still subject to change, but under more stable conditions.

Figure 1. Evolution of the gully morphology through time: after (Kosov et al. 1978).

Sidorchuck’s model of gully erosion considers the following fundamental processes during the dynamic initial stage: Water flow induced by substantial precipitation events excavates a rectangular channel in the upper soil. Destabilisation of the vertical incision walls triggers shallow mass movements and transforms the rectangular cross-section into a trapezoidal form during the time period between relevant precipitation events. The gully erosion rate is affected by multiple factors including surface runoff flow velocity, incision depth, turbulence of water flow, air and water temperature, soil texture, soil mechanical characteristics and the influence of a protective vegetation cover.

(1) $\frac{\mathrm{\partial }Z}{\mathrm{\partial }t}-a\frac{dZ}{dx}=0\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{s}\mathrm{p}\mathrm{o}\mathrm{r}\mathrm{t}\mathrm{e}\mathrm{q}\mathrm{u}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}$(1)

With

(2) $a=k_E\mathrm{H}\frac{Q}{W}$(2)

Where,

$k_E$ = longitudinal erosion coefficient

W = flow width

Q = water discharge

H is 1 when flow velocity U is more than critical velocity of erosion initiation U_cr and H = 0 when U ≤ Ucr

The complex interaction between flow and soil cohesiveness in the gully section is modelled with a rather simple transport calculation given by Equation 1(1) $\frac{\mathrm{\partial }Z}{\mathrm{\partial }t}-a\frac{dZ}{dx}=0\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{s}\mathrm{p}\mathrm{o}\mathrm{r}\mathrm{t}\mathrm{e}\mathrm{q}\mathrm{u}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}$(1) for bed elevations Z along the main flow line in direction x. This transport equation (Equation 1(1) $\frac{\mathrm{\partial }Z}{\mathrm{\partial }t}-a\frac{dZ}{dx}=0\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{s}\mathrm{p}\mathrm{o}\mathrm{r}\mathrm{t}\mathrm{e}\mathrm{q}\mathrm{u}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}$(1) ) considers the feedback phenomenon of the relief transformation on the erosion rate and it takes into account a rather limited number of input characteristics such as the initial longitudinal profile for each flow line, the specific water discharge distribution along this line, the erosion coefficient k_E and the associated critical velocity of erosion initiation for each soil texture (Equation 2(2) $a=k_E\mathrm{H}\frac{Q}{W}$(2) ). Gully channel transformation is very intensive, and morphometric characteristics of the gully (length, depth, width, area, volume) are far from stable and changing rapidly. This model was described in more details; see Sidorchuk (1999, 2015, 2021).

3 Study area and environmental data

The model framework described above was tested in the Drakensberg Foothills situated in the KwaZulu-Natal Province of South Africa (Figure 2). This region was selected because large parts of KwaZulu-Natal show a moderate to extremely high risk of erosion (90% of its area), while the currently intact grass cover keeps the actual erosion risk at much lower levels at around 18% of the area (Le Roux et al. 2008). However, the combination of erodible soils and extended subsistence agriculture around informal settlements threatens the current state of stability (Rooseboom et al. 1992, Fluegel et al. 2003). Our study area, thus, represents a typical example of donga erosion in such an environment with a pronounced research history.

Figure 2. Location of the gully Kwa Thunzi – South Africa.

Following the geomorphic classification by Partridge et al. (2010), the area belongs to the Eastern Coastal Hinterland with a Cwb climate following the Koeppen-Geiger classification, which is characterised by a warm temperate climate, with warm summers and dry winters (Conradie 2012). Mean monthly temperatures measured at Himeville (distance 20 km) range between 7°C and 19°C, whereas mean annual rainfall amounts to 935 mm and is mainly falling in the austral hemisphere’s summer months (Nel 2009). The KwaThunzi Gully (Lat −29.6164°, Lon 29.6480°) is located at an elevation of ~1200 m a.s.l. on a slope exposed to southwest covered by grassland in the upstream area of the Mkomazi River. The gully developed in colluvial sediments of the Masotcheni formation that is of Late Pleistocene age (Botha et al. 2016, Bosino et al. 2021). The underlaying basement belongs to the Beaufort Group shale stones that is disturbed by a network of dolerite dykes in turn forming the local erosion bases. The gully walls expose a sequence of paleosols and colluvial deposits of the Masotcheni formation up to a thickness of 9 m comprising a basal unstratified layer, followed upwards by a transitional layer and several paleosols with columnar structures, iron-manganese and carbonate concretions, as well as leaching zones. The overlying recent Solonetzic soil is developed to a thickness of 30 cm characterised by a columnar structure and a Clay-Loam texture (USDA: ClLo).

The main model input data comprise 1) the terrain morphology (digital elevation model), 2) the soil and lithological composition (substrate and soil mechanical parameters such as texture and hydraulic conductivity) and 3) meteorological information. The topographic data, in terms of a digital terrain model (DTM) of the KwaThunzi gully, was created using a structure-from-motion approach. The images were acquired in the field on 13 January 2018 by a UAV equipped with a Canon PowerShot S110 camera. For the 3D processing the Mavis (v.17.09.00) software by BitMapping was used. Initially, the camera images were matched with their GPS-coordinates based on the time tag in both datasets. Subsequently, images exceeding a maximum threshold in yaw, pitch, or roll were excluded. This allowed the alignment of images and tie-point generation. The majority of images overlap with at least 10 neighbouring images. Subsequently, a gridded DTM has been generated with a resolution of 10 cm and an RGB orthomosaic that features a spatial resolution of 5.5 cm. The DTM was resampled to a spatial resolution of 1 m as a compromise between the level of detail and respective computational costs.

Surface runoff data have been derived from discharge measurements at a gauging station close to the study area at Lundy’s Hill (Lat −29.9098°, Lon 29.7398°) with a data time series of more than 60 years. The lithology and substrates were investigated in the field using selected sections that were chosen based on the geological map (see Bosino et al. 2021). The lithological section of the study area was finally classified into four soil/substrate units which are, topsoil, paleosol, transitional colluvial layer and the colluvial base layer. The main soil parameters that are used in the model are calculated for each layer such as: i) critical velocity values following Bogomolov and Mikhaylov (1972), ii) stable slope that was determined in radians for each unit and iii) the erodibility coefficient (K) after Renard et al. (1997).

4 Methodology

Implementing the gully erosion model in a GIS system is quite complex and involves the conceptualisation of a workflow that includes the following steps: 1) processing of input data, 2) model calculation, 3) analysis and 4) visualisation of results. The workflow is then implemented using a visual programming tool that is integrated in GIS (model builders). We used python as programming language that supports fully automated simulations (Schaller et al. 2009). Complex physical-based models such as the chosen gully models usually are based on partial differential equation solutions. However, these equations can be fully implemented in a GIS framework using the system’s geospatial data management application programming interface (API).

The rate of gully incision is controlled by water flow velocity, channel depth, as well as by soil texture, soil mechanical pattern and the degree of protection by vegetation (Sidorchuk 1999).

Based on the elevation and distance of specific points from the gully mouth along the longitudinal profile of the flow line (including existing gullies) gully morphology is calculated.

Therefore, also the change in water discharge over time was estimated for all these points using surface-specific runoff derived from the discharge data available at Lundy’s Hill. The empirical relationships (see Equations 1(1) $\frac{\mathrm{\partial }Z}{\mathrm{\partial }t}-a\frac{dZ}{dx}=0\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{s}\mathrm{p}\mathrm{o}\mathrm{r}\mathrm{t}\mathrm{e}\mathrm{q}\mathrm{u}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}$(1) and 2(2) $a=k_E\mathrm{H}\frac{Q}{W}$(2) ) are used to calculate flow width, flow depth and flow velocity. Soil or substrate properties are used in the model comprising critical velocity, slope stability and erodibility coefficient as described above. Table 1 shows the input parameters used in the modelling.

Table 1. Modelling parameters.

Parameters	Value/Formula	Source
Soil Parameters
Critical Velocity (m/s)	0.7	Maerker (2001)
Slope Stability (Rad)	0.59	Maerker (2001)
Erodibility Coefficient (K)	0.004	Maerker (2001)
Hydrological Parameter
Daily Specific Discharge (q)	’Specific discharge’: Conversion to drainage per area [mm d^−1 m^−2]:-$q=\frac{Q}{A}$q = runoff per area (l s^−1 km^−2) Q = runoff (m^3s^−1) A = catchment area (km^2)	(Baumgartner & Liebscher 1996)
Channel Parameter
Channel Width	$W=3.0\ast Q^{0.4}$	Sidorchuk (1999)
Channel Depth	$\mathrm{D}=0.48\ast Q^{0.45}$	(Sidorchuk 1999)
Topographic Parameters
Flow lines extraction	ModelBuilder	(Schröder and Omran 2013)
Flow Accumulation	ModelBuilder & Python Script	(Schröder and Omran 2013)
Elevation for each point on flowline	Python script developed in step (4.1)
Coordinates (x,y,z) for each point along flow line

The complete simulation can be divided into four main steps (Figure 3): 1) identification and delineation of the flow path network from the digital elevation model (DEM), 2) preparation of the topographic data in terms of selecting the flow path that is representing the main stream flow of the gully, 3) calculation of the main parameters characterising the resulting gully morphology such as length, width, depth and volume for different time steps (e.g. yearly time step) and 4) creation of 2D and 3D representations of gully incision for the established time period.

Figure 3. Workflow scheme of the GIS Model of Gully erosion.

4.1 Flow path extraction

The aim of this step is to extract the flow line paths from the digital elevation model. Therefore, the concept of Strahler’s stream ordering technique was applied to perform a qualitative and quantitative assessment of the drainage system. The flow lines can be extracted in a GIS-environment using a specific automated tool for ordering, rasterisation and vectorisation flow lines (Schröder and Omran 2013, Omran et al. 2016). This tool is organised in the following steps as illustrated in Figure 4: 1) Creation of a DTM (if not available for analysis), 2) Application of a fill sink procedure to preprocess the DTM, 3) calculation of local flow directions, 4) calculation of flow accumulation, 5) selection of the optimum threshold value for stream line extraction, 6) rasterisation of flow lines based on a Strahler order segmentation and 7) vectorisation of the flow lines.

Figure 4. Automated stream network extraction using the GIS model (Schröder and Omran 2013).

4.2 Derivation of topographic input information

In the second step, the required topographic data is prepared that is used in the model simulation step. The topographic files contain the elevation data along the flow lines at discretised points as well as the corresponding values of the flow accumulation that are extracted in the first step together with the thickness of each soil/substrate layer measured in the field. Additionally, a file containing the global coordinates of each segmentation point along the selected flow line is generated including latitude (x), longitude (y) and elevation (z) for each point. The model uses Esri’s ArcGIS Model Builder, including several customised Python tools as shown in Figure 5.

Figure 5. Full model for topographic data preparation, which include three parts.

5–1: Procedure to convert flow lines to points and attributing flow accumulation. 5–2: Procedure for routing process and extracting the input table. 5–3: Procedure extracting the coordinates of flow line points

For the elevation data of each layer along the flow line, two digital elevation models were used. The first DTM represents the initial situation of the area before gully erosion starts and the second DTM which reflects the situation of the study area after gully erosion, or at a specific time step during gully evolution. The initial DTM can be reconstructed by refilling the gully trench using the fill module applied on the present-day DTM of the study area or creating an initial DTM, e.g. from stereo aerial photos with high resolution. For the data preparation, we implemented two Python scripts within full ModelBuilder including different tools. The first part of the model discretised the selected flow line into point geometries and calculates the distance of each point along the flow line to the gully mouth (i.e. the outlet of the gully basin) (Figure 5–1).

The second part of the model includes the first script (Figure 5–2) to make routing process for all the points along the flowline and to prepare the tables including the geographic coordinates (Easting and Northing values) and the elevation points from the initial digital elevation model along the flow line. This information will be used later to embed the simulation results given as local coordinates into the DTM and thus the real-world spatial reference. The number of soil/substrate layers and their thicknesses can be added by the user in order to calculate the elevation of each layer along the point geometry. Finally, the third part of the model is a script (Figure 5–3) to generate the topographic file including the flow accumulation as well as the elevation of each point and geographic coordinates along the surface of each layer in the form of a text file as input for the simulation.

4.3 Gully simulation

This step aims at the calculation of the main gully morphometry. The script computes for given time steps the main geometrical parameters of the area that are affected by gully erosion such as the change in incision depth, the change in top- and bottom width of each layer, which is represented by the point features along the gully flow line, as well as the volume of the gully erosion and the cross-sectional area along the flow lines. All the mandatory input data for the gully model such as runoff data and the soil parameters (critical velocity, slope stability, and erodibility coefficient) for each layer, as well as the topographic information from the previous step before are passed to the Python script of the model.

Altogether four Python functions were used: 1) a function defining the calculation environment, i.e. the input parameters such as time steps for the numerical solution of the differential equations as well as for the output interval, topographic parameters, runoff and soil parameters; 2) a function implementing the numerical solution of the main equations based on Sidorchuk’s model; 3) a function calculating the soil parameters such as slope stability, critical velocity and erodibility coefficient for each point along the selected flow lines as well as the wall correction equation from the cuboid form to trapezoidal form; and, finally, 4) a function that writes the output parameters in respect to the defined interval.

4.4 Gully visualisation

The last step regards the visualisation of the model output in a GIS environment. The aim of this stage is to embed the calculated gully depths and widths values along the flow line into a coordinate system in order to be visualised in a 2D and/or a 3D visualisation format, together with the related topography. This embedding is achieved by calculating 3D polygon features representing the gully bottom and the gully walls between adjacent flow line points for each time step, e.g. for each year. These polygons were created based on the defined coordinate file of the study area calculated in the data preparation step, and the geometrical parameters that were calculated in the simulation. Using simple geometric calculations implemented in a Python script, the local coordinates along the flow line for each layer are transformed into coordinates (Easting and Northing values) related to a geographic coordinate system. Additionally, the script creates 3D polygon geometries based on the GIS API for each year. The polygons can be projected and merged to define the upper boundary of gully erosion. As an output of this script, we achieve sequences of polygons that represent the bottoms and side walls of the stream flow lines for each layer after the erosion process at each time steps (e.g. yearly). The creation of these sequences of polygon geometries is used in the further 2D and 3D output visualisation of gully erosion in the GIS environment.

5 Results and discussion

Before starting the model processing steps, all relevant geo-data were projected in the spatial reference system UTM zone 35 South related to the geodetic datum WGS84. The flow lines have been extracted from the original DTM of the study area as shown in Figure 6. The resulting flow lines were evaluated visually by overlaying the RGB orthomosaic of the study area. The extracted main flow line coincides with the main path of the gully development. Therefore, this main flow line was selected to calculate the gully’s geometry. The model preparation step was run on the selected main flow line. This flow line was transformed into a sequence of point features. We assigned to each point the elevation of each soil or substrate layer, the flow accumulation and the coordinates of each point.

Figure 6. Extraction of streamflow lines from high resolution DEM (15 cm).

In the gully simulation step, different parameters were calculated such as the volume of gully erosion, the change of the cross-sectional area, and the change of the length of the gully corresponding to the number of years used in the simulation. During precipitation events, water runoff excavates a rectangular shaped channel cross-section in the topsoil or at the gully bottom. The incision is mainly dictated by the imbalance of detachment of particles from and the sedimentation on the gully bottom (Sidorchuk 1999).

The gully model in the end allows the calculation of the final gully morphology. The mechanism of incision is mainly related to the actual flow velocity that must be higher than the critical velocity for sedimentation. The shape of the longitudinal profile of the gully is decoupled from the initial slope, and the evolution of the gully is controlled only by discharge parameters like critical velocity and soil texture.

Particularly, the critical velocity of the soils and substrates, as well as the extent of runoff discharge, primarily affect the speed of incision of the gully (Sidorchuk et al. 2001). For the implemented model, the erosion rate of the first year was assumed to be of 0.5 m depth corresponding to the Solonetzic soil characteristics showing high erodibilities in the Btz horizon at 0.5 m soil depth. Figure 7 shows that the length of the gully evolves very fast following the findings of Kosov et al. (1978). Already 80% of gully length is developed after ca. 20 years. The area values show an almost stable behaviour after 60% of gully lifetime. Instead, the evolution of the gully depth and gully volume is much slower and continues up to the end of gully’s lifetime.

Figure 7. The results of the simulation step, the left side represents the text files of the calculated parameters for gully erosion (total depth, total volume, change in top width, change in bottom width and the change in area for each year) and the right side represents the changes of gully in volume, length, depth and area over time.

In the last modelling step, the results are visualised in 2D and 3D visualisations utilising a GIS environment (Figure 8) As Sidorchuk (1999) noted, during the dynamic phase of gully formation, runoff can erode a rectangular channel into the topsoil or gully floor. Thereafter, the vertical sidewalls of the gully may be unstable, resulting in small mass movements and further erosion. Hence, the rectangular cross-sections transform into a trapezoidal shape. In this stage, the simulation assumes straight alignment of finite elements for the flow line. However, the translation of the modelled flow line into the real flow line topography using the geographic coordinates is quite challenging for the visualisation in a GIS environment due to curved bottom features and side wall characteristic.

Figure 8. Simulated Gully bottom evolution during a 50-years period – 2D and 3D view.

The problem with the bending flow line was tackled by calculating line segments perpendicular to the flow line based on the modelled top width. Points with overlapping line segments were skipped in the calculation. If the segments are not overlapping, the coordinates of the points defined by the top and bottom width of the layer are calculated as well as the coordinates of the points defining both side walls at the bottom of the respective layer. These points are used as vertices for polygons representing the boundaries of erosion areas for each time step (e.g. year). The width of the gully and the area of incision were simulated over a 50-years period. Both side walls of the primary water flow can be visualised in a 3D-view.

For a detailed validation of the results and for a sensitivity analysis of the model, different parameters were tested. Particularly, we tested the sensitivity of flow discharge and critical velocity over a 50 years model run, keeping constant all other model settings. Two discharge magnitudes were simulated by just changing the units from l/s to m³/s corresponding to an increase by a factor of 1000. Concerning the critical velocity, we tested two different values representing two different soils: i) a resistive soil with high critical velocity (0.57 m/s) and ii) an erodible soil with low critical velocity (0.32 m/s). We observed that discharge values have a significant effect on gully area and volume evolution. The volume of the gully increases significantly with the increase in the discharge magnitude as shown in Figure 9(a). Additionally, it has also been observed that the resistant soil with high critical velocity (0.57 m/s) has a significant lower total gully volume than the soft soil with low critical velocity (0.32 m/s). Consequently, the highest potential for gully erosion is close to the gully outlet since we have the largest discharge values (see Figure 9(a)). This may lead to incision for low resistive substrates or a broadening of the gully if resistive layers are reached, as shown in Figure 9(b). However, as shown, the simulation model is very sensitive to changes in critical velocity and discharge; therefore, these parameters should be selected carefully.

Figure 9. (a) Effect of changed discharge units on gully erosion. (b) Effect of the changed critical velocity values on gully erosion.

6 Conclusion

Gully erosion is associated with a variety of on-site and off-site damages, particularly affecting agricultural lands. Governments and research institutions around the world have been trying for years to evaluate the hazards associated with gully erosion and to predict their spatio-temporal distribution. The implemented gully simulation code presented in this paper is mainly based on Sidorchuk’s gully model using Python libraries along with GIS APIs. Thus, we propose a modelling tool implemented in a GIS environment to assess the changes in gully morphology over time. Consequently, gully evolution over the full gully lifetime can be simulated. The application uses relevant input data of the study area such as the topographic characteristics, soil and substrate properties as well as surface runoff information. This data is preprocessed and organised to be easily implemented and visualised in the chosen GIS environment. The application is designed to model annual changes in gully morphology. The complete simulation can be divided into four main steps: Step 1) establishment of the flow lines from a digital elevation model; Step 2) preparation of the topographic data relevant to select flow paths that represent gully erosion; Step 3) calculation of the main morphometric parameters such as length, width, depth and volume per each segment along the flow path that finally allows for the calculation of soil erosion for each year; and Step 4) visualisation of the main morphometric features and dynamics of gully erosion in a GIS.

The application is able to handle many different layers of soil or substrates. It must be noted that processing time depends on the number of flow lines, number of years for discharge data and the erodibility coefficient of the soil layer types. The erodibility coefficient is used as a Courant number in the numerical solution of the underlying differential equations defining the internal time step. Thus, the higher the value, the more time is needed to incise or erode a certain layer or to shape the bottom width of a gully channel. The quality of the generated model depends mainly on the accuracy of the input parameters, e.g. soil characteristics, discharge data and topographic data derived from DTM. The resulting morphology can be analysed to evaluate the amount of soil loss per year.

The proposed tool was tested for a gully system of the KwaThunzi area located in the upper Mkomazi River basin in KwaZulu Natal, South Africa. The results were very promising through the visualisation of the time series of the gully’s development. However, there are still challenges to be faced, especially, related to the use of empirical parameters in the model. Especially, critical velocity and discharge are two sensitive parameters that need a careful calibration. Additionally, the simulation algorithm assumes straight alignment of finite elements. Improvement of the algorithm should consider the bending of the real flow line. Finally, the application should be tested in different geographical regions and climatological conditions, so that a better understanding of the role of the input parameters may be achieved.

Authorship contribution

Adel Omran concept, model developer, analysis, writing, editing; Dietrich Schröder model developer, writing, editing; Christian Sommer Data provider, fieldwork, writing and editing; Volker Hochschild data provider, and editing; Michael Maerker concept, data provider, writing, editing and fieldwork.

Acknowledgments

The authors would like to thank profoundly Prof. Dr Alexey Sidorchuk for providing the code of the gully models and his fundamental support and cooperation throughout the entire process of implementing the model and writing this paper. Moreover, the authors would like to thank the Department of Photogrammetry and Geoinformatics of Stuttgart University of Applied Sciences, the Department of Geography of University of Tübingen, the ROCEEH project (The Role of Culture in Early Expansions of Humans) of the Heidelberg Academy of Sciences as well as the Department of Earth and Environmental Sciences of Pavia University for providing relevant data and support in conducting this study.

Disclosure statement

Conflicts of interest the authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.