The Potential of Gully Erosion on the Yamal Peninsula, West Siberia

: The Yamal Peninsula occupies the northern part of the West Siberian Plain in Russia. This territory has rapidly developed due to the exploitation of several gas ﬁelds. At the same time, the Yamal Peninsula is one of the most severely gullied landscapes in the Arctic. The potential risk of damage to the environment or structures and the cost of such damages are very high there. The erosion potential is the cumulative erosion by runo ﬀ above critical, calculated for each point at a catchment. Calculations take into account the geomorphic, lithological, and vegetation cover thresholds, realized in the form of critical runo ﬀ depth of erosion initiation. It also takes into account action of all ﬂows between the critical and maximum runo ﬀ . The calculations for several gullied catchments on the Yamal Peninsula show the uneven distribution of erosion potential level with the maximum of gully erosion on the steep banks of the river valleys and on gully heads with bare soil. The area with potential erosion in these catchments varies within the range of 17–33%. The erosion on the Yamal Peninsula is mainly of natural origin. It occurs on steep slopes and at the heads of gullies. These landforms are not used for exploitation camps and settlements. Nevertheless, the linear structures, such as railways, roads and pipelines, can cross these unstable landforms with the risk of damage. Erosion potential increases at the spots with bare soil, which appear due to both construction work and natural processes, such as slumping.


