Models of Gully Erosion by Water

Abstract

The type of modelling of gully erosion for the projects of land management depend on the targets and degree of details of these projects, as well as on the availability of input data. The set of four models cover a broad range of possible applications. The most detailed information about predicted gullies, change of their depth, width, and volume throughout the gully lifetime is obtained with the gully erosion and thermoerosion dynamic model. The calculation requires the time series of surface runoff, catchment relief, and lithology and the complex of coefficients and parameters, some of which can be estimated only by model calibration on the measurements. The difficulty in obtaining some of these coefficients makes it necessary to use less complicated models. The stable gully model predicts final gully depths and widths and is useful for projects where only stable gully geometry is used. The modified area–slope approach is used in the two simplest models, where the position on the slopes of possible gullies is calculated without details of the gully geometry. One of these models calculates total erosion potential, taking into account all water runoff transforming a gully. The second calculates gully erosion risk, using the information about slope inclination, contributing area and maximum surface runoff.

1. Introduction

The need to assess possible gully position, depth, width and volume in agricultural areas and areas of new development is well known [1,2,3,4]. This practical necessity of land management causes the emergence of a wide variety of methods for assessing gully erosion potential (see reviews [5,6]). However, in this diversity, there is an obvious bias towards using the same theoretical approach, determining the critical slope and catchment area for the formation of a gully (area–slope approach), proposed in [7,8]. Often it is enough for practical purposes to estimate the points on the initial slope, where flow achieves threshold conditions and gully erosion initiates. The area–slope approach contains a significant empirical component (empirical coefficients) and the use of such models is restricted to the region of empirical data collection ([9,10,11] and many others, see bibliography in [6]). Most existing models of gully erosion, more or less empirical, find such threshold points and predict a maximum length of gully. The limitation of this approach and ways to overcome them are discussed in this paper.

It is not possible to estimate other measures of the gullies, such as depth, width, area and volume within the area–slope approach. Therefore, additional empirical models are used to calculate, for example, a gully width [12,13].

Maximum morphometric measures are not achieved simultaneously during the formation of a gully. Observations and experiments show [14] that a gully formation has two main stages. The first is the stage of rapid gully erosion, when during about 5% of the gully lifetime, more than 90% of the gully length, more than 80% of its depth, and about 60% of the area is reached. The second is the stage of gully stabilisation, when the area and last of all, the volume increases. Therefore, models of gully erosion calculation naturally fall into two groups for the prediction of gully morphology at these two stages.

It is possible to use the area–slope approach to calculate threshold conditions of flow not only in the gully head points but also in all points along the gully longitudinal profile. This way led to the calculation of the longitudinal profile of the stable gully, at each point of which flow velocity (controlled by slope, contributing area, resistance to flow and runoff depth) is equal to the critical velocity of erosion initiation. The processes of gully head and wall transformation by gravitation become important in this stage. Zorina [15] proposed one of the first models of this type, and its modification is described further. The full development of the stable gully takes many hundreds of years [16].

Most sophisticated dynamic models of gully evolution are based on the solution of an equation of mass continuity [17,18]. Usually, dynamic models describe the first stage of rapid gully formation, although such models can be run up to the stage of gully stabilization.

Different types of gully erosion models are used for different purposes, depending on the targets of such calculations and on the existence of input data. We offer in this paper the set of models, which cover the range of such targets in land management projects. This set includes four models of gully erosion prediction, where all known approaches to solving this problem are applied.

Two models calculate gully depth and width along its longitudinal profile:

(1)

The dynamic model of gully formation in “real” time. Not only mechanical erosion by water on a gully bed but also erosion by thermal action of water (if frozen soil is present) and gravitation processes on the gully walls are included in this model.

(2)

The model of finite stable gully describes gully geometry when water flow and gravitation already do not deform the gully bed and walls.

In the next two models, the novel modification of the area–slope approach is used, which gives the most probable position of possible gullies:

(3)

The model of gully erosion potential, which takes into account the effect of all runoff from the catchment.

(4)

Express model of the risk of gully erosion, which takes into account only the maximum surface runoff.

The main goal of this research article is to describe these four gully erosion models of different levels of simplification, to show their similarities and differences, and to compare the results of calculations of the characteristics of the same gully using these four models. The sequential comparison of the calculation results obtained with different models for the same object is novel and allows identifying the advantages and disadvantages of the proposed models and developing recommendations for their use in practice. This comparison makes it possible to choose which of the models gives the best results for available input data and for the level of complexity of the land management project.

2. Methodology—Model Descriptions

