PERIGLACIAL GULLY EROSION ON THE EAST EUROPEAN PLAIN AND ITS RECENT ANALOG AT THE YAMAL PENINSULA

. The net of dry valleys, gullies and shallow hollows is typical for the East European Plain. Dense vegetation usually covers their bottoms and slopes, so the modern erosion there is negligible in the pristine conditions. This erosion landscape formed in periglacial conditions during the terminations of the last two glaciations. The same kind of the erosion landscape is typical for the Arctic regions, especially for the Yamal, Gydan, and Tazovsky peninsulas. The size and the density of such valleys and gullies are quite similar to those existing on the East European Plain, but these erosion features are active there, especially in the conditions of natural or anthropogenic deterioration of the vegetation cover. As the density of dry valley network is an indicator of hydrological conditions in the river basin, the landscapes of the Arctic regions can be used as the modern analogs of the territories with the past periglacial erosion. The recent hydrological characteristics of the west-central Yamal Peninsula were used to estimate the parameters of erosion network at the Khoper River basin, formed in periglacial conditions. For these purposes gully erosion and thermoerosion model GULTEM was verified and calibrated based on the observation of the modern processes on the Yamal Peninsula. The meteorological characteristics were taken from ERA-Interim Reanalysis grid. To calculate the flow characteristics a synthetic hydrological model was used. These verified and calibrated models were used to find the most suitable characteristics of climate and vegetation cover, which can explain the structure and density of the Perepolye dry valley in the Khoper River basin. This dry valley with the main trunk length of 6400 m was formed at the end of the Late Valdai Glaciation (MIS 2). The conditions required for the formation of a periglacial gully of such length were estimated with the GULTEM model. The critical velocity of erosion initiation was within the range 0.8-0.9 m/s, and the surface runoff depth was close to the recent one on the Yamal Peninsula (330 mm). The system of shallow hollows in the Perepolye catchment (the gullies formed at the end of the Moscow Glaciation, MIS 6) is denser and longer than the dry valley system, and the modelling estimates showed that the surface runoff during that period was almost 3.3 times more than the recent one on the Yamal Peninsula.


INTRODUCTION
One of the most interesting relief features of the East European Plain is the net of dry valleys (Fig. 1A). These dry valleys typically have lengths in the range of 1-10 km. At present, these systems drain water mostly during snow thaw period or after heavy rains. Dense vegetation usually covers their bottoms and slopes, so that in the pristine conditions the modern erosion there is negligible.
The main part of these dry valleys (the Russian term is «balka») was formed during erosional events of the Late Quaternary (Eremenko et al. 2011;Matoshko 2012) in the periglacial conditions with the permafrost. This type of ancient erosion features, some buried, some renewed by the modern erosion, is widespread in Europe (Larsen et al. 2016;Hošek et al. 2019), Northern America (Bettis et al. 1985), and Australia (Prosser et al. 1994). The reconstruction of the environmental conditions of Quaternary gullies initiation and evolution are critical for understanding their morphology and history. The erosion landscapes in the Arctic regions are quite similar to those on the East European Plain in their morphology, but these erosion features are currently active in the Arctic (Fortier et al. 2007), especially in the conditions of natural or man-induced deterioration of vegetation cover (Sidorchuk 2015). The size and the density of such valleys and gullies, especially at Yamal, Gydan and Tazovsky peninsulas in Russia (Fig. 1B) are nearly the same as those of the net of «balkas» in the lowland regions of the East European Plain, such as the upper Dnepr and Don River basins. As the erosion network density depends on the hydrological conditions on the watershed (Carlston 1963;Gregory 1966;Gregory et al. 1968), the erosion processes at the gully catchments in the Arctic regions can be used as the modern analogs of the former periglacial erosion.
The main goal of this paper is to determine the most suitable characteristics of surface runoff and vegetation, which can explain the density of the ancient dry valleys systems in the Don River basin on the East European Plain. For this purpose, the modern erosion features on the Yamal Peninsula were used as the closest morphological analog. Calculations of erosion net density were performed with the gully erosion and thermoerosion model (GULTEM). The model was verified and calibrated using observations on the modern hydrological and erosion processes in the Arctic. The climatic and meteorological characteristics for such verification are available through reanalysis.