Introduction
Gully erosion has been and remains one of the most dangerous processes of land damage [1,2]. Although it mainly occurs on agricultural lands [3], in recent years more and more gullies are formed in urbanized and industrial areas [4,5]. A typical example is the territory of the Yamal Peninsula where lands previously used for pastures only by local residents [6] are actively used at present in the development of gas and oil fields. It is important to find methods for the development of these lands that would not lead to significant changes in the vegetation cover, and in unstable landforms.
The Yamal Peninsula occupies the northern part of the West Siberian Plain (Figure 1). The peninsula is surrounded by the waters of the Baydaratskaya Bay and Kara Sea from the west and north, and of the Ob Gulf in the east. The area of the peninsula is 122,000 km 2 , and it is 750 km long from north to south and 140-240 km wide. The territory of the peninsula is open to the air masses from Arctic Ocean and the continent. The summer is short and cold, the winter is windy and frosty. The average annual temperature in the south of Yamal is −6.6 • C, and in its northern part −10.2 • C [7]. The peninsula is located within the permafrost zone. The depth of the seasonally thawed layer varies from 0.2 m in the northern part of Yamal to 2.0 m in the south. The entire Yamal Peninsula is a lowland plain, composed of marine loam and alluvial sands. The predominant altitudes are 40-70 m above sea level, the highest point is 102 m. The Yamal peninsula is an area of rapid technogenic development due to gas field exploitation. There are several fields already prepared for exploitation, and more will appear in future. At the same time, the Yamal Peninsula is one of the most severely gullied landscape in the Russian Arctic ( Figure 2) with the density of gullies is up to 1-2 km/km 2 . Therefore, natural environment, buildings The Yamal peninsula is an area of rapid technogenic development due to gas field exploitation. There are several fields already prepared for exploitation, and more will appear in future. At the same time, the Yamal Peninsula is one of the most severely gullied landscape in the Russian Arctic ( Figure 2) with the density of gullies is up to 1-2 km/km 2 . Therefore, natural environment, buildings and structures potentially can be damaged by gully erosion and the cost of such damage can be very high.
Sustainability 2020, 12, x FOR PEER REVIEW 3 of 16 and structures potentially can be damaged by gully erosion and the cost of such damage can be very high. The methodology of evaluation of gully erosion potential is related to the ideas of geomorphic thresholds and magnitudes and frequencies of erosion events. In the classic paper "Geomorphic thresholds: the concept and its applications". Schumm [8] argued: "geomorphic thresholds can be of two types, and they can be redefined in the following way. A geomorphic threshold is a threshold of landform stability that is exceeded either by intrinsic change of the landform itself, or by a progressive change of an external variable". Schumm [8] presented a gully system development as an example of both intrinsic and extrinsic threshold conditions. Gully head erosion and gully length growth occur in the condition of exceeding of an intrinsic geomorphic threshold. After drainage area and/or slope steepness above a gully head decreases below a threshold values, a gully head erosion and gully length growth stops. The value of this threshold can increase or decrease due to external variable changes: surface runoff or vegetation cover transformation. In the case of threshold decrease gully can become unstable again.
In the other classic paper, Wolman and Miller [9] wrote: "Almost any specific mechanism requires that a certain threshold value of force be exceeded. However, above this threshold or critical limit, there occurs a wide range in magnitude of forces which results from variations in intensity of precipitation, wind speed, etc.". "Available evidence indicates that evaluation of the effectiveness of a specific mechanism and of the relative importance of different geomorphic processes in molding specific forms involves the frequency of occurrence as well as the magnitude of individual events".
According to this statement, a measure of the potential risk that the natural environment or structure (if it exists) will be damaged by gullies can be taken proportional to the combination of the values and frequencies of all erosion events that exceed a certain geomorphic threshold above which this form of relief becomes unstable [3]. This combination is the amount of cumulative erosion during various events, or the quantity of sediment transported.
Patton and Schumm [10] proposed the formula to calculate threshold conditions for gully initiation that includes "intrinsic change of the landform itself" and "a change of an external variable".  The methodology of evaluation of gully erosion potential is related to the ideas of geomorphic thresholds and magnitudes and frequencies of erosion events. In the classic paper "Geomorphic thresholds: the concept and its applications". Schumm [8] argued: "geomorphic thresholds can be of two types, and they can be redefined in the following way. A geomorphic threshold is a threshold of landform stability that is exceeded either by intrinsic change of the landform itself, or by a progressive change of an external variable". Schumm [8] presented a gully system development as an example of both intrinsic and extrinsic threshold conditions. Gully head erosion and gully length growth occur in the condition of exceeding of an intrinsic geomorphic threshold. After drainage area and/or slope steepness above a gully head decreases below a threshold values, a gully head erosion and gully length growth stops. The value of this threshold can increase or decrease due to external variable changes: surface runoff or vegetation cover transformation. In the case of threshold decrease gully can become unstable again.
In the other classic paper, Wolman and Miller [9] wrote: "Almost any specific mechanism requires that a certain threshold value of force be exceeded. However, above this threshold or critical limit, there occurs a wide range in magnitude of forces which results from variations in intensity of precipitation, wind speed, etc.". "Available evidence indicates that evaluation of the effectiveness of a specific mechanism and of the relative importance of different geomorphic processes in molding specific forms involves the frequency of occurrence as well as the magnitude of individual events".
According to this statement, a measure of the potential risk that the natural environment or structure (if it exists) will be damaged by gullies can be taken proportional to the combination of the values and frequencies of all erosion events that exceed a certain geomorphic threshold above which this form of relief becomes unstable [3]. This combination is the amount of cumulative erosion during various events, or the quantity of sediment transported. Patton and Schumm [10] proposed the formula to calculate threshold conditions for gully initiation that includes "intrinsic change of the landform itself" and "a change of an external variable".
where S is the gully catchment slope gradient, A is the contributing area, and a and b the coefficient and exponent. The regional implications of Equation (1) showed broad variability of the coefficient a, and more stable value of the exponent b. Vandaele et al. [11] compiled the parameters of Equation (1) for 25 gullied basins and obtained the coefficient a in the range 0.02-0.5 and exponent b in the range 0.15-0.47. With the increasing of the number of investigated gullied basins, these ranges increase also ( [12][13][14] and many others). It is obvious that the applicability of empirical Equation (1) for gully erosion potential prediction is limited due to uncertainty of the parameters. This paper describes a novel methodology of gully erosion potential evaluation in the conditions of high Arctic and presents several examples of calculation for the gullied catchments on the Yamal Peninsula.