2.1. The Models with Gully Longitudinal Profile Calculations

2.1.1. The Dynamic Model of Gully Erosion

The dynamic gully erosion and thermoerosion model GULTEM (or DYNGUL model for erosion only) was proposed for calculations at the first stage of gully evolution. At this stage, the erosion and thermoerosion are predominant at the gully bottom, and rapid mass movement occurs on the gully walls. 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 is described in detail elsewhere [17,18], therefore, only basic information is given here.

The two main processes to be described are:

(a)

Formation of a rectangular cut in the initial slope by mechanical and thermal action of the flowing water.

(b)

Transformation of the rectangular cut into a trapezoidal shape by shallow landslides during the period between adjacent runoff events.

The rate of gully incision is controlled by water flow velocity, depth, temperature, as well as by the soil mechanical pattern and the level of protection by vegetation. These characteristics are combined in the equations of mass conservation:

$$\frac{\partial Q_s}{\partial X}=C_w{q}_w+EW+E_{b}D-C{U}_{f}W$$

(1)

and deformation:

$$\left(1-\epsilon \right)\frac{\partial V}{\partial t}=\frac{\partial Q_s}{\partial X}$$

(2)

where Q_s = Q C is sediment transport rate (m³ s⁻¹), Q—discharge (m³ s⁻¹); X—longitudinal coordinate (m); t—time (s); C—mean volumetric sediment concentration; C_w—sediment concentration of the lateral input; q_w—specific lateral discharge (m² s⁻¹); E—gully bed erosion rate (m s⁻¹); E_b—channel banks erosion rate (m s⁻¹); V—gully “empty” volume (m²); W—flow width (m); D—eroded bank height (m); U_f—particles fall velocity in the turbulent flow (m s⁻¹), ε—soil porosity.

During the episode of erosion, the accumulation of sediments on the gully bed is assumed negligible. Therefore, Equations (1) and (2) can be transformed to:

$$\frac{\Delta V}{\Delta t}=EW+E_{b}D$$

(3)

The analysis of the experiment results in the natural gullies in different environments and in experimental flumes [17,18] shows that in the conditions of steep slopes and cohesive soils, common for gullies, the mean rate of bed erosion is linearly correlated with the product of bed shear stress τ = g ρdS and the mean flow velocity U (i.e., stream power):

$$E\sim \tau U=k_{E}HqS$$

(4)

where S is flow surface slope, d is flow depth (m), q = Q/W is specific discharge, g is acceleration due to gravity, ρ is water density and k_E is the combined erosion coefficient, which depends on soil properties. H is Heaviside step function, equal to 0 when flow velocity U is less than the critical velocity of erosion initiation U_cr, and equal to 1 when U ≥ U_cr. Critical velocity depends on soil and vegetation cover properties [19].

Equation (4) links the dynamic model with the methods of gully potential estimation based on the threshold slope-contributing area approach [8,9].

In the case of gully erosion in the frozen soil (thermoerosion), the water temperature becomes the main factor of erosion. Field and laboratory experiments [18] showed that as the first approximation, the soil detachment rate is linearly correlated with water temperature T °C and equal to the rate of thawing of frozen soil with the open surface:

$$E_T=k_{TE}T$$

(5)

where k_TE is the coefficient of thermoerosion. Equation (5) is used for calculation only for the case of direct contact of water with frozen soil, without any thaw layer. Therefore the inequalities

$$E<E_T\mathrm{or}E>E_T$$

(6)

are checked in the calculation algorithm. If the rate of thermal erosion is less than the rate of mechanical erosion $\left(E>E_T\right)$ and a thaw layer does not form. Then Equation (5) is used for estimation of the gully bed lowering by thermoerosion. On the contrary, if the rate of the thermal front movement in soil is more than the rate of mechanical erosion $(E<E_T)$ then a thaw layer is formed, and Equation (4) is used in the model for the mechanical erosion rate calculation.

The erosion rate of the gully banks can be estimated only in a first approximation as some function of the rate of gully bed erosion, controlled by the ratio of lateral velocity v and longitudinal velocity U. Using this assumption, an approximate formula for the rate of eroded bank erosion can be proposed:

$$E_b=k_{b}E$$

(7)

where k_b is the bank erosion coefficient, the details see in [18].

Flow width W and depth d in the gullies usually is calculated via regime equations as the power functions of discharge:

$$W=pw{Q}^{mw}$$

(8)

$$d=p{Q}^m$$

(9)

Flow velocity U is

$$U=\frac{Q}{Wd}$$