The dry valley systems in the Don River basin
The Don River basin on the East European Plain was not glaciated at least for the last 250 ka. Therefore, mostly erosional processes determined the landscape there, and dry valley systems are of the main importance for the long-term palaeohydrological reconstructions. There were at least two major hydrological and erosional events, which governed the morphology of erosional landscape: the catastrophic surface runoff at the end of the Moscow Glaciation (MIS 6), about 130 ka ago, and the one at the end of the Late Valdai (MIS 2), 16-19 ka ago. The dry valleys formed during these events were described by Eremenko and Panin (2011). We used the dry valley Perepolye as a typical example for palaeohydrological reconstructions (Sidorchuk et al. 2018).
The Perepolye dry valley with a catchment area of 41.7 km 2 is situated on the left bank of the Khoper River, the left tributary of the Don River, near Povorino settlement (Fig. 2). The modern dry valley with the main trunk length of 6400 m is characterized by steep slopes (10-15°) and a concave longitudinal profile; the valley depth is 6-8 m in the middle part and 8-10 m in the lower part. Together with the tributaries, the dry valley forms a tree-like network with a total length of 11.3 km, and the density of the modern network is 0.27 km / km 2 .
In the upper part of the catchment, the modern dry valley and its tributaries are continued by shallow hollows, clearly visible on the plowed fields by their brighter tone on the space images. The heads of such hollows almost reach the watershed (see Fig. 2). The length of the entire dry valley-hollow network is 40.5 km and its density of 0.97 km / km 2 ; its longest element (from the mouth to the top of the main dry valley) reaches 8400 m.
To study the sediments within this dry valley, the coring along nine transverse profiles was performed. The sedimentary structures uncovered by drilling indicate the presence of two major buried erosion features (Fig. 3). The first (most ancient) erosion feature is embedded into the sands of the third terrace of the Khoper River. Its bed is incised 5-6 m be-low the bottom of the modern dry valley in its lower and middle parts and 2-3 m below it in the upper part. This buried gully has an average longitudinal slope of about 2.5 m/ km. On its bed buried soil layers with well-defined humus and illuvial horizons, with numerous filled animal burrows («krotovinas») exist. This buried soil (Fig. 4) probably formed during the Mikulino Interglacial (=Eemian, MIS 5e), since its structure corresponds to the Mezin soil complex (Velichko et al. 2007). Therefore, the incision event most probably took place at the end the Moscow Glaciation (MIS 6).
This first erosion feature was filled later almost to the top by the loess-like yellow-brown silt deposits with the maximum thickness over 13 m. These deposits are dated by OSL method to about 65 ka ago . The hollows in the upper part of the catchment were also formed during the late Moscow (MIS 6) erosion event and subsequently filled with similar loess-like yellow-brown sediments.
The second, younger erosion feature has a steeper longitudinal profile and is narrower than the first one. The average slope of its bottom is 3.7 m/km; the width is 20-30 m at the head of the dry valley and up to 100 m near its mouth. In the lower part of the dry valley, the loess-like deposits filling the older incision were completely eroded along the thalweg; they preserved only at the sides of the younger incision. Here the depth of this incision is about 20 m.
The second incision is filled with gray and black loams and clays with sandy interlayers. The top of these deposits forms the main part of the bottom of the modern dry valley. The AMS radiocarbon dating showed that the age of black heavy loams and clays which overlie the sands of the third terrace is 11441 +/-60 14C yr BP (AA104015). This agrees with the date 11900 +/-120 (Ki-5305), obtained for the deposits filling a similar incision in the nearby dry valley Khoprets (Sidorchuk and Borisova 2000).
The second erosional feature can be attributed to the powerful surface runoff event at the beginning of the degradation of the last glaciation, 16-19 ky ago, which left numerous traces in the relief of the periglacial zone of the Northern Hemisphere (Panin et al. 2006). The modern dry valley network was formed during this erosion event. The network of hollows, which was formed during the previous erosion event, was not affected by erosion at this stage. The modern dry valleys and hollows functioned as active deep gullies during the periods of catastrophic surface runoff. They were mainly in general inactive in the climatic and landscape conditions of the Holocene. The morphology of these systems including the maximum lengths of their main trunks can be regarded as the 'memory' of the fluvial system (Knighton 1984). This memory can be used to find the most probable characteristics of the hydrological regimen for the period of the dry valley systems formation. To resolve such inverse problem of geomorphology, an optimisation of the input characteristics of gully erosion model, which leads to the given output morphology of erosion landscape, is required.
The vegetation on the East European Plain during the periods of catastrophic surface runoff represented a so-called 'tundra-steppe' or 'periglacial forest-steppe' . The closest modern analogs of such past vegetation (and climate) on the East European Plain are found in the Altai and Sayan mountains (Borisova et al. 2006;Sidorchuk et al. 2011). As the erosion processes in the mountain regions cannot be applied as the analogs for erosion on lowlands, different indicators have to be used to find such analogs. At present, the lowland territories with the periglacial climate are situated in the Arctic and Subarctic regions of European Russia and Siberia. Therefore, the first approximation of input hydrological regimen and vegetation cover was taken from the west-central Yamal Peninsula. The main reason to use this analog is the similarity of erosion relief of the watersheds in the Don River basin and on the Yamal Peninsula, especially the similarity in drainage density, which is the main indicator of hydrological conditions (Gregory 1966;Sidorchuk et al. 2018).

