The Theoretical Approach to the Modelling of Gully Erosion in Cohesive Soil

The stochastic gully erosion model (STOGEM) is based on a combination of deterministic mechanics and a stochastic description of the erosion control factors. The main proposition in the model is that the depth of the active surface layer of eroded cohesive soil is equal to one particle diameter, and the deposition of eroded particles is negligible. The erosion rate at the gully bed is calculated directly from the equation of the balance between driving and resistance forces acting on soil particles in flowing water using the probability density functions of stochastic variables: flow velocity, soil aggregate size and cohesion. Probability density functions of cohesion in the model vary through time and space during the erosion event due to the changes in soil composition—armoring and loosening. This theory is still far from achieving practical application, but opens up a new way for better understanding the experimental results of gully erosion and shows the direction for future investigations.


Introduction
Gullying is a form of linear erosion of loose and cohesive soils (and of rocks) by a concentrated water flow [1]. It is possible to classify the forms of linear erosion by their size and stage of development. The sequence of such linear forms begins with the smallest ones-the rills on the slopes. With increasing depth, the larger linear forms are ephemeral gullies [2], which are usually destroyed by plowing, followed by typical gullies. This quantitative difference in the depth causes the main qualitative difference between active typical gullies and smaller linear erosion forms-the instability of gully walls. The processes of gully bed incision and bank erosion cause the increase in gully wall inclination and soil slumping and falling of different types. These processes are typical for the first stage of gully development. Gully formation is very intensive and its geometry (length, depth, width, area, volume) is far from stable and changes rapidly [3]. The first stage lasts from 4 to 10 years in loose sands and frozen loams, and up to 100-150 years in typical agricultural landscapes [4]. The gullies cut slopes to their whole potential (possible) length during this period, and their walls eventually stabilize.
After gully wall stabilization, usually due to vegetation, gully development moves into the second stage [3]. The gully width increases slowly due to slow soil movement (creep). Slumping can still occur at the gully head, and so gully length can also slowly increase through gravitation processes and suffusion, not by fluvial erosion [5]. At this stage, the gully is transformed into a dry valley. Often sediments washed from the catchment are deposited at the dry valley's bottom. The walls of a dry valley become more and more gentle, its depth decreases, and the dry valley is transformed to a shallow linear hollow on the slope. The formation of hollows due to sediment deposition in dry valleys takes a significant amount of time and is discovered in paleo-geographical investigations [6]. Gullies can be also classified by their position: gullies on the hillslopes, bank gullies on the valley slopes, and bottom gullies at the valley bed. Gullies can be discontinuous and continuous, differing by soils and type of vegetation on gully slopes [7].
Mathematical models of different levels of simplification [7][8][9][10] describe all of the abovementioned forms, stages and processes. All of these models represent significant empirical components in process description. The goal of the current paper is to formulate the general theoretical principles of only one, but the main, process of gully development in cohesive soil-erosion by water. The main approach is a stochastic description of the synergy of the main destructive and stabilizing forces, which causes soil particle detachment.

Gully Erosion by Water
The rate of gully incision is controlled by water flow velocity, depth, turbulence, temperature, as well as soil texture, soil mechanical pattern, and the level of protection by vegetation. These characteristics are described by equations of mass conservation: where Q s = Q C is the sediment transport rate (m 3 s −1 ), Q is the water discharge (m 3 s −1 ); X is the longitudinal co-ordinate (m); C is the mean volumetric sediment concentration; C w is the sediment concentration of the lateral input; q w is the specific lateral discharge (m 2 s −1 ); E is the erosion rate or the mean particle detachment rate (m s −1 ); E b -bank erosion rate (m s −1 ); W is the flow width (m); d is the bank height (m); and U f is the sediment particles' falling velocity in the turbulent flow (m s −1 ). During the first stage of gully development in cohesive soil, the accumulation on the gully bed of soil particles, eroded inside the gully or from the lateral input, is usually negligible. Flow velocity of erosion initiation in cohesive soils is higher than the critical flow velocity when the settling of detached particles occurs. The rate of bank erosion E b is related to the rate of erosion E of the gully bed [10]. Therefore, the main term in Equation (1) is the erosion rate E.
The soil erosion rate E is the mean rate of the lowering of the soil surface Z 0 or the increase in volume ∆V eroded from the area S during the time interval T.
The erosion volume ∆V is the sum of the volumes of individual eroded soil particles (aggregates) V ai with correction on soil porosity ε: After multiplying and dividing each component of the sum (3) by the product of the soil particle bottom area s i and the duration of particle detachment τ i , Formula (2) takes the form: where the term V ai /s i τ i represents the instant and local rate of particle detachment α, and the term s i τ i /ST represents the probability density of this particle detachment. Therefore [11]: where p α is the spatial-temporal probability density function (PDF) for α. Equations (4) and (5) show that the mean soil erosion rate is the result of spatial (on the area S) and temporal Earth 2022, 3