Theoretical Meanings of Empirical Equations
There is a large number of theoretical justifications of Equation (1). There are several ways to derive Equation (1) from the expression for the critical tangential stress in the flow τ, at which erosion begins [15][16][17][18], or from the expression for the critical flow velocity U cr [19,20]. The way used by Begin and Schumm [15,16] is following. If shear stress in flow is equal or more than its threshold value then the basin is gullied or potentially will be gullied in future. Tangential stress in a flow is the product of water density ρ, acceleration due to gravity g, flow depth d and flow slope S: and discharge Q is where M is runoff depth, uniform for a gully catchment, and k is the coefficient of transformation of units. Note that Begin and Schumm [15] used non-uniform M, which depended on contributing area A. These formulas transform coefficients in Equation (1) into Inequality (7) shows the critical spatial-temporal conditions in the basin to be potentially eroded at a given point with inclination S cr and contributing area A cr at a given time when surface runoff depth of some probability is equal to M cr . It is also possible to use the critical velocity of erosion initiation as a threshold condition.
Quantitatively, expression (8) is easy to derive from Equation (2). Nevertheless, critical velocity better describes erosion initiation, as the driving forces (lift and drag) depend on flow velocity, not on shear stress [19]. Flow velocity is represented by Chezy-Manning formula Here, n is roughness coefficient. Then coefficients in Equation (1) can be calculated as and inequality (8) takes the form There are several other ways to derive Equation (1) from the formulas of hydraulics. For example, Zorina [20] transformed Equation (9) into Width/depth ratio is the function of discharge Therefore, The changes in exponent b is controlled by the method used to derive Equation (1), as well as by regional variability of the exponent m or f [21]. The factors of the changes in coefficient a are more complicated. This coefficient depends on the critical shear stress or velocity and on the surface runoff depth in addition to the factors of variability of the exponent b.

The Novel Approach to Erosion Potential Estimation
One of the main problems of calculations with Equation (7) or Equation (11) is an indefinite value of surface runoff depth for estimation of critical velocity or shear stress. The proposed methodology aims to avoid this uncertainty. According to the proposed novel approach, the inequality (11) (or any other formulation of Equation (2) or Equation (8)) can be transformed in such way that the indefinite critical surface runoff depth becomes the function of known arguments: Inequality (15) shows that a gully catchment will be eroded at a given point with inclination S cr and contributing area A cr when the surface runoff depth M is equal or more than the critical value. The ground and vegetation cover conditions on this catchment correspond to critical velocity of erosion initiation U cr .
The potential erosion R (depth of erosion at a given point) is the cumulative erosion during subcritical runoff.
Here E M is erosion rate under all given runoff depths M during time interval dt, T-T 0 is duration of surface runoff depth on the catchment above critical M cr for a given point. According to [22], Here, k e is erodibility coefficient, q-specific discharge Q/W. The channel width can be represented as W = p w Q mw , then Gully erosion potential is calculated with the Equation (16) for those lines on the catchment, where flow velocity is more than critical, and erosion by water can be initiated. Outside these lines, erosion by water cannot be initiated at a given morphological, hydrological, lithological and vegetation cover conditions. All the processes of gully formation related to water flow can occur only within the area with initially critical trigger factors, although these factors (flow velocity and slope) can change through time. It is in this area that there is a risk of linear erosion by a concentrated flow.
The erosion by water in the Arctic conditions develops in two ways: mechanical erosion described by Equation (17) and thermal erosion, when frozen soil is thawed by more warm water, and the thawed layer is removed by flow. Investigation showed [23] that thermal erosion rate by a given flow is always less than values calculated with Equation (17).
The second important process, which is not described with Equation (17), is slumping of thawed soil from the gully walls to the bottom. This process in general decreases the rate of gully incision as the gully bottom occasionally rise due to accumulation of slumped sediments. Therefore, Equation (17) describes the upper limit of the combination of erosion rate and erosion factors.
Slumping from gully walls increases the gully top width and, to some extent, the gully length. This process is slower than erosion by water. It becomes significant at the late stages of gully development. With the general stabilization of a gully, its slopes reach the angle of repose φ. Then the gravitational processes stop. Nevertheless, if we take into account gravitational processes, this increases the estimate of the area with gully erosion potential, calculated with Equation (16). This increase should be taken into account at the late stages of the development of the gully. For example, a gully can become longer by ∆L = d H cot(φ), where d H is gully depth at the head. For the stable gullies on the Yamal Peninsula the typical angle of repose φ is 0.6 radians, gully depth at the head is in the range 2-4 m [23]. Therefore, an additional length of stable gullies compares to the calculated with Equation (16) is around 3-6 m.
Potential erosion in Equation (16) for a given catchment is controlled by seven local parameters (M max , U cr , n, m, p, m w , p w , k e ), three variables (S cr , kA cr , M cr ) distributed in space, and one variable (M) changing through time t between T 0 and T. There are different ways to estimate these variables and parameters for the Yamal Peninsula. Two of these variables (slope S and contributing area A) were taken from a digital elevation model, six parameters (U cr , n, m, p, m w , p w , k e ) were derived from field experiments, and other parameters and variables (M max , M cr and M) were calculated with appropriate models.
The digital elevation models of the catchments were extracted from ArcticDEM [24] with the horizontal resolution 2 m. Modules of geographic information system QGIS with GIS SAGA [25] were used sequentially to fill sinks (SAGA Terrain Analysis-Hydrology-Fill sinks), to calculate flow accumulation (SAGA Terrain Analysis-Hydrology-Catchment area) and to calculate slope (SAGA Terrain Analysis-Morphology-Slope, aspect, curvature).
Field measurements in 1990-97 and in 2007 showed [23] that for the gullies on the Yamal Peninsula the exponent m is 0.3 and coefficient p is 0.21, the exponent m w is 0.45 and coefficient p w is 3.4 ( Figure 3). A typical roughness coefficient n is 0.06 (see Table 4 in ref. [23]). Erodibility coefficient vary from 5.3 to 18.6, depending on soil type (see Figures 2 and 3 in ref. [23]).  The conditions of the erosion initiation for cohesive soils, which are typical for the gully erosion, are well investigated for a broad variety of soil types, both in laboratory and in the field [19]. It is possible to select an appropriate value for a given environment. According to field investigations on the Yamal Peninsula, the erosion under natural vegetation cover is negligible. Only bare ground is eroded by mechanical and thermal action of water. Measurements in [26] showed that erosion initiation in the loam deposits typical of the most part of Yamal Peninsula begins at the velocity 0.2 m/s. Therefore, critical daily surface runoff in mm is: and the erosion rate in kg/(m 2 s) is