The gully erosion in the west-central Yamal
Natural processes of the gully erosion are typical for the arctic tundra of the Yamal Peninsula (Sidorchuk 1996(Sidorchuk , 2015. The main natural factors of erosion in this region are: 1) sufficient annual precipitation (up to 350-400 mm); 2) low soil permeability due to the permafrost and therefore high runoff coefficients (up to 0.9-1.0); 3) high erodibility of frozen and thawing soils; 4) natural triggers of gully erosion, such as thermokarst processes and shallow landslides, which often cause vegetation cover deterioration.
This natural high erosion risk has been greatly increased recently by human impact. In the areas of gas production and

Fig. 4. Mezin-type buried soil complex on the erosional bed of the Perepolye dry valley
transportation facilities the erosion potential increases due to: 1) deterioration of the vegetation cover due to industrial development; 2) increased snow depth on the slopes due to excessive snow accumulation near buildings and roads; 3) an increase of the runoff coefficient on impermeable surfaces of the urbanized territories and roads; 4) local industrial and urban sources of warm water; 5) exploitation of sand-pits, gas and oil fields, and construction of pipe lines and ditches. The combination of high natural erosion potential and human interference causes extremely intensive gully, rill and sheet erosion (Fig. 5). Erosion follows the old drainage lines, mostly formed by the former shallow slumps, and new gullies cut into previously gently sloping elongate depressions. On the marine terrace in the west-central Yamal Peninsula, composed from frozen loams and sands, it takes anthropogenic gullies 4-10 years to reach their potential (possible) length.
The investigations of 1990-97 within the Scientific Program «Yamal» of RSC «GAZPROM» gave the opportunity to collect the information about hydraulic characteristics of the flows and morphological changes in the gullies, both natural and anthropogenic. The main feature is the alternating processes of quick intensive incision and rapid sidewall slumping in the gullies. Repeated instrumental measurements of the gully longitudinal profiles showed relatively low resulting mean rates of deepening. For the period of 1990-1995, the increase of the mean depth of gully K_1 was 0.9 m, and that of gully K_2 was 0.6 m. Note that during one rainfall on 8-9 Aug. 1990 the bed of the gully K_2 was incised by 0.58 m (Sidorchuk 2015). The discharge from the catchment with area 64300 m2 reached 21.3 l/s; the maximum sediment concentration was 74 kg/m3. This incision was already partly filled by sidewall slumping only 2-3 days after the flood. The mean rate of gully deepening is about 10 times less than the peak rate of incision during snowmelt and rainfall runoff, due to backfilling by material from slumping walls.
The gully increase in length is also a complicated process. The main head of gully K_1 remains stable since 1970, but in the 1990s a long and deep rill was formed in its middle part, and the second active head developed in the convex upper part of the slope. The length of gully K_2 was 165 m in 1988, 190 m in 1989, 210 m in 1990, 230 m in 1991, and 280 m in 1995. The growth rate of the main gully decreased in time, but in 1990-95, a new additional gully head was formed 120 m upstream from the old head, so that the gully became discontinuous. Nevertheless, both the photos taken in 2007 and space images of 2016 showed a general stabilisation of these gullies.