230
(during the period T) averaging of instant and local detachment rates α. The detachment rate α possess all positive values, including zero ("detachment" of stable particles).

Instant and Local Rate of Soil Particle Detachment
Erosion of cohesive soil is a rather slow process, even in gullies, as water mainly only comes into contact with, and thus only affects, the surface layer. The main proposition in the further consideration in this paper is that the depth of the active surface layer equals one particle diameter. The instant and local detachment is a discrete process: the surface becomes lower when and where a soil fragment is detached. A soil fragment becomes unstable, and its movement begins when and where the driving forces F dr exceed the maximum of the resisting forces F res and the resultant force Θ is more than zero [12]: The acceleration of such an unstable fragment within the soil surface layer in the direction of the resultant force Θ (along the axis Z) is described by Newton's second law: where ρ s is soil density. A soil fragment becomes detached when all links between particles are broken, i.e., when the soil fragment is removed from its initial position by a distance greater than its height D. The duration of detachment equals.
The local and instant (average for the period τ) detachment rate is Equations (8) and (9) are valid for the simplest case when Θ is independent of Z.
Equations (8) and (9) describe the behavior of a single soil particle. It is evident that is impossible in general to predict the exact force balance or to foresee the geometry and mass of a soil particle or aggregate at a given point at a given time. The main variables in these equations are of a stochastic nature. Therefore, it is possible to describe these variables with probability density functions and use these functions to calculate the mean rate of detachment (erosion rate) with Equations (5) and (9) [13].

Probability Density Function for the Rate of Detachment
The detachment rate calculated with Equation (9) is the function of multiple factors: soil particle geometry and density, as well as driving and resistance forces. Each of these factors (often compound) is a stochastic variable characterized by a probability density function (PDF). To obtain the PDF for the function of random variables (in our case, for the particle detachment rate), it is necessary to use the appropriate calculation techniques. The main technique is the formula for calculating the mean E(Y) of a function of random variables Y = g(x 1 , x 2 , . . . x n ) [14]: For independent variables, the multivariable PDF p(x 1 , . . . x n ) is equal to the product of the PDFs of each of the random variables.

The Main Resistance and Driving Forces
The process of a single particle's detachment from the surface of cohesive soil caused by water was described in detail by Mirtskhulava [15]. We shall follow this work with additions.
The main resistance forces are: submerged weight of soil fragments, with slope inclination γ taken into account, projected in the direction of driving forces β: and the geo-mechanical force of cohesion: The latter is a combination of electro-chemical, capillary and friction forces, which is parameterized by the product of cohesion C (in terms of Coulomb law) and contact area between and within soil aggregates S b .
Hydrodynamic driving forces are the gradients of static and dynamic pressure, acting on the soil surface, depending on geometry of the soil surface, the shape and size of soil aggregates, the flow depth and the velocity distribution. The field of these gradients can be calculated with 3D hydrodynamic models. More often [15], these pressure gradients are calculated with the use of flow velocity and aggregates geometry as drag (F d ) and lift (F l ) forces, which are parameterizations of longitudinal and vertical pressure gradients: where C R is the coefficient of drag resistance; C y is the coefficient of uplift; U is the instant near-bed flow velocity; ρ s and ρ are the soil aggregate density and water density, respectively; S d is the area of aggregate cross-section exposed perpendicular to flow; and S a is the cross-section area of the soil aggregate parallel to the flow (vertical projection). The direction of the sum of drag and lift forces (the angle with the mean soil surface) is: With this parametrization, the local detachment rate is [13] (with modification): 2ρ s cos γ sin β and k C = DS b 2ρ s V a . If the expression under the square root is equal to zero or negative, then α = 0. For this case, Equation (16) is a generalized form (with cohesion) of a well-known expression for incipient motion criteria (see, for example, [16]).
The rate of detachment is the stochastic variable, which is the function of six other independent stochastic variables: flow velocity, particle (aggregate) size, soil cohesion and three geometric-kinematical coefficients k U , k D and k C . In the case of turbulent flow, C y , C R and k U are constants [15]. If we assume a simple shape and composition of soil particles (aggregates), then the coefficient k D is also constant. Then, Equation (16) contains three independent stochastic variables: flow velocity U, soil particle size D and force of cohesion k C C, and Equation (10) takes the form: where p( √ k U U), p(k D D) and p(k C C) are probability density functions for flow velocity, soil particles vertical size and cohesion, respectively. In numerical calculations, the continuous quantities are replaced by discrete one, infinitesimal by finite, and integrals by sums:

Probability Density Function for the Factors of Soil Erosion
When geometric-kinematical coefficients and PDFs for flow velocity, soil aggregate size and cohesion are known from the theory or experiment, the PDF for instant and local erosion rate α and, therefore, the mean rate of erosion can be calculated with Equation (17). Various combinations of the input mean values of the factors, as well as different PDFs, lead to a great variability of the resulting relationships between erosion and controlling factors.
Measurements in different environments show [7,10,17] that the flows in the gullies are usually shallow, turbulent and often supercritical. A large number of measurements of turbulent characteristics using different techniques in such shallow flows [18,19] showed Gaussian PDF of fluctuations of instant longitudinal velocity U: where U mean is the mean velocity and σ is the standard deviation of velocity fluctuations. For small mean velocities, the function p(U) can be asymmetrical. The ratio of the standard deviation of velocities and kinematic velocity U * increases from the flow surface to the bottom [19]; where H is the flow depth and z equals the roughness height. According to Mirtskhulava [15], in shallow turbulent flow C y = 0.1 and C R = 0.42; then k U = (0.1S a + 0.42S d ).
The PDF of the aggregate diameters can be expressed as a lognormal distribution [20,21]: The parameters in Equation (20), calculated with arithmetical values, namely the mean of soil particle height D m and standard deviation σ D , are The direction of sum of drag and lift forces is β = arctan 0.238 S a S d , therefore The details of the influence of soil cohesion on erosion rate are not properly investigated. Mean values, measured with existing techniques, such as tore-vane and penetrometers of different kinds, are not well correlated with aggregate stability tests and measured rates of soil erosion (see the results of experiments in [22]). Even less is known about PDFs of aggregates cohesion. Mirtskhulava [15] determined a normal distribution function for the variability of soil cohesion, measured by a round stamp. The same PDF fit to the textural tensile strength of one hundred 2-3 mm rounded aggregates [23], but the initial data, published there, show a distinct asymmetrical distribution. Our experiments, with the cutting of the surface of soil sample with a blade [24], showed the lognormal distribution of soil strength (Equation (20)). We shall use lognormal distribution in further calculations, keeping in mind that this distribution changes over time during an erosion event.

The Algorithm of Erosion Rate Calculation
Two opposite processes transform the initial spatial distribution of soil cohesion: soil armoring and soil loosening ( Figure 1). The armoring due to the removal of unstable soil particles and aggregates leads to the following transformation of cohesion PDF through time: Earth 2022, 3, FOR PEER REVIEW 6 of aggregates cohesion. Mirtskhulava [15] determined a normal distribution function for the variability of soil cohesion, measured by a round stamp. The same PDF fit to the textural tensile strength of one hundred 2-3 mm rounded aggregates [23], but the initial data, published there, show a distinct asymmetrical distribution. Our experiments, with the cutting of the surface of soil sample with a blade [24], showed the lognormal distribution of soil strength (Equation 20). We shall use lognormal distribution in further calculations, keeping in mind that this distribution changes over time during an erosion event.

The Algorithm of Erosion Rate Calculation
Two opposite processes transform the initial spatial distribution of soil cohesion: soil armoring and soil loosening ( Figure 1). The armoring due to the removal of unstable soil particles and aggregates leads to the following transformation of cohesion PDF through time:

1.
The probability density p(k C C) in the part of the cohesion PDF where resistance forces are less than driving forces decreases due to the erosion E iC of soil with particular cohesion C iC The initial PDF transforms into intermediate PDF (p*) 2.
Simultaneously, the intermediate PDF of cohesion is transformed due to the exposition of fresh initial soil in the "windows" of the eroded surface layer to PDF of armored soil (p a ) where the index 'iC' indicates a particular soil cohesion and related erosion rate, the index 'í' indicates a sequence in time, and p 0 (k C C) is the PDF of cohesion in the initial soil. This process leads to the increase in the proportion of the surface covered by stable soil aggregates, and an increase in the mean soil cohesion. The rate of erosion decreases through time when armoring is the prevailing process. Soil loosening takes place due to the destabilization of stable aggregates when links between and/or within soil aggregates are weakening. This is complicated process, including decrease of electro-chemical, capillary and friction forces in saturated soil, the removal by flow of unstable soil particles and aggregates, which previously stabilized other aggregates, the widening of pores and cracks due to changes in water pressure in soil. Destabilization also occurs due to the vibration of stable aggregates caused by turbulent flow [15].
The failure of links decreases the contact surface area S b between and within soil aggregates over time. The rate of failure of links in soils can be described with the common exponential failure function [25].
The parameter λ > 0 must be estimated from experiments. The remaining contact area at time T i+1 is the product of the contact area at time T i and the failure function (Equation (26)) Loosening leads to an increase in the proportion of the surface covered by unstable soil aggregates, and leads to a decrease in the mean soil cohesion. The rate of erosion increases through time when loosening occurs at a higher rate than soil armoring. Soil loosening is the main cause of erosion at the final stage of the erosion event, when the soil surface is nearly completely stabilized by armoring.
The final PDF of cohesion after the armoring and loosening cycle, as the PDF of function (27) [14], is equal to: The rate of gully erosion at a given moment is calculated from the equation of the balance between the driving and resistance forces acting on a soil particle in flowing water (Equation (17)) after soil armoring and after soil loosening (Figure 1). The armoring-loosening cycle exists only in the numerical representation of the process of erosion. In reality, soil armoring and loosening is a unified and simultaneous process.

The Materials for Comparison of Calculations with Measurements
The calculations with any model must be compared with the measurements for possible validation and calibration. Field experiments were performed in the low-altitude hills of Ballantrae Hill Country Research Station near Woodville, New Zealand ( Figure 2). The local soil belongs to Pallic Soils of the Wainui series formed on loess [26]. The entire depth of the A + AB-horizon is 40-50 cm. The organic carbon content in the topsoil is 2.4%, and 0.4% in the parent loess. The topsoil is well-structured and highly water-stable: wet sieving of 2-4 mm aggregates for 30 min showed that 88% was retained for this class. The aggregate stability in the parent loess is much lower-the proportion of aggregates retained on the sieve was 7.6%. The mean size of soil aggregates obtained by dry sieving D s ≈ 4S a π was 1.83 mm, with a standard deviation σ D = 0.88 mm. The aggregates were flattened, to a plate or ellipsoidal shape, with D m D s ≈ 1/3. The mean soil strength, measured with a tor-vane, was nearly equal for parent loess and topsoil: 52 and 51 kPa, respectively. The soil density ρ s was 1460 and 1230 kg/m 3 , respectively.