Hydrological Modelling
The systematic hydrological measurements are absent for the Yamal Peninsula, only some information from a few short expeditions is available [27,28]. Therefore, hydrological modelling is the main source of information on the changes in surface runoff depth through time and its frequency distribution in this region. Meteorological data for hydrological modelling are available from a few stations on the peninsula or from the reanalysis grid. In this paper, Era Interim reanalysis [29] with the grid step of 0.75° was used for the hydrological calculations.
In the conditions of Arctic climate in the region with deep permafrost, hydrological modelling bears specific characteristics. The main simplification of such modelling is the possibility to neglect water infiltration into frozen or highly saturated soil. Two main sources of water flow are typical for this environment: snow thaw during the spring and rainfall, mostly during the summer. The surface runoff depth Ms for the period of snow thaw was calculated with the models from [30][31][32][33] using the depth of snow pack (water equivalent), precipitation, air temperature and evaporation. The The conditions of the erosion initiation for cohesive soils, which are typical for the gully erosion, are well investigated for a broad variety of soil types, both in laboratory and in the field [19]. It is possible to select an appropriate value for a given environment. According to field investigations on the Yamal Peninsula, the erosion under natural vegetation cover is negligible. Only bare ground is eroded by mechanical and thermal action of water. Measurements in [26] showed that erosion initiation in the loam deposits typical of the most part of Yamal Peninsula begins at the velocity 0.2 m/s. Therefore, critical daily surface runoff in mm is: and the erosion rate in kg/(m 2 s) is