The modelling of gully erosion
The most typical for the gully erosion is to develop in two stages. At the first stage, gully channel formation is very intensive, and morphological characteristics of the gully (its length, depth, width, area, and volume) are changing rapidly. This stage takes about 5% of the whole gully lifetime, but about 90% of its length is formed during this stage (Kosov et al. 1978). To describe this first stage of gully evolution, we developed the gully thermoerosion-erosion model GULTEM. Water flow is calculated with the appropriate hydrological model. Meteorological data are taken from the reanalysis, as the station net is sparse in the Arctic.

Meteorological data from the reanalysis
Reanalysis is a grid data of surface and upper-air meteorological parameters. It spans a period over 30 years and can be used for climate analysis. The reanalysis combines observational data assimilation system and a forecast model. The atmospheric forecast model uses various hydrodynamic equations solved using finite-difference methods of forecast. Therefore, grid data of reanalysis are formed step-by-step using the assimilation data system and forecast model.
ERA-Interim Reanalysis is one of the best global atmospheric reanalyses (Lindsay et al. 2014). It was produced by the European Centre for Medium-Range Weather Forecasts (ECMWF). The atmospheric model in ERA-Interim has a horizontal resolution approximately 79 km (a grid of 0.75°×0.75°). The vertical resolution is inconstant and includes 60 layers with the highest level at 0.1 hPa. ERA-Interim covers the period from 1979 to the present. Time resolution is 3 hours for surface parameters and 6 hours for upper-air parameters. Atmospheric model of ERA-Interim assimilates many observation data: various types of satellite data (Dee et al. 2009), surface observations from land stations, ships, drifting buoys, reports from radiosondes and pilot balloons launched from land stations and ships, drop sondes, aircraft reports, wind profiler, and METAR airport weather reports.
The ERA-Interim reanalysis quite successfully reproduces most of the atmospheric parameters, including air temperature and precipitation (Dee et al. 2011;Donat et al. 2014;Sun et al. 2018). Over the Yamal Peninsula, observational data of weather stations Marresale and Novy Port are assimilated in the ERA-Interim. In this study we used daily 2 meter temperature and daily precipitation to compare station observations with reanalysis data. As observational data, we used the Marresale weather station (69.72°N, 66.80°E) records from RIH-MI-WDC archive (Bulygina et al. 2008). The data sampling of this station is representative because missing data amounts to not more than 10%. Selected reanalysis grid point for comparison has coordinates 69.75°N, 66.75°E.
There is a consistent overestimation of the daily 2 meter temperature in winter, spring and autumn and an underestimation in the summer by the ERA-Interim. The grid cell, which contains the Marresale weather station, includes both land and sea. Therefore, the temperature in winter, spring and autumn according to reanalysis is higher than according to station observations, and in summer, it is lower. The standard deviation of the difference between daily 2-meter temperature data series from ERA-Interim and from the station is 3.2°С. This is the effect of the sea area included in this reanalysis grid cell.
At the same time, ERA-Interim reanalysis underestimates the daily precipitation data compared to Marresale station. A decrease of precipitation compared with observational data is noted in most reanalyses and grid archives for the entire globe. Also, the reanalysis significantly underestimates the values of extreme precipitation, one of the reasons for that being the local character of the extreme values (Donat et al. 2014;Matveeva et al. 2015), while reanalysis data are averaged over the cell. The standard deviation of a series of daily precipitation is 1.8 mm.
Comparison of hydrological parameters (including runoff ) from ERA-Interim reanalysis data with observational data showed significant inconsistency (Betts et al. 2009). To increase the accuracy of gully erosion calculation, we apply our own hydrological model.