The Materials for Comparison of Calculations with Measurements
The calculations with any model must be compared with the measurements for possible validation and calibration. Field experiments were performed in the low-altitude hills of Ballantrae Hill Country Research Station near Woodville, New Zealand ( Figure 2). The local soil belongs to Pallic Soils of the Wainui series formed on loess [26]. The entire depth of the A + AB-horizon is 40-50 cm. The organic carbon content in the topsoil is 2.4%, and 0.4% in the parent loess. The topsoil is well-structured and highly water-stable: wet sieving of 2-4 mm aggregates for 30 min showed that 88% was retained for this class. The aggregate stability in the parent loess is much lower-the proportion of aggregates retained on the sieve was 7.6%. The mean size of soil aggregates obtained by dry sieving  Two erosion plots were organized on two steep slopes (at straight parts), which accurately represent the gully beds. The abandoned road on the parent loess was used to organize plot 1, at 17.4 m long and 0.55 m wide with an inclination of 0.235 (13.2°). Plot 2 was cut into the topsoil by the removal of the upper 5-7 cm of soil with grass cover to avoid grass roots' influence on the erosion rate. The residual content of thin roots in the soil was less than 0.1% by weight (the initial content in the upper 5-7 cm was about 4%). Plot 2 was 11.5 m long, 0.32 cm wide and had an inclination of 0.483 (25.8°).
Experiments were conducted during November 2001. Five runs, each lasting half an hour, were performed at each plot with discharges of 1.48-11.2 L/s controlled with a Vnotch weir (Table 1). Flow width was measured at 8-10 cross-sections and the mean flow velocity was measured using salt injections. The mean flow depth was calculated from discharge, flow width and velocity. All flows were turbulent and supercritical. Samples of water with sediment particles, with a total volume of about 4 L, were taken at the end of the plot and at the end of the half-an-hour run. The sediment concentration was estimated by means of the standard procedure of filtering, drying, weighting and dividing based on water volume. Two erosion plots were organized on two steep slopes (at straight parts), which accurately represent the gully beds. The abandoned road on the parent loess was used to organize plot 1, at 17.4 m long and 0.55 m wide with an inclination of 0.235 (13.2 • ). Plot 2 was cut into the topsoil by the removal of the upper 5-7 cm of soil with grass cover to avoid grass roots' influence on the erosion rate. The residual content of thin roots in the soil was less than 0.1% by weight (the initial content in the upper 5-7 cm was about 4%). Plot 2 was 11.5 m long, 0.32 cm wide and had an inclination of 0.483 (25.8 • ).
Experiments were conducted during November 2001. Five runs, each lasting half an hour, were performed at each plot with discharges of 1.48-11.2 L/s controlled with a V-notch weir (Table 1). Flow width was measured at 8-10 cross-sections and the mean flow velocity was measured using salt injections. The mean flow depth was calculated from discharge, flow width and velocity. All flows were turbulent and supercritical. Samples of water with sediment particles, with a total volume of about 4 L, were taken at the end of the plot and at the end of the half-an-hour run. The sediment concentration was estimated by means of the standard procedure of filtering, drying, weighting and dividing based on water volume. Key: Q-discharge, U-mean velocity, W-flow width, d-flow depth, U * -kinematic velocity, σ U -velocity fluctuation standard deviation (from Equation (19)), E-erosion rate, Re-Reynolds number, Fr-Froude number.
The coefficient k U for ellipsoidal aggregates is 3ρ 16ρ s 0.1 + 0.42 S d S a and varied in the range 0.026-0.04 for the first plot and 0.031-0.047 for the second (depending on S a /S d ratio in the range 2-4). The standard deviation of velocity fluctuations was σ U ∼ = 2.1U * (according to [19]).
was 2.04 × 10 −4 for the first plot and 1.54 × 10 −4 for the second. The variation coefficient for cohesion was estimated as 0.59 in the experiments with the cutting of the surface of the soil sample with a moving blade [24].
The relationship between the erosion rate and the flow velocity shows a higher erosion rate at plot 1 than at plot 2 for the same velocity ( Figure 3). This is explained by the higher soil aggregate stability and organic matter content in the soil at plot 2. The power-law functions describe these relationships well. The exponents are different for these two plots (5.2 for the first and 4.1 for the second) and differ from the function commonly used in the majority of shear-stress erosion models (see, for example, [27]).   (19)), E-erosion rate, Re-Reynolds number, Fr-Froude number.
The relationship between the erosion rate and the flow velocity shows a higher erosion rate at plot 1 than at plot 2 for the same velocity ( Figure 3). This is explained by the higher soil aggregate stability and organic matter content in the soil at plot 2. The powerlaw functions describe these relationships well. The exponents are different for these two plots (5.2 for the first and 4.1 for the second) and differ from the function commonly used in the majority of shear-stress erosion models (see, for example, [27]).

General Numerical Experiments
The first set of calculations with the algorithm described above (Figure 1) show the erosion rate change through time (Figure 4). In the numerical experiments for different flow velocities, soil was "eroded" for 1800 or 3600 s. During this period, the mean cohesion of the soil surface layer and the instant erosion rate changed through time due to the action of soil armoring and soil loosening. The temporal trend changed its sign depending on the sign of the cumulative effect of armoring and loosening. The most common scenario was a general decrease in soil mean cohesion during the first 500 s, then a slight increase and stabilization after about 1000 s in the experiment. The erosion rate decreased during the first 5-20 s, then significantly increased and stabilized after about 1000 s in the experiment. The temporal changes in the percentage of the area of erodible soil P E , where the driving forces were greater than the resistance forces, corresponded with the trends in erosion rate. The influence of the soil particle size on the erosion rate was much smaller than that of cohesion, and is thus not discussed further.