(10)

or, according to the Chezy–Manning formula

$$U=\frac{\sqrt{S}}{n}{d}_{}^{2/3}$$

(11)

where n is Manning’s roughness coefficient.

Gully walls become practically straight after rapid sliding, following the incision. In this case, a model of straight slope stability can be used for the prediction of gully walls inclination. It is possible for the practical needs to measure the gully walls inclinations at the stable sections and use the measured angle φ for the further calculations.

When the bottom width, wall inclination φ and whole volume of the gully V are known, the shape of the gully cross-section can be represented as a trapezoid with bottom width W_b, depth

$$D_g=\left[\sqrt{W_b^2+\frac{4V}{\tan \left(\phi \right)}}-W_b\right]\frac{\tan \left(\phi \right)}{2}$$

(12)

and top width

$$W_t=W_b+2\frac{D_g}{\tan \left(\phi \right)}$$

(13)

The result of the calculations was gully depth, bed and top width, volume for each cross-section of the gully along each flowline and for each time step. The calculations required information about initial relief of the gully catchment and about boundaries of all lithological units, including topsoil with vegetation in the form of DTM (Digital Terrain Model), the sequence of surface runoff values, and all parameters and coefficients used in Equations (4)–(13). The coefficients are empirical and local and must be estimated for a given territory. The calibration of the model parameters, especially the erosion/thermoerosion coefficients and critical velocity on measured data, is highly recommended.

The dynamic model of gully erosion is the only one with which gully evolution can be simulated in “real” time, and all details of this process (within postulated assumptions) can be adequately described. The main problem was the simplification in the model of the real process (such as straight slopes assumption) and, nevertheless, a large number of parameters and coefficients, some of which required model calibration. The level of simplification of all other models, described further, was significantly greater, and the result of the modelling was less definite. However, these models require less initial information and can be useful for preliminary investigations.

2.1.2. Stable Gully Model

The stability of the gully refers to such a stage when the gully bottom is not eroded by water, and the gully walls do not change their shape. This requires two conditions: (1) flow velocities U in the gully are less than critical velocities of erosion initiation U_cr for soils composing the gully bed; (2) gravitation processes become negligible and the gully walls have reached the “angle of repose” for the given soil texture, cohesion and wetness.

In the STABGUL model, the critical velocity of erosion initiation U_cr was calculated using the Chezy–Manning formula (Equation (11)), and the flow depth d was represented as a power-law function of discharge (Equation (9)), which was calculated as the catchment area F multiplied by the runoff depth M:

$$d=p{\left(kMF\right)}_{}^m$$

(14)

where k is the coefficient for changing the dimension of the quantities included in the formula. Then

$$U_{cr}=\frac{\sqrt{S_{}}}{n}{d}_{}^{2/3}=\frac{\sqrt{S}}{n}{p}^{2/3}{\left(kMF\right)}_{}^{2m/3}$$

(15)

Therefore, the stable gully bed inclination is:

$$S=\frac{{\left(n{U}_{cr}\right)}^2}{p^{4/3}{\left(kMF\right)}_{}^{4m/3}}$$

(16)

Equation (16) links the model of the stable gully to the methods of gully potential estimation, based on the threshold contributing area-slope approach [8,9]. It can be written in the form of an ordinary differential equation:

$$\frac{dZ}{dL}=a\left(L\right)F{\left(L\right)}^{-b\left(L\right)}$$

(17)

where the contributing area F, coefficient a and exponent b are functions of length L from the mouth of the gully with the initial altitude z₀₀ to i-th point on the flowline with initial altitude z_0i. Equation (17) is solved numerically along all flowlines with the known altitudes of the mouth of the stable gully z₀₀, runoff depth M, critical velocity, bed roughness, parameters in Equation (16), and the required functions of length. The result of the solution is the partial altitudes z_ij of gully bed along flowlines for a given j-th runoff depth M_j. The partial gully depth D_ij at i-th point is

D_ij = z_0i − z_ij.

(18)

The partial longitudinal profiles differ substantially from each other. To calculate the most probable shape of the stable gully longitudinal profile, all discharges passing through a channel must be taken into account. The channel deformations during a given discharge are proportional to the magnitude of the sediment transport rate and duration of this discharge [20]. Following this assumption, the following sequence of calculations was proposed to calculate the most probable longitudinal profile of a stable gully:

(a)

The entire range of runoff depths is divided into N equal intervals, and probability p_k of each j-th interval of runoff depth is calculated.

(b)

