Soil erosion is recognized 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 are eroded by wind and water related processes often impacting agricultural land . Soil erosion by water, such as gullying, is acknowledged as a global main driver of land degradation . Gully erosion is typically defined as a deep channel that has been eroded by concentrated water flow, removing surface soils and materials . 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 [3, 8]. Gully erosion is related with 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 infrastructure such as roads, power lines, and communication networks [15, 13, 14]. 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 20 to 80% of a catchment’s mean sediment yield (e.g., [19, 46]). 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. ) 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. ).
Even though gully erosion has a major quantitative impact on a catchment’s sediment yield, (e.g., ), it is rarely taken into consideration by erosion models [8, 30, 46]. One main reason is the high spatial and temporal heterogeneity of gully erosion processes dynamics [40, 38]. 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, modeling, and prediction of gully erosion events (see ). The prediction and modeling of soil erosion is an important basis for decision-making for soil conservation and landscape planning . 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 assessing gully erosion quantitatively, different modeling approaches have been proposed. In general, the potential location of a gully can be identified with probabilistic models [21, 22, 1]. Empirical models offer an alternative approach by predicting retreat rates of a gully headcut erosion (e.g., [45, 30, 11, 18]). However, only a few studies focused on gully widening and depth erosion . 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 approach 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 models such as DIMGUL/STABGUL [37, 40, 38] are capable of simulating gully evolution dynamics. DIMGUL/STABGUL models allow the simulation of gully morphology considering the causes for their development. Hence, stable gully morphology can be predicted [49, 23], as well as the conditions supporting the incision of ephemeral gullies , gully head growth rate [43, 44], and the transformation of the length profile .
The process of gully formation needs a 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 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 mechanism of the 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; 2) the static model that calculates the final morphometric parameters of stable gullies covering the rest of the gully lifetime (see also ). Hence, the two models are used to predict both the development of gullies and their final morphologies.
Traditionally, the modeling of gully and channel erosion in GIS environments 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 visualization of model results, incorporates sophisticated statistical analytical tools and offers computational solutions for model calibration and validation [32, 24]. Recently GIS is becoming 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 within the user interface and the input and output is stored in standardized geodatasets, so there is no need for transformation between different, model-specific formats. Further advantages for such mode development include customizable development tools, integration of built-in raster calculators  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 (1998, 1999, 2003) [41, 37, 40]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 for this work is to simulate gully erosion in 2D and 3D within a GIS-environment, either using a commercial or an open-source GIS-environment. The GIS-based gully erosion model was set up and run based on a study area located in KwaZulu-Natal, South Africa characterized by large gully systems developed in colluvial material also called Masotcheni formation . Subsequently, the simulation results were compared with the existing gully morphology. Finally, the gully model modelling framework allows also the reconstruction of the initial environmental conditions before gully incision.