General Numerical Experiments
The first set of calculations with the algorithm described above (Figure 1) show the erosion rate change through time (Figure 4). In the numerical experiments for different flow velocities, soil was "eroded" for 1800 or 3600 s. During this period, the mean cohesion of the soil surface layer and the instant erosion rate changed through time due to the action of soil armoring and soil loosening. The temporal trend changed its sign depending on the sign of the cumulative effect of armoring and loosening. The most common scenario was a general decrease in soil mean cohesion during the first 500 s, then a slight increase and stabilization after about 1000 s in the experiment. The erosion rate decreased during the first 5-20 s, then significantly increased and stabilized after about 1000 s in the experiment. The temporal changes in the percentage of the area of erodible soil PE, where the driving forces were greater than the resistance forces, corresponded with the trends in erosion rate. The influence of the soil particle size on the erosion rate was much smaller than that of cohesion, and is thus not discussed further.  The temporal trends in the mean cohesion of the surface layer and in the erosion rate are controlled by changes in cohesion PDFs ( Figure 5). A very small part of the PDF of the cohesion force at the segment, overlapping the PDF of the driving force ( Figure 5A), controls the rate of erosion and the percentage of erodible soil surface P E ( Figure 4C). The values of probability density (PD) are very small here ( Figure 5B). Therefore, numerical experiments focusing on this segment require precise calculations. The temporal trends in the mean cohesion of the surface layer and in the erosion rate are controlled by changes in cohesion PDFs ( Figure 5). A very small part of the PDF of the cohesion force at the segment, overlapping the PDF of the driving force ( Figure 5A), controls the rate of erosion and the percentage of erodible soil surface PE ( Figure 4C). The values of probability density (PD) are very small here ( Figure 5B). Therefore, numerical experiments focusing on this segment require precise calculations. The PD at the segment around the mode of PDF increases through time and the mode position shifts to lower values of cohesion due to the loosening of soil at the surface (Figure 5C). The PD at the main part of PDF to the right of the mode mostly decreases until ~500 s, and then stabilizes ( Figure 5D). The part of the PDF with a soil cohesion force greater than the driving force controls the temporal changes in the mean cohesion of the soil surface layer.

The Comparison of Calculated Erosion Rates with the Measured
In the situation where one of the main factors of erosion, i.e., the resisting forces of soil cohesion, are indefinite, it is possible to find these indefinite characteristics by comparing the measured rates of erosion and calculating them with Equations (17)- (28).
Numerical experiments were performed with input characteristics from Table 1. The fields of erosion rate E at given flow velocities were calculated for a variety of initial soil cohesions kCC and parameters λ in the failure function. Examples of such fields, represented by isolines of erosion rate E, are shown in Figure 6 for two runs in plot 1. The soil The PD at the segment around the mode of PDF increases through time and the mode position shifts to lower values of cohesion due to the loosening of soil at the surface ( Figure 5C). The PD at the main part of PDF to the right of the mode mostly decreases until~500 s, and then stabilizes ( Figure 5D). The part of the PDF with a soil cohesion force greater than the driving force controls the temporal changes in the mean cohesion of the soil surface layer.