Hydrological Modelling
The systematic hydrological measurements are absent for the Yamal Peninsula, only some information from a few short expeditions is available [27,28]. Therefore, hydrological modelling is the main source of information on the changes in surface runoff depth through time and its frequency distribution in this region. Meteorological data for hydrological modelling are available from a few stations on the peninsula or from the reanalysis grid. In this paper, Era Interim reanalysis [29] with the grid step of 0.75 • was used for the hydrological calculations.
In the conditions of Arctic climate in the region with deep permafrost, hydrological modelling bears specific characteristics. The main simplification of such modelling is the possibility to neglect water infiltration into frozen or highly saturated soil. Two main sources of water flow are typical for this environment: snow thaw during the spring and rainfall, mostly during the summer. The surface runoff depth Ms for the period of snow thaw was calculated with the models from [30][31][32][33] using the depth of snow pack (water equivalent), precipitation, air temperature and evaporation. The algorithm of calculating the runoff during the summer rains is based on precipitation depth and evaporation from Era Interim reanalysis with additional losses due to filling the depressions on the catchment surface calculated after Popov [34].
Hydrological modelling shows that the duration T D (days) of daily surface runoff depth ≥ M more than 1 mm is well approximated by an exponent: Therefore, the module of time interval dt of daily surface runoff depth M is Using this approximation for surface runoff frequency, the formula for calculating the gully erosion potential level (GEPL) takes the following form: where γ is the lower incomplete gamma function. The logarithms in Equation (23) Sustainability 2020, 12, x FOR PEER REVIEW 8 of 16 evaporation from Era Interim reanalysis with additional losses due to filling the depressions on the catchment surface calculated after Popov [34]. Hydrological modelling shows that the duration TD (days) of daily surface runoff depth ≥ M more than 1 mm is well approximated by an exponent: Therefore, the module of time interval dt of daily surface runoff depth M is Using this approximation for surface runoff frequency, the formula for calculating the gully erosion potential level (GEPL) takes the following form: where γ is the lower incomplete gamma function. The logarithms in Equation (23) Value -R0min is added to Equation (23) to avoid negative erosion potential levels. The minimum value of R 0 (when M cr = M max ) is minus infinity. It is necessary to assign some maximum M cr slightly lower than M max . If M cr is 99% of M max then R 0min = ln[γ(1.55; B N M max ) − γ(1.55; 0.99B N M max )]. (25) Value -R 0min is added to Equation (23) to avoid negative erosion potential levels. It is obvious that at the points with critical runoff depth M cr exceeding the annual maximum of the surface runoff M max the erosion rate is zero. Therefore the mean gully erosion potential (MGEP) R for a given territory or a catchment can be calculated as the ratio of number of points (pixels) with M cr < M max to the total number of points (pixels).

Model Validation
It is not easy to validate the results of erosion potential calculations, as this potential almost never occurs in the natural conditions, except in the badlands. Usually erosion potential is calculated to prevent dangerous results of erosion by soil conservation measures and in practice gully erosion stop at some stage of its development. The second difficulty in validation procedures is the inaccuracy of existing DEMs used for erosion potential calculations. Linear erosion can follow along such lines on the catchment as vehicle tracks, locally damaged vegetation, the local sources of water, etc., which cannot be implemented into DEM of any resolution.
The common definition of validation offered by Society for Computer Simulation Technical Committee on Model Credibility [35] 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 difficulties mentioned above. The territory of Yamal Peninsula is not badlands, although in several parts with erosion net density about 2 km/km 2 it is close to such condition. ArcticDEM, the base for our calculations, depends on the resolution of the images used for DEM construction. Therefore, the model validation can be only approximate. The catchments selected for the validation procedure (point 2 in Figures 1 and  5A) are situated close to the severely gullied part of the peninsula. The human activity at the KEKH exploitation camp (Bovanenkovskoye gas field) in 1986-1997 led to a nearly complete destruction of the vegetation cover there. Images from 1986-1995 and up to 2016 show that several bank gullies [23], which already became stable in virgin conditions, resumed their growth in length, mostly in the form of deep rills and ephemeral gullies ( Figure 5A). At these catchments, the distribution of erosion potential level ( Figure 5B) calculated with Equation (23) and U cr = 0.2 m/s is quite similar to the erosion net, which already existed in 1995 ( Figure 5C). Although there is no complete spatial overlapping of these two nets, the qualitative picture shows a consistence between the model results and measurements and "a satisfactory range of accuracy". The mean erosion potential for this territory is 70%, which corresponds to the degree of the damage to the vegetation cover.