Hydrological model
Any erosion modelling begins with water flow calculations. In the conditions of the Arctic climate with deep permafrost and high density of linear erosion features hydrological modelling requires specific characteristics described further. Two main sources of water flow are typical for this environment: snow thaw during the spring and rainfall, mostly during the summer. The runoff depth Hrunoff for the period of snow thaw was calculated with the models of Komarov (Komarov et al. 1969), Vinogradov (1988), Vinogradov et al. (2014, and Gelfan (1999), using the maximum depth of snow fall Hsnow (water equivalent), precipitation P, air temperature T, and evaporation E from reanalysis. Snow thaw rate was calculated with the air temperature and snow thaw coefficient k sm . Runoff occurs when the content of water in snow exceeds the maximum potential water content, which was calculated with current snow density. Snow thaw period ends when Hsnow=0.
The initial snow pack depth on the catchment is irregular. In the natural conditions, snow driven by wind during the winter accumulates in depressions (river valleys and gullies), at the base of steep slopes and in dense vegetation. In the conditions of human impact, snow accumulation around the buildings and roads becomes significant. The measurements on the gullied territory of the west-central Yamal Peninsula (Bobrovitskaya et al. 1999) show the typical distribution of snow depth/mean snow depth ratio (Nratio) at the beginning of snow thaw period. The distribution is bimodal, the first mode describing the snow pack variability in the areas near the water divide and on gentle slopes, and the second one in the gullies. The relative weight of these two modes depends on the areas proportion of gentle slopes and gully cuts in the area. To take into account the uneven thickness of the snow cover, the calculations of the runoff depth Hrunoff for each time period were performed with different values of Nratio, weighted with their probabilities and then integrated to get the resulting value. The runoff depth for the period of summer rains was assumed to be equal to rainfall depth P(t) from reanalysis. Evaporation E for the periods with rainfall was assumed to be 0, and the losses were calculated as the amount of water required to fill depressions on the catchment and water evaporated from them during the periods between rainfalls using the formula of Popov (1956).
The sensitivity analysis of the hydrological model shows that the most important are absolute values of input climatic characteristics, type of snow pack depth distribution, and snow thaw coefficient k sm . The climatic characteristics were taken from reanalysis. They are quite valid for air temperature, but differ at random from measured values by up to -50% to +100% for snow pack depth at the beginning of snow thaw period and for precipitation depth. The coefficient of snow thaw k sm was calibrated using the measurements over the spring period of 1993, when both snow pack depth distribution and runoff were measured (Bobrovitskaya et al. 1999). For k sm =0.75 10 -4 the measured and calculated runoff depth are well correlated (Fig. 6A). Runoff sequence during the spring 1992, calculated with this calibrated value of k sm , have the pattern similar to measurements (Fig. 6B). The measured snow depth (in water equivalent) at the beginning of snow thaw period in 1992 was nearly two-fold that from reanalysis, and the measured and calculated runoff depths are well correlated only when snow depths from reanalysis were corrected according to measurements. Without such correction, the calculated maximum runoff depth is considerably underestimated (65 instead of 80 mm), and the date of this calculated maximum was one day earlier than the real one.

Erosion model
Gully thermoerosion and erosion model GULTEM was previously described in detail in (Sidorchuk 2015). Therefore, only the major formulae are listed below. According to the model, two main processes occur along each flowline, as in nature: 1) formation of erosion cut in the topsoil or at the gully bottom by flowing water during snowmelt or a rainstorm event; 2) transformation of gully walls by shallow landslides during the period between adjacent water flow events.
In the first process, the rate of erosion is linearly correlated with the product of bed shear stress τ=gρDS and the mean flow velocity U: Here S is gully bottom mean slope, q is specific discharge, i.e. the ratio of discharge and flow width Q/W, g is acceleration due to gravity, ρ is water density and k E is erosion coefficient, which depends on soil properties. If flow velocity is less than critical value for erosion initiation U cr , then erosion is negligible.
In the case of gully erosion in the soil with permafrost, another type of bed lowering develops -the combined action of thawing of frozen soil by running water with positive temperature and mechanical removal of thawed layer -socalled thermoerosion. Thermoerosion rate is equal to the rate of soil thawing and linearly correlated with water temperature T 0 C: The linear relationship of thermoerosion rate with water temperature, as well as k TE values for different types of frozen soil were obtained in the course of laboratory and field experiments (see details in Sidorchuk 2015). The same type of relationship follows from solution of Stefan's problem for the case of thermal erosion (Yershov et al. 1982).
The erosion rate of gully banks can be estimated only in a first approximation as a function of the rate of gully bed erosion: Here V is the lateral flow velocity. For k b estimates see Eq. 42 and 43 in Sidorchuk 2015. Together, erosion and thermoerosion produce gully volume V 0 .
The second process, rapid sliding and slumping following the incision, leads to formation of stable straight slopes of the gully with an inclination equal to the angle of repose φ for a given soil type. This process transforms the erosion cut of volume V 0 into trapezium with bottom width W b , depth and top width The altitudes Z of gully longitudinal profile are transformed in time t along flowline length x according to the following partial differential equation for the process of erosion and according to ordinary differential equation for the case of thermoerosion If the rate of thermoerosion is less than erosion rate, the thaw layer is removed, and Eq. 7 is used for calculating the gully bed lowering. If the rate of thermoerosion (i.e., the rate of the thermal front shifting in soil) is greater than erosion rate, a thaw layer is formed, and calculation with Eq. 6 is used in the model.