The Comparison of Calculated Erosion Rates with the Measured
In the situation where one of the main factors of erosion, i.e., the resisting forces of soil cohesion, are indefinite, it is possible to find these indefinite characteristics by comparing the measured rates of erosion and calculating them with Equations (17)- (28).
Numerical experiments were performed with input characteristics from Table 1. The fields of erosion rate E at given flow velocities were calculated for a variety of initial soil cohesions k C C and parameters λ in the failure function. Examples of such fields, represented by isolines of erosion rate E, are shown in Figure 6 for two runs in plot 1. The soil erosion rate is characterized by equifinality. The calculation for each pair of k C C and λ at a given flow velocity leads to a unique value of the resulting erosion rate E. At the same time, calculations for different pairs of initial soil cohesion and parameter λ at a given flow velocity lead to the same resulting erosion rate, forming the isoline in Figure 6. The same measured rate of erosion at a given flow velocity, shown by bold dashed isolines, can also Earth 2022, 3 239 appear with different combinations of soil properties. This set of numerical experiments shows that the result of a single run with measured characteristics of flow and erosion rate cannot be used to find indefinite soil properties (initial soil cohesions k C C and parameters λ), which are used in the proposed model. erosion rate is characterized by equifinality. The calculation for each pair of kCC and λ at a given flow velocity leads to a unique value of the resulting erosion rate E. At the same time, calculations for different pairs of initial soil cohesion and parameter λ at a given flow velocity lead to the same resulting erosion rate, forming the isoline in Figure 6. The same measured rate of erosion at a given flow velocity, shown by bold dashed isolines, can also appear with different combinations of soil properties. This set of numerical experiments shows that the result of a single run with measured characteristics of flow and erosion rate cannot be used to find indefinite soil properties (initial soil cohesions kCC and parameters λ), which are used in the proposed model. The shape of the erosion rate isolines is different for different flow velocities. Therefore, a combination of such isolines can probably narrow the region of indefinite soil properties.
Keeping the initial soil cohesion constant in the numerical experiments for a given flow velocity, it is possible to minimize the difference between measured erosion rates and calculated ones with Equations (24)  The shape of the erosion rate isolines is different for different flow velocities. Therefore, a combination of such isolines can probably narrow the region of indefinite soil properties.
Keeping the initial soil cohesion constant in the numerical experiments for a given flow velocity, it is possible to minimize the difference between measured erosion rates and calculated ones with Equations (24)-(28) by varying parameter λ. The points in Figures 7 and 8 shows the optimal erosion rates for each flow velocity from the measurements at Ballantrae Hill Country for different combinations of k C C and λ. The intersections of erosion rate isolines narrow the range of such combinations for the set of experiments on the same soil with different flow velocities (green borders at Figures 7 and 8). The soil at plot 1 is characterized by a cohesion force k C C in the range of 12-14 m 2 /s 2 (or C = 60-70 kPa), a parameter λ in the range of 0.0025-0.004 and the percentage of the area, occupied by erodible soils after 3600 s of the armoring-loosening cycles, in the range of 0.55-0.7%.
The same types of calculations were performed for the flow velocities and erosion rates, measured at plot 2 ( Figure 8). They show that the soil is characterized (see green boarders) by an initial cohesion force k C C in the range of 27-31 m 2 /s 2 (or C = 180-200 kPa), a parameter λ in the range of 0.0015-0.002, and the percentage of the area, occupied by erodible soils after 3600 s of the armoring-loosening cycles, in the range of 0.35-0.55% per second.  a parameter λ in the range of 0.0015-0.002, and the percentage of the area, occupied by erodible soils after 3600 s of the armoring-loosening cycles, in the range of 0.35-0.55% per second. These results are in general accordance with the measured characteristics at these two plots. The cohesion, measured with a tore-vane, was about 50 kPa at both plots, which is rather close to the estimates provided by modelling for plot 1. The much higher estimate for the cohesion for plot 2 is explained by the anti-erosion effect of the organic matter in the soil (2.4%). The parameter λ in the failure function, with controls the rate of soil loosening during erosion and the rate of increase in the area of eroded soil, which controls the rate of erosion, are both two-fold larger at plot 2 than at plot 1. This is also due to the anti-erosion effect of organic matter. These observations show that erosion at plot 1 was due to the destruction of within-aggregate links, and so small soil particles were moved from the soil surface. The erosion at plot 2 was due to the destruction of between-aggregate links, and so rather large particles were collected at the end of the plot (Figure 9). This difference in the erodibility of the soil at plots 1 and 2 is also well- These results are in general accordance with the measured characteristics at these two plots. The cohesion, measured with a tore-vane, was about 50 kPa at both plots, which is rather close to the estimates provided by modelling for plot 1. The much higher estimate for the cohesion for plot 2 is explained by the anti-erosion effect of the organic matter in the soil (2.4%). The parameter λ in the failure function, with controls the rate of soil loosening during erosion and the rate of increase in the area of eroded soil, which controls the rate of erosion, are both two-fold larger at plot 2 than at plot 1. This is also due to the anti-erosion effect of organic matter. These observations show that erosion at plot 1 was due to the destruction of within-aggregate links, and so small soil particles were moved from the soil surface. The erosion at plot 2 was due to the destruction of between-aggregate links, and so rather large particles were collected at the end of the plot (Figure 9). This difference in the erodibility of the soil at plots 1 and 2 is also well-illustrated by the aggregate stability tests: wet sieving aggregates for 30 min showed that 88% was retained for soil from plot 2 and only 7.6% was retained for soil from plot 1. illustrated by the aggregate stability tests: wet sieving aggregates for 30 min showed that 88% was retained for soil from plot 2 and only 7.6% was retained for soil from plot 1. The numerical experiments (Figures 7 and 8) show the possibility to estimate the main soil characteristics used in the model through gully erosion simulation in flume or direct measurements of the gully erosion rate in the field with a variety of flow velocities on the same soil. The number of such experiments is limited now, and there is a large field for future investigations.