Results
The distribution of erosion potential level was calculated for seven small gullied catchments on the Yamal Peninsula (see Figure 1, Table 1). Hydrological characteristics were calculated for the nearest nodes of Era Interim net for the period 1979-2018. The maximum of the daily surface runoff Mmax vary at those points in the range of 74-118 mm. For investigated catchments the value of -R0min = 4 was used in calculations.
The distribution of the GEPL calculated with Equation (23) for the gully basin at the headwaters of a small river Nekhedeyvar-Yakha in the north of the Yamal Peninsula (72.85° N, 70.87° E) is highly uneven (Figures 6,7).
A combination of critical surface runoff, local slope and contribution area for a given critical velocity of erosion initiation develops a complex pattern of erosion potential level. The maximum GEPL is associated with steep banks of river valleys and with gully heads. These landforms are vulnerable to shallow landslides which damage the vegetation cover [36]. In this case, we used in the

Results
The distribution of erosion potential level was calculated for seven small gullied catchments on the Yamal Peninsula (see Figure 1, Table 1). Hydrological characteristics were calculated for the nearest nodes of Era Interim net for the period 1979-2018. The maximum of the daily surface runoff M max vary at those points in the range of 74-118 mm. For investigated catchments the value of -R 0min = 4 was used in calculations.
The distribution of the GEPL calculated with Equation (23) for the gully basin at the headwaters of a small river Nekhedeyvar-Yakha in the north of the Yamal Peninsula (72.85 • N, 70.87 • E) is highly uneven (Figures 6 and 7).
A combination of critical surface runoff, local slope and contribution area for a given critical velocity of erosion initiation develops a complex pattern of erosion potential level. The maximum GEPL is associated with steep banks of river valleys and with gully heads. These landforms are vulnerable to shallow landslides which damage the vegetation cover [36]. In this case, we used in the calculations critical velocity for erosion initiation 0.2 m/s typical for bare ground conditions ( Figure 6).
The risk of gully erosion on the same catchment in the conditions of well-preserved tundra vegetation (U cr = 0.5 m/s) is significantly lower (Figure 7). Still, the heads of gullies bear a high potential of erosion.
Sustainability 2020, 12, x FOR PEER REVIEW 11 of 16 calculations critical velocity for erosion initiation 0.2 m/s typical for bare ground conditions ( Figure  6). The risk of gully erosion on the same catchment in the conditions of well-preserved tundra vegetation (Ucr = 0.5 m/s) is significantly lower (Figure 7). Still, the heads of gullies bear a high potential of erosion.  The calculations for seven small catchments situated in different parts of the Yamal Peninsula (Table 1) show similar, highly uneven distributions of GEPL with the same maximums on the steep banks of the river valleys and in gully heads. Table 1. The main characteristics of gully erosion potential level calculated for the gullied catchments on the Yamal Peninsula (for positions see Figure 1).  The calculations for seven small catchments situated in different parts of the Yamal Peninsula (Table 1) show similar, highly uneven distributions of GEPL with the same maximums on the steep banks of the river valleys and in gully heads. Table 1. The main characteristics of gully erosion potential level calculated for the gullied catchments on the Yamal Peninsula (for positions see Figure 1).

Discussion
The gully erosion potential level R L for a given catchment is the sum of two logarithms: logarithm of product of slope S and contributing area A powered on 0.55 and logarithm of weighted duration of surface runoff depth above M cr . The first component depends on paired distribution of slope and contributing area and is invariant for a given catchment. The second component depends on M cr , which is also the function of slope and contributing area, as well as on parameters M max , U cr , n, m, p, m w , p w , k e . There are several options of estimation M cr . We shall discuss only three of them that use Chezy-Manning formula (9) and critical flow velocity U cr .
In option 1, the regime equation used in Equation (9) If in Equation (9) Using m w = 0.45 and p w = 3.4, we obtain Calculations for the gully basin at the headwaters of a small river Nekhedeyvar-Yakha in the north of the Yamal Peninsula (72.85 • N, 70.87 • E) with these three options show that the maximum differentiation of erosion potential is achieved if using option 1 (Equation (27) for M cr ) ( Figure 8). Therefore, this option was used in all calculations described in the paper. All three options are very sensitive to variations in critical velocity, so that possible error of Mcr calculation is 4.55-5-times larger than the error of Ucr estimation. The error affects both the mean gully erosion potential and distribution of its level. Sensitivity of calculations to the value of critical velocity must be taken into account in the situations when the quality of vegetation cover becomes one of the main factors of gully erosion. This is important when the restoration of territories with industrial damage is achieved by the improvement of vegetation cover. In such cases, the calculations of gully erosion potential should be accompanied with careful estimations of critical velocity of erosion initiation for a given quality of vegetation cover [37,38].