Fig. 6. Measured discharges and calculated runoff depths during the spring of 1993 (A) and 1992 (B)
(1) (2) (3) (4) The input morphological and hydrological characteristics in the model are: longtitudinal profile Z(x) of initial slope, discharge Q(x,t), flow width W(x,t), and water temperature T(x,t). As the gully basin is relatively small, discharge can be calculated as product of runoff depth and contributing catchment area Hrunoff(t)F(x). The flow width, as well as other flow hydraulic characteristics, are calculated using the empirical relationships with the discharge (see Eq. 55 and 56 in Sidorchuk 2015).
The coefficients and parameters used in the model are k E , k TE , k b , U cr , and φ. Some of these values are used for the model calibrations, and others must be known. Usually k E and/or U cr are used for the calibration. Field and laboratory experiments were used for k TE estimations (see Table 6 in Sidorchuk 2015) and k TE =0.000025 m/( o C s) was taken for the loam deposits of the west-central Yamal. The measurements of bank erosion in the gullies show that the coefficient of bank erosion depends on the ratio between gully flow width W and gully bottom width W b , as well as on the flow depth in the gully D (see Eq. 42 and 43 in Sidorchuk 2015). The angle of repose φ is estimated from the stable gully sections or calculated with the straight slope stability model (see Eq. 46 and 47 in Sidorchuk 2015). The angle of repose φ for the gullies investigated in 1990-97 was about 30 o (the mean value for 11 stable cross-sections).
Calibration procedure was two-staged. At the first stage, the critical velocity of erosion initiation was calibrated according to the information about gully length growth, as this morphological characteristic mostly depends on critical velocity in given hydrological conditions. Two main scenarios were used: 1) initial catchment topography, which existed prior to the gully incision, with undisturbed natural vegetation cover; 2) catchment topography in the year 1990 with vegetation cover completely disturbed due to human impact. Critical velocity U cr was 0.8 and 0.12 m/s, respectively. At the second stage, the optimal erodibility coefficient k E was found according to the information about gully mean depth change during the calibration period 1990-95 (Fig. 7); optimal erodibility coefficient k E value was 0.0025 1/s. Finally, with the values of all coefficients and parameters k E =0.0025 1/m, k TE =0.000025 m/( o C s), U cr =0.12 m/s, and φ=30 o the gullies evolution calculated for the period 1986-2016 was compared with the surface photos of 2007 and space images of 2016 (see Fig. 5).