Discussion
The combination of classical deterministic mechanics-Newton's second law-and the stochastic description of the erosion control factors (driving and resistance forces) has led to the development of the stochastic gully erosion model (STOGEM).
The mean erosion rate is calculated with the algorithm in Figure 1 using the equation of the force balance acting on a soil particle in a water flow with the PDFs of flow velocity, soil aggregate size and cohesion, including several model parameters, such as kU, kD and kC. The mathematical methods for the calculation of the PDF of the function of stochastic variables, knowing the PDFs of these stochastic variables [14], are used. The relationships between the main factors of erosion (for example, flow velocity) and erosion rate are theoretically determined within the model for a given combination of input data. This theory opens up a new way for better understanding the experimental results of soil erosion, and demonstrates the direction for future investigations. It is possible to estimate the main soil characteristics used in the model by gully erosion simulation in flume or direct measurements of the gully erosion rate in the field with a variety of flow velocities in the same soil. The results of numerical experiments are consistent with the measurements in the laboratory flumes [15] and show the same stages of erosion development in cohesive soils.
The relationships between the rate of erosion and control factors are not pre-installed in this model. This is the main advantage and its main difference from empirical USLE-type [28] and "process-based" shear stress-type models [27]. STOGEM also contains novel elements compared to previously published stochastic models of soil erosion [29][30][31][32].
The expression for cohesive soil erosion rate calculation (Equation (17)) can be naturally included into the equation of deformation in models such as GULTEM or DYNGUL [10], or any other process-based model of gully erosion instead of using empirical shear-stress type formulas. The numerical experiments (Figures 7 and 8) show the possibility to estimate the main soil characteristics used in the model through gully erosion simulation in flume or direct measurements of the gully erosion rate in the field with a variety of flow velocities on the same soil. The number of such experiments is limited now, and there is a large field for future investigations.

Discussion
The combination of classical deterministic mechanics-Newton's second law-and the stochastic description of the erosion control factors (driving and resistance forces) has led to the development of the stochastic gully erosion model (STOGEM).
The mean erosion rate is calculated with the algorithm in Figure 1 using the equation of the force balance acting on a soil particle in a water flow with the PDFs of flow velocity, soil aggregate size and cohesion, including several model parameters, such as k U , k D and k C . The mathematical methods for the calculation of the PDF of the function of stochastic variables, knowing the PDFs of these stochastic variables [14], are used. The relationships between the main factors of erosion (for example, flow velocity) and erosion rate are theoretically determined within the model for a given combination of input data. This theory opens up a new way for better understanding the experimental results of soil erosion, and demonstrates the direction for future investigations. It is possible to estimate the main soil characteristics used in the model by gully erosion simulation in flume or direct measurements of the gully erosion rate in the field with a variety of flow velocities in the same soil. The results of numerical experiments are consistent with the measurements in the laboratory flumes [15] and show the same stages of erosion development in cohesive soils.
The relationships between the rate of erosion and control factors are not pre-installed in this model. This is the main advantage and its main difference from empirical USLE-type [28] and "process-based" shear stress-type models [27]. STOGEM also contains novel elements compared to previously published stochastic models of soil erosion [29][30][31][32].
The expression for cohesive soil erosion rate calculation (Equation (17)) can be naturally included into the equation of deformation in models such as GULTEM or DYNGUL [10], or any other process-based model of gully erosion instead of using empirical shear-stress type formulas.
The theoretical approach to gully erosion calculations is far from achieving direct application in soil conservation practice. There are many unknown processes and parameters, mostly related to changeable and complicated matters, such as the soil. The list of possible improvements and clarifications of the model is long. It is obvious that most of the parameters taken as constants in the above-described model (soil particle shape, porosity, etc.) are stochastic variables. The exponential failure function is the simplest option to model the process of soil loosening, and the wet-sieving test of aggregate stability shows that the Weibull function better explains this process. The different resistance to erosion of the links within and between soil aggregates must be investigated. The effect of raindrops on erosion by water in the gullies must also be taken into account.
The general limitation of this model is the requirement of erosion of the surface layer particle by particle. This requirement may not be respected at high rates of soil loosening and high velocities. The processes of transport and the deposition of eroded soil particles are not taken into account, and the influence of such processes can be important [33]. The further comparison of the measured rates of gully erosion performed in the frame of this model with calculated ones will reveal other details necessary for stochastic modelling.

Conclusions
The stochastic gully erosion model (STOGEM) is based on a combination of deterministic mechanics-Newton's second law-and the stochastic description of the erosion control factors. The main proposition in this model is that the depth of the active surface layer of eroded cohesive soil is equal to one particle diameter, and erosion develops particle by particle. The processes of transport and deposition of eroded soil particles are not taken into account. The erosion rate at the gully bed is calculated directly from the equation of the balance between the driving and resistance forces acting on a soil particle in flowing water. All factors of erosion-flow velocity, soil aggregate size and cohesion-are regarded as stochastic variables, described by probability density functions. The probability density function of cohesion in the model varies through time and space during an erosion event due to the changes in soil composition, namely armoring and loosening. This theory is still far from achieving practical application, but opens up a new way for better understanding the experimental results of gully erosion and demonstrates the direction for future investigations.