Conclusions
The novel measure of gully erosion potential can be taken proportional to the combination of magnitudes and frequencies of all erosion events exceeding some threshold, above which a given landform becomes unstable [9]. This combination was implemented into Equation (16) and represents the erosion depth. This equation takes into account the geomorphic, lithological and vegetation cover thresholds, realized in the form of critical runoff depth Mcr of erosion initiation (Equation (15) or (27)). It also takes into account action of all flows between the values of this threshold and maximum runoff depth Mmax. The logarithm of potential erosion depth was assigned as gully erosion potential level (GEPL) at a given point (Equation (23)), and the ratio of number of points (pixels) with Mcr < Mmax to the total number of points (pixels) was assigned as mean potential of gully erosion (MPGE) for entire territory or catchment. Model validation shows the consistence of calculations with the All three options are very sensitive to variations in critical velocity, so that possible error of M cr calculation is 4.55-5-times larger than the error of U cr estimation. The error affects both the mean gully erosion potential and distribution of its level. Sensitivity of calculations to the value of critical velocity must be taken into account in the situations when the quality of vegetation cover becomes one of the main factors of gully erosion. This is important when the restoration of territories with industrial damage is achieved by the improvement of vegetation cover. In such cases, the calculations of gully erosion potential should be accompanied with careful estimations of critical velocity of erosion initiation for a given quality of vegetation cover [37,38].

Conclusions
The novel measure of gully erosion potential can be taken proportional to the combination of magnitudes and frequencies of all erosion events exceeding some threshold, above which a given landform becomes unstable [9]. This combination was implemented into Equation (16) and represents the erosion depth. This equation takes into account the geomorphic, lithological and vegetation cover thresholds, realized in the form of critical runoff depth M cr of erosion initiation (Equation (15) or (27)). It also takes into account action of all flows between the values of this threshold and maximum runoff depth M max . The logarithm of potential erosion depth was assigned as gully erosion potential level (GEPL) at a given point (Equation (23)), and the ratio of number of points (pixels) with M cr < M max to the total number of points (pixels) was assigned as mean potential of gully erosion (MPGE) for entire territory or catchment. Model validation shows the consistence of calculations with the measurements and "a satisfactory range of accuracy".
The calculations for seven gullied catchments on the Yamal Peninsula (Table 1) shows the highly uneven distribution of erosion potential level with the maximums on steep banks of the river valleys and in gully heads. The same landforms are vulnerable to shallow landslides, which damage the vegetation cover. Therefore, critical velocity for erosion initiation 0.2 m/s, which is suitable for bare ground, was used in calculations. The histograms of erosion potential level above zero are unimodal with probability maximums at 9.7-10.4. The mean erosion potential for these territories vary in a rather narrow range 0.17-0.33. This erosion potential is mostly natural, linked to the landforms with steep slopes and unstable vegetation cover. Such landforms are not used for exploitation camps and settlements. Nevertheless, the linear structures, such as railways, roads and pipelines, can cross these unstable landforms with high level of risk to be damaged. This erosion potential increases due to appearance of the spots of bare soil during construction works. This was shown for the territory of KEKH exploitation camp, where MPGE is 0.7. Such spots of bare soil can appear even on flat territories with an initially low level of erosion potential. Further work should be focused on calculations of erosion potential level for the territories of future development of gas exploitation to avoid possible damage to the environment and infrastructure.
Funding: This research was funded by RFBR grant 18-05-60147 "Extreme hydrometeorological phenomena in the Kara Sea and the Arctic coast".