RESULTS AND DISCUSSION
Before discussing the results it is useful to consider two questions: 1) the interconnection of thermokarst and erosion processes in formation of the gullies in the regions with permafrost, and 2) the use of morphological analogs as hydrological ones.
The first question was raised by Kosov and Konstantinova (1970) during their investigation of erosion processes in the Arctic. Among other types they distinguished so-called thermokarst -erosion gullies, which were formed mostly due to in situ melting of polygonal ice wedges without significant influence of running water from the catchment. These features were rather small and narrow, and their pattern resembles the relief of ground hummocks of conical shape (so-called «baydjarakhs»). This type of gullies was only briefly mentioned in the subsequent works of the same authors (see Yershov et al., 1982); instead, after more extensive investigations they described typical thermo-erosion gullies formed by running water (see, for example, Konstantinova 1973). Special works devoted to gully erosion processes in the Arctic (Poznanin 2012), as well as general works in geocryology (e.g., Yershov and Williams 2004) did not discuss the formation of such gullies without running water, and this hypothesis was nearly forgotten.
Interest in the formation of such ravine-like features formed in the absence of water has been renewed due to planetary studies. Numerous attempts to find terrestrial analogs for the gullies on Mars (see, for example, review by Conway et al. 2019) showed that water-lacking processes (including ground ice melting) can produce only small gully-like features, which cannot be compared with large gully nets formed by water erosion. A general degradation of permafrost on the East European Plain with former periglacial conditions left numerous remnants of thermokarst features in recent relief morphology (Velichko 1965). However, that is the relief of small flat hummocks and shallow pits, but not the linear erosion pattern. In the modern conditions in the Arctic, some small tributaries not more than 10-20 m long can be formed as the result of ice wedges melting, but only running water prevents them from filling by thaw ground slumps from the slopes. The general shape of gully network can follow the system of cracks in frozen ground. Often in this case, not thermokarst initiates erosion, but erosion and thermoerosion by running water leads to ice wedges melting in the gully cuts (Godin et al. 2012). In some cases, an opposite process was described (Levy et al. 2008): filling of the former gullies by ice wedges. Certainly, the processes of ground ice melting influence the morphology of gullies in permafrost conditions, but this influence remains secondary to the main process of erosion and thermoerosion by water from the catchment even in the Antarctic (Dickson et al. 2019).
The second question arose after Horton (1945) formulated the basic laws of channel network morphology. Gregory (1966) came to conclusion that for hydrological indication it is preferable to use a channel network, which includes both permanent and temporary watercourses. He was one of the first researchers to appreciate the structure of erosion networks as a source of paleohydrological information. Gregory pointed out that the presence of dry valleys not occupied by permanent watercourses may reflect a significant decrease in water flow and pointed out that paleohydrological estimates based on changes in the shape of river channels and cross sections can be supplemented by studying the patterns of fluvial systems.
Gartsman (1968) found a linear relationship between the total length of watercourses ΣL in a system of order K and the mean output discharge Q of this system. He introduced a hydro-morphological coefficient γ Q (HMC) as an indicator of the length of the river network necessary for the formation of a unit discharge The hydro-morphological coefficient represents a way of using morphological analogs as hydrological ones, since the networks of channels of the same density are hydrologically similar if certain geomorphological and lithological characteristics coincide. This opens a possibility to use a morpho-paleohydrological analogy, that is, correspondence between the structure of modern channel networks and their hydrological characteristics can be applied to estimate the hydrological characteristics of the ancient channel systems by their structure.
The relationship between the annual runoff depth (mm) and the averaged density of the dry valley-river network (km/ km 2 ) was established for nine large lowland regions with tundra and sparse northern taiga vegetation and permafrost (Sidorchuk et al. 2018). This relationship is linear, with the coefficient of determination R 2 =0.85 (Fig. 8). A relatively high network density of dry valleys and rivers (1 km/km 2 ) is formed by the annual runoff depth of ~270 mm, which indicates a high erosion efficiency of water flow in such landscapes. The density of the network of dry valleys and rivers K in the Khoper River basin (a tributary of the Don River) is 0.73 km/km 2 . This erosion landscape was formed during the powerful surface runoff event at the beginning of the degradation of the last glaciation 16-19 ka ago. The closest modern morphological and hydrological analogs of this landscape are the Northern Siberian Lowland (K=0.7 km/km 2 ) and Yamal Peninsula (K=0.77 km/km 2 ), where drainage networks were formed in the Holocene. For further calculations, we used the hydrological data on the Yamal Peninsula, as the region is much better investigated than the Northern Siberian Lowland.