For each j-th runoff depth M_j in the middle of the j-th interval, the partial longitudinal profiles with bottom altitudes z_ij are calculated at each i-th point on the flowlines with Equation (17).

(c)

For each j-th runoff depth M_j the partial magnitude of the sediment transport rate E_ij is calculated with Equation (4), written without erosion coefficient as

$$E_{ij}\sim \frac{M_j{F}_i}{W_{ij}}{S}_{0i}$$

(19)

where S₀ is the inclination of the initial slope. Flow width is calculated with Equation (8).

(d)

The most probable stable gully bottom altitudes Z_pi are calculated weighted with partial magnitudes of the sediment transport rate E_ij and runoff depth probability densities p_j

$$Z_{pi}=\frac{{{\displaystyle \sum }}_{j=1}^N{z}_{ij}{E}_{ij}{p}_j}{{{\displaystyle \sum }}_{j=1}^N{E}_{ij}{p}_j}$$

(20)

The result of these calculations is the most probable stable gully profile along each flowline. The calculation requires information about the initial relief of the gully catchment and about boundaries of all lithological units, including topsoil with vegetation in the form of the DTM (Digital Terrain Model); the probability density function of surface runoff depths, all parameters and coefficients used. As was already mentioned, these coefficients are empirical and local and must be estimated for a given territory.

The most probable stable gully profile can be close to one of the partial longitudinal profiles calculated with particular surface runoff. However, it did not coincide completely with any of them. It was obvious that flow velocities along the most probable profile for larger surface runoff values could be greater than the critical velocity of erosion initiation. The ultimate stable gully profile was, therefore, lower than the most probable one.

2.2. Area–Slope Approach

The product of the critical slope S and the power-law function of contributing catchment area F determine the condition of erosion initiation within this approach [7,8]:

$$a=S{F}^b$$

(21)

It is also can be written in the form:

$$S=a{F}^{-b}$$

(22)

In several works [9,21,22] Equation (21) is written in the form

$$a_m{S}^{-\propto }=C$$

(23)

where C is critical contributing area of the catchment. It was shown in [22] that C was similar to Horton’s [23] critical distance for channel initiation.

Equations (21)–(23) bear equal information, as α = 1/b and a_m = $a^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$b$}\right.}$.

Equation (21) is an empirical analogue of the formula for the critical shear stress [24] or the square of the critical flow velocity U_cr [25], at which erosion begins. The flow in gullies is usually of turbulent type. Therefore, flow velocity U was expressed in terms of the slope and catchment area using the Chezy–Manning formula (Equation (15)):

$$U^2=\frac{S}{n^2}{d}_{}^{4/3}=\frac{S}{n^2}{p}^{4/3}{\left(kMF\right)}_{}^{4m/3}$$

(24)

The parameter a in Formula (21) will be written as

$$a=\frac{{\left(nU\right)}^2}{p^{4/3}{\left(kM\right)}_{}^{4m/3}}$$

(25)

and the exponent b as 4m/3.

According to [6], the empirical information for a set of investigated regions showed that exponent b varied across a rather narrow range and its mean was about 0.4 (0.38). Therefore, m in Equation (24) must be 0.3, and this value fits the measurements in the gullies of the Yamal peninsula [26].

The coefficient a in Equation (21) was much more variable even within one catchment area [6] since it depends on many factors. These were the shape of the gully channel, the roughness of its bottom, the critical velocity of erosion initiation and the water runoff depth. The last characteristic was the most difficult to determine for calculation with Equation (21). The surface runoff depth varies through time (and less, in space) and is a probabilistic variable. The choice of the value (probability) of this quantity was largely subjective and at least should be justified in each specific case.

The way to avoid this uncertainty was to transform Equation (24) thus that the indefinite surface runoff depth is regarded not as a predictor but as a response:

$$M_{cr}=\frac{{\left(n{U}_{cr}\right)}^{3/2m}}{S_{}^{3/4m}{p}^{1/m}{\left(kF\right)}_{}^{}}$$

(26)

It is worth noting that Equation (26) is the version of Equation (23), where both parts are divided by the critical catchment area F (C in Equation (23)) and multiplied by the runoff depth:

$$M=\left(a_{m}M\right)S^{-\propto }{F}^{-1}$$

(27)

This version is novel and was not discussed in the literature. All parameters and variables at the right side of Equation (26) can be determined unequivocally for every point on a catchment. The calculated value of surface runoff depth is also unique for each point at the catchment and has the meaning of critical. If at a given catchment point with some slope and area, the actual runoff depth is greater than or equal to the critical value calculated by the Equation (26)

$$M\ge M_{cr}$$

(28)

the erosion initiation is potentially possible at this point. This critical runoff depth M_cr and risk of erosion initiation have determinable probability (duration) P_M. Equation (26) is used further in two following models of gully erosion potential estimation.

2.2.1. Total Erosion Potential (TEP) Model

The critical runoff depth M_cr is the minimal runoff depth, which initiates erosion. Therefore, erosion occurs during all flows with runoff depths more than M_cr. The influence of all flows with runoff depths more than M_cr must be taken into account to calculate the gully erosion potential. In the TEP model, the similar procedure, as in the STABGUL model, is proposed:

(a)

The entire range of flow depths is divided into N intervals, and probability p_j of each j-th interval is calculated.

(b)

For each j-th runoff depth M_j in the middle of the j-th interval, the partial erosion rate $E_{ij}$ is calculated at each point i-th on the catchment with Equation (4).

$$E_{ij}={\left(k_E\frac{Q_j}{W_j}S\right)}_i={\left[\frac{k_E}{pw}{\left(M_{j}F\right)}^{1-mw}S\right]}_i$$

(29)

(c)

Erosion potential $E_{pi}$ by flows with runoff depths more than critical is calculated as the sum of partial erosion rates $E_{ij}$, weighted with runoff depth probability p_j, beginning from the interval j_cr with critical runoff depth M_cr:

$$E_{pi}={{\displaystyle \sum }}_{j=jcr}^N{E}_{ij}{p}_j$$

(30)

(d)

To remove coefficients in Equation (29), relative erosion potential R_Ei is estimated by dividing E_pi by its maximum value E_Ni:

$$E_{TPi}=\frac{{{\displaystyle \sum }}_{j=jcr}^N{E}_{ki}{p}_j}{E_{Ni}}=\frac{{{\displaystyle \sum }}_{j=jcr}^N{\left(M_j\right)}^{1-mw}{p}_j}{{{\displaystyle \sum }}_{j=0}^N{\left(M_j\right)}^{1-mw}{p}_j}$$

(31)

The influence of morphometric characteristics S and F at the i-th point of the catchment is reached through the value of M_cr from Equation (26), which determines the number of intervals of M from jcr to N, used in the calculation of E_TPi.

The value of E_TPi varies from 1 (at M_cr = 0) to 0 (at M_cr ≥ M_max). At M_cr = 0 (j = 0), erosion begins at a given point for any, the smallest values of the runoff depth and E_TPi is maximum. Accordingly, at M_cr ≥ M_max (j = N), the runoff depth never exceeds the critical depth, and erosion will not occur at this point.

The input information to run TEP is simpler than for the previous models, as morphometry for each flowline is not required. Instead, the digital models for flow accumulation and slope are used, which can be calculated from DEM with most GIS. Empirical coefficients and parameters in Equation (26) are used in this model. Hydrological information is represented in the form of the probability density function of runoff depths, as in the STABGUL model.

The output is the values of the relative erosion potential E_TPi, calculated for each pixel on the catchment area DEM. The probability density function of E_TPi is usually negatively skewed [26], it is better to use the logarithms of E_TPi.

2.2.2. Express Estimation of Gully Erosion Risk

Probability density functions of runoff depth for a small catchment can often be approximated with an exponential relationship [27]. For this case, E_TPi in Equation (31) is an inverse power-law function of M_cr. This is the base to formulate the express model of gully erosion risk GER, (the term risk is used here instead of the term potential to distinguish these two models) when power-law functions are used instead of Equation (31). Using the exponent 2m/3 in Equation (26), the expression for such function of critical runoff depth takes the form:

$$M_{cr}^{2m/3}=\frac{n{U}_{cr}}{\sqrt{S} p^{2/3}{\left(kF\right)}_{}^{2m/3}}$$

(32)

Therefore, Equation (31) is transformed to the expression of gully erosion risk E_GR:

$$E_{GR}=1-{\left(\frac{M_{cr}}{M_{max}}\right)}^{2m/3}$$

(33)

As in Equation (31), the value of E_GR varies from 1 (at M_cr = 0) to 0 (at M_cr ≥ M_max), which corresponds to a change in the degree of potential erosion risk. At M_cr = 0, erosion begins at a given point for any, the smallest values of the runoff depth and E_GR is maximum. Accordingly, at M_cr ≥ M_max, the runoff depth does not exceed the critical depth, and erosion will not occur at this point.