Estimations of the past periglacial erosion rate
After successful calibration of GULTEM for the gullies in the west-central Yamal, we used this model for palaeohydrological reconstructions. Calculations were carried out for the development of the erosion network within the catchment area of the Perepolye dry valley in the Khoper River basin. The initial relief of the catchment was reconstructed by complete "filling" the existing erosion forms, while the modern structure of the drainage network was preserved in the calculations. The main task was to estimate the total length of the gully system and the maximum length of the main trunk for two main periods of catastrophic surface runoff -the Late Moscow time (MIS 6) and the Late Valdai time (MIS 2). As up to 90% of a gully length forms during the first 5% of its lifetime, a 30-year simulation is sufficient for the purpose. The length of the main trunk of the active paleogully was calculated with critical velocity of erosion initiation U cr varying within the range 0.2-1.6 m/s and the ratio M/M Yamal from 0.2 to 4.1 due to increased or decreased surface runoff depth M (Fig.  9). The recent mean annual runoff depth M Yamal at the reanalysis grid point 70.25°N, 68.25°E on the Yamal Peninsula is 340 mm, the maximum runoff depth for six-hour period is 50 mm for the snow thaw and 24 mm for rainfall period. The main trunk length of the modern Perepolye dry valley system inherited from the Late Valdai time is 6400 m. The critical velocity of erosion initiation U cr corresponding to this value is 0.87 m/c, which is more than the value obtained for the Yamal Peninsula with natural (undisturbed) vegetation. Since the hydrological conditions in the Yamal Peninsula (M Yamal ) are not completely analogues to the conditions in the Perepolye catchment in the Late Valdai, an uncertainty of the M/M Yamal ratio of ±20% was included in the calculations. This uncertainty interval leads to variability in critical velocity of erosion initiation U cr within the range 0.8-0.9 m/s. If considering the difficulties of estimating this characteristic for the natural objects, such a range is not especially broad.
The value of U cr =0.87 m/s was then used to reconstruct the hydrological and erosion conditions in the late Moscow time (MIS 6). The system of shallow hollows (that were gullies at the end of the Moscow Glaciation) is denser and longer than the system of dry valleys. The length of the main trunk of an active paleogully was then 8400 m. The calculations with GULTEM showed a high annual and maximum depth of the surface runoff (see Fig. 9). The calculated value of M (1090 mm), which corresponds to the network of hollows in the case of U cr =0.87 m/s, is almost 3.3 times greater than M Yamal . The uncertainty of this estimate of M/MYamal is very high and is in the range 2.3-4.1. Such high errors are explained by the very weak relationship between the surface runoff and the growth of the gully when its top is close to the watershed.

CONCLUSION
The same type of erosion landscape is characteristic of the Arctic regions, especially for the Yamal, Gydan and Tazovsky peninsulas with active gullies, and for various regions of the East European Plain, where ancient active gullies were later transformed into stable dry valleys. In the both cases, this type of landscape develops due to erosion in periglacial conditions. Thus, the Arctic regions can be used as modern analogs to reconstruct the conditions for the development of periglacial erosion features in regions with a modern temperate climate. For these purposes, gully erosion and thermoerosion model GULTEM was verified and calibrated The most suitable characteristics of climate and vegetation cover, which can explain the structure and density of dry valley systems in the East European Plain were then estimated with this model. The main input characteristics, which control calculation results are critical velocity of erosion initiation U cr and surface runoff M. These characteristics were varied in numerical experiments to obtain the best similarity with the morphology of existing dry valleys. The topography of the modern dry valley and hollow system in the basin of the Perepolye dry valley (the Khoper River basin) was used to calculate the length of the active paleogullies. The length of the main trunk of this erosion system was 8400 m in the end of the Moscow Glaciation (MIS 6) and 6400 m in the end of the Late Valdai Glaciation (MIS 2). The critical velocity of erosion initiation U cr corresponding to these values is within the range of 0.8-0.9 m/s for both time intervals. For the first period, the surface runoff, which fits the network of hollows, is almost 3.3 times greater than the recent one in the Yamal Peninsula, or more than six-fold the modern value in the Khoper basin. For the second period, the mean annual surface runoff, which could produce these large gullies, was close to the present annual surface runoff in the Yamal Peninsula, approximately 330 mm. This is about twice the modern runoff on this catchment. The uncertainty of these estimates is relatively low for the middle parts of the catchments and increases significantly near the watershed.