The input information to run GER comes from digital models for flow accumulation and slope, which can be calculated from DEM in most of GIS. Empirical coefficients and parameters p, m, U_cr and n are used in this model. Hydrological information is restricted to the maximum daily runoff depth. The outputs are the values of the gully erosion risk E_GR, calculated for each pixel on the catchment area DEM.

3. Results

Since each of the presented models of gully erosion was characterized by a certain degree of simplification of the real process, it is rational to compare calculation results with measurements of actual erosion. Of greatest interest is the comparison of calculations based on the most complicated model GULTEM with actual data. If successful, the calculation results for simpler models can be compared with calculations from GULTEM.

3.1. GULTEM Model Validation and Calibration

The common definition of validation offered by the Society for Computer Simulation Technical Committee on Model Credibility [28] is “Substantiation that a computerized model within its domain of applicability possesses a satisfactory range of accuracy consistent with the intended application of the model”. In our case, “a satisfactory range of accuracy” depends on the quality of the measurements on the natural object and on the accuracy of the input data necessary for calculations. The catchment selected for the validation and calibration procedure (Figure 1) is situated on the Yamal Peninsula, West Siberia, close to the most severely gullied part of the peninsula. The investigated gully P-1 did not exist in 1986. After construction in 1986–1987 of the exploitation camp at the top of the catchment, erosion and thermoerosion began due to the increased water supply and vegetation cover damage. In 1988, the gully P-1 length was 450 m, in 1989—740 m, in 1990–1991—940 m and in 1995—965 m. The gully head reached the buildings of the exploitation camp, and repeated filling of the gully head with heavy loam by bulldozers stopped the gully growth. The observations showed that gully P-1 was still active in 2007 [18] and further (Figure 1), increased in depth and volume, though the gully length did not increase further. Therefore, the period 1991–1995 of not regulated gully activity was used for model validation and calibration with two longitudinal profiles of the gully bed measured in July 1991 and 1995 (Figure 2).

Empirical coefficients and parameters, required to run GULTEM, were estimated from morphometric and hydraulic measurements on selected and adjacent catchments: the coefficient of thermoerosion k_TE = 2.5 × 10⁻⁵ m s⁻¹ (°C)⁻¹; the gully walls inclination φ = 0.6 rad; bank erosion coefficient k_b was calculated via [18]; coefficients in Equation (8) are pw = 3.0 m (m³/s)^−mw and mw = 0.4; coefficients in Equation (9) are p = 0.21 m (m³/s)^−m and m = 0.3; Manning’s roughness coefficient n = 0.06 m⁻¹ s; critical velocity of erosion initiation U_cr = 0.2 m s⁻¹. Meteorological information, in form of changing through time air temperature and precipitation, was taken from ERA5 reanalysis. The sequences of surface runoff depths were calculated with the hydrological model, validated on the hydrological measurements on the selected catchment [27].

The combined erosion coefficient k_E in Equation (4) varied in calibration procedure in the range 6.5 × 10⁻⁵–6.5 × 10⁻⁴ m⁻¹. The difference between calculated (Z95c) and measured (Z95m) gully bed altitudes in 1995 was estimated with

$$RMSE=\frac{\sqrt{{{\displaystyle \sum }}_1^N{\left[\left(Z95c-Z95m\right)dx\right]}^2}}{L}$$

(34)

where N is the number of gully segments with the measured altitudes, dx—the distance between segments, L—gully length in the year 1995. The best fit of Z95c and Z95m was at k_E = 3.6 × 10⁻⁴ with RMSE = 0.057 m (Figure 3). This value shows a consistency between the model results and measurements and “a satisfactory range of accuracy”.

3.2. The Comparison of GULTEM and STABGUL Models

The calculations with GULTEM and STABGUL models were performed for the gully P-1. The same empirical coefficients were used in both models. The probability density function for surface runoff, used in STABGUL (Figure 4), was derived from a 30-year long sequence of surface runoff depths [27], which was used in GULTEM calculations. The GULTEM model was run for a 1200-year period (40 repetitions of a 30-year long sequence) to obtain a stable longitudinal profile, suited for comparison with the stable gully calculated from the STABGUL model.

The results of calculations with GULTEM in Figure 5 show a gradual increase of the gully depths and the length of the stable segment of the gully bed through time. At these segments of the gully, its depth does not change anymore. After about 1000 years in calculations, almost the entire gully bed becomes stable, and flow velocities for any discharge become less than critical at all points.

Calculations with STABGUL with the above-described algorithm resulted in the most probable stable longitudinal profile with depths (line 3 in Figure 5), close to the depths from GULTEM after about 120–150 calculation years. This period is close to the period of stabilization of most human-induced gullies on the East European plain [29]. These depths are about 40% less than from GULTEM after 1000 years of calculated erosion (line 2 at Figure 5). The ultimate stable longitudinal profile calculated with STABGUL for the maximum runoff depth (line 4 at Figure 5) fits well with the stable profile after 1000 years from GULTEM. This result shows that the most probable stable longitudinal profile from the STABGUL model can accurately simulate a stable gully after 120–150 years of its lifetime, which can further slowly be transformed to its ultimate stability. It also points out the significance of high discharges on the process of gully erosion.

3.3. The Comparison of the GULTEM and TEP Models

The total erosion potential was calculated in the TEP model along the same flowline and with the same empirical coefficients and hydrological data as gully depth in GULTEM. The similarities and differences in the results (Figure 6) are explained by the similarities and differences in the models’ formulations. Both models operate with the rate of erosion E (Equation (4)), which depends on local slope inclination and specific discharge.

At a given point, the gully depth in GULTEM and erosion potential in TEP are the sum of the product of the erosion rate values and their durations. In GULTEM the products are included to the sum if flow velocity is more than critical, in the TEP model, if runoff depth is more than critical. Therefore, we can expect a correlation of the calculation results. The differences shown are due to the significant simplification in the TEP model of the erosion process. The result of calculations in the TEP model shows the potential of erosion, its position at the catchment, and potential intensity, not gully geometry. In the TEP model, the influence of the altitude of the basis of erosion at the gully mouth and the processes of gully walls slumping are not taken into account. The main difference is the absence in the TEP of the gully longitudinal profile self-forming evolution. Slope inclination in the TEP model is the inclination of the initial slope. Therefore, the correlation between the results of calculation with the two models quickly decreases with time (Figure 7).

3.4. The Comparison of TEP with GER Model

As in previous cases, gully erosion risk was calculated in the GER model along the same flowline and with the same empirical coefficients and hydrological data. There is a close relationship (power-law or logarithmic) between calculation results from the TEP model of total erosion potential E_TP and of gully erosion risk E_GR in GER (Figure 8). Therefore, the comparison of calculation results from the GER model and GULTEM is more or less the same as in the previous paragraph.

4. Discussion

The problems to be discussed are the availability to possible users of input information, required for calculations with detailed models, and the limitations in simplified models.

The GULTEM is the only model of the four described in this article in which gully evolution can be simulated in “real” time and gully geometry change though time can be adequately simulated. In this model, it is possible to take into account the temporal trends in hydrological information, common in the conditions of global climate change. The temporal and spatial changes in land use, for example, properties of vegetation cover, are also possible to include in the calculations. The large number of input empirical coefficients, morphometric and hydrological data required to run GULTEM is a limitation of the model. The calculations require information for all or particular flowlines in the form of digitized 2D longitudinal profiles of the altitudes of top catchment surface Z (see Figure 2) and of the surfaces Z_j of other underlying lithological units along these flowlines. The contributing catchment areas F and coordinates x, y are also collected for the same points. If a digital elevation model (DEM) is available for the catchment, most of this information can be collected using standard tools included in any geographical information system (GIS). The lithological composition of the territory may require investigations in the field.

Parameters and coefficients, used in Equations (4–13), are also better to obtain via measurements at the selected territory, as in [18]. Values of some of these coefficients can be found in special literature. The coefficients in regime Equations (8) and (9) for gullies are discussed in [6,13], the coefficient of thermoerosion in [18], the critical velocity of erosion initiation in [19,30], Manning’s roughness coefficient in [31,32]. The gully walls inclination φ is assumed to equal the angle of repose for a given lithology (see discussion in [33]).

Experimental data on the values of the erosion coefficient k_E for cohesive soils, common for gullies, is limited. Experiments and model calibration were provided on soils of the Yamal peninsula [17,18], for silt loams in the gullies of the coast of George Lake in Australia [34], loess soils in the Manawatu River basin and rendzina in the Waimakariri River basin, New Zealand [35] and granodiorite saprolites at the basin of the Mbuluzi River (Swaziland) [36]. The variability of this coefficient is very high. In the direct measurements of sediment budget in the gullies of the Yamal peninsula during the snow thaw and summer rains the value of k_E changes from 0.0008 to 0.01 for loams and silt loams with cohesion from 10 to 40 KPa, estimated with a torvane shear tester, and up to 1.3 for loose silt sand. The measurements in Australia showed a similar range of k_E: for the loams with cohesion 30–40 KPa 0.0028–0.0034, with cohesion 50–70 kPa—0.0008. The measurements in New Zealand for the loess soil with cohesion 51 kPa showed k_E = 0.004–0.005 for soil with organic carbon content 2.4% and k_E = 0.021 for the same soil with organic carbon content 0.44%. For rendzina with cohesion 22 kPa and 3.5% content of organic carbon, k_E was 2 × 10⁻⁵. Model calibration for the gullies in Swaziland showed k_E value 0.006 for sandy saprolites with cohesion 4.5–9 kPa. The listed information was not enough to find the relationship between k_E and eroded cohesive soil properties, and the large range of its values does not permit the recommendation of values for calculations on not-investigated objects. The only way to use GULTEM for gully erosion prediction is to calibrate the model on the measurements (see Section 3.1). Usually, such calibration is possible in the framework of detailed projects of land management [36,37].

The complexity in erosion coefficient estimations is the reason to use simpler models for preliminary applications. The models STABGUL, TEP and GER do not require the erosion coefficient. Other parameters and coefficients used in these models, if not estimated in the field measurements, can be found in special literature.

Probability density functions (PDF) of surface runoff values are used in STABGUL and TEP models but for different erosion form lifetimes. The minimum prediction duration in the STABGUL model is about 150 years. Therefore, it is assumed that the time series of hydrological characteristics are stationary. This is not true in the current global climate change, and predictions with STABGUL will not include the trends in the time series of surface runoff. Keeping in mind this limitation, STABGUL is a powerful tool to find the longitudinal profiles, lengths, depths, widths, areas and volumes of the gullies close to the stage of stabilization [16]. This can be useful for the classification of the gullies in the intensively gullied territory to separate stable gullies and those that are still active [38].

The total erosion potential model is based on the modification of the area–slope approach in gully erosion investigations when critical runoff depth is calculated. The TEP model estimates possible gully position in the catchment and erosion intensity. This information can be enough in preliminary projects of land management. The absence of the need for the flowline geometry is the significant simplification: the calculations are performed directly for given points in the catchment, for which slope and contributing area are measured. Calculations of the total erosion potential with the TEP model are more flexible in the conditions of trends in hydrological time series. Much smaller prediction periods are needed in the TEP model (one or a few decades), allowing it to use step functions of PDFs to take into account temporal changes in hydrological statistics. The main limitation of the TEP model is the absence of prediction of the gully geometry.

The GER model requires very limited hydrological data—only maximum daily surface runoff depth for the prediction period. At the same time, the GER model output for a certain catchment contains nearly the same information as total erosion potential from the TEP model, which uses more hydrological data. This similarity is explained by the more or less close relationship between parameters of PDFs of surface runoff values and their maximums. Such relationships are different for different PDFs, as other parameters of PDFs must be taken into account. Therefore, transfer from TEP to GER output will also be different for different hydrological conditions, as unnamed parameters of PDFs are implicitly used in such transfer.

5. Conclusions

The four models described here encompass all now existing approaches to gully erosion prediction. Comparison of the results of calculations for the same gully using these four models makes it possible to choose the model that is most suitable for use in a land management project depending on its objectives and the available information. The set of gully erosion models discussed here can meet the requirements of land management projects with different levels of detail.

Use of the dynamic gully erosion model GULTEM is recommended for calculations of gully geometry transformation in time and space for the most detailed projects of land management. Flowline geometry, soil texture, parameters and coefficients, used in GULTEM, can be obtained via measurements at the selected territory and from DEMs. The calculations take into account the temporal trends in hydrological information, important for conditions of global climate change. The main difficulty in modelling is the estimation of the erosion coefficient, which usually requires calibrating the model.

The input data for the stable gully model STABGUL is the same as for GULTEM, except for the erosion coefficient, which is the most difficult to estimate. Therefore, STABGUL is easier to apply. The STABGUL model is useful for projects where only final gully geometry is required. Only stationary hydrological time series are used in the model.

Possible gully position in the catchment and erosion intensity are estimated with the total erosion potential (TEP) and gully erosion risk (GER) models. The novel type of area–slope approach allows effective use of hydrological information in the form of critical surface runoff for erosion initiation. These models can be useful for preliminary estimates of the possible effects of different land uses.

The presented models, with proper selection of initial data and necessary calibration, can be used for territories with different climates, soils, and land use. They are most effective for areas with possible linear disturbance of the natural vegetation cover, which is typical for the use of land for pastures, for the initial stages of urbanization and the industrial development of new territories.