Thermo‐erosional valleys in Siberian ice‐rich permafrost

Thermal erosion is a major mechanism of permafrost degradation, resulting in characteristic landforms. We inventory thermo‐erosional valleys in ice‐rich coastal lowlands adjacent to the Siberian Laptev Sea based on remote sensing, Geographic Information System (GIS), and field investigations for a first regional assessment of their spatial distribution and characteristics. Three study areas with similar geological (Yedoma Ice Complex) but diverse geomorphological conditions vary in valley areal extent, incision depth, and branching geometry. The most extensive valley networks are incised deeply (up to 35 m) into the broad inclined lowland around Mamontov Klyk. The flat, low‐lying plain forming the Buor Khaya Peninsula is more degraded by thermokarst and characterized by long valleys of lower depth with short tributaries. Small, isolated Yedoma Ice Complex remnants in the Lena River Delta predominantly exhibit shorter but deep valleys. Based on these hydrographical network and topography assessments, we discuss geomorphological and hydrological connections to erosion processes. Relative catchment size along with regional slope interact with other Holocene relief‐forming processes such as thermokarst and neotectonics. Our findings suggest that thermo‐erosional valleys are prominent, hitherto overlooked permafrost degradation landforms that add to impacts on biogeochemical cycling, sediment transport, and hydrology in the degrading Siberian Yedoma Ice Complex.


| INTRODUCTION
Climate warming in the Arctic is occurring at a much faster rate than the global average, 1 impacting polar permafrost regions. Permafrost warming and erosion of permafrost deposits have been reported throughout the northern high latitudes. [2][3][4] Thermokarst and thermal erosion are two major processes of ice-rich permafrost degradation in periglacial landscapes. They result in thawing of permafrost-stored organic carbon, which then can decompose and be released to the atmosphere and the hydrological system. [5][6][7][8] Landscapes affected by thermokarst and thermal erosion are estimated to contain 330 Pg carbon in the upper 3 m subsurface, constituting 30% of the total 0-3 m permafrost region soil organic carbon storage (1,073 Pg C), highlighting their importance for the global carbon cycle in a rapidly warming Arctic. 9,10 The two processes may also substantially alter the water and energy balances of affected regions. 11,12 Thermal erosion is defined as the erosion of ice-bearing permafrost by the combined thermal and mechanical action of moving water. 13 Whereas thermokarst is an in situ process including melting of ground ice followed by surface subsidence but without hydraulic transport of earth materials, thermal erosion is a dynamic process involving the wearing away by thermal means (i.e., the melting of ice) and by mechanical means (i.e., hydraulic transport). 13 Two types of thermal erosion can be distinguished: linear or vertical thermal erosion, which acts to depth, and lateral thermal erosion, which acts sideways. 14,15 Thermokarst processes occur mostly in flat lowland terrain with low hydraulic gradient; the resulting characteristic landforms are thermokarst lakes, thermokarst depressions (alasses), and thermokarst mounds. Thermal erosion takes place along river banks [16][17][18][19][20] and coastlines, 21,22 at the shores of lakes, 23,24 but also in lowlands with abundant ground ice and sufficient gradients to allow for channelized surface water flow. 25 Here, it can result in the formation of thermoerosional gullies [26][27][28][29] or even large thermo-erosional valleys and valley systems. Thermal erosion, thermokarst, and the landforms resulting from these processes interact with each other, as thermo-erosional gullies and valleys can supply water to thermokarst lakes and basins 30,31 and enlarge thermokarst depressions. 32 At the same time, they may also inhibit thermokarst activity by drainage of flat uplands and breaching of thermokarst lakes. [33][34][35][36] In contrast to thermokarst lakes, which have been investigated in numerous studies, for example as sources of carbon release to the atmosphere, 5,6 with respect to spatial distribution patterns, 37 or as indicators of a changing water balance in permafrost regions, 38,39 only a few studies are available on thermo-erosional gullies and valleys or discuss the interaction of thermal erosion with thermokarst under differing regional relief conditions. Thermo-erosional valleys have been described in the Lena River Delta 33 and mapped in two regions of the East Siberian coastal lowlands in the context of an overall quantification of thermokarst-affected terrain types. 30,40 Other studies have reported rapid formation and growth of thermo-erosional gullies due to thawing permafrost. [26][27][28][29]41 Gullies and valleys can deeply erode underlying deposits in a short time and therefore increase sediment and nutrient delivery to rivers, lakes, and the sea. 32,41,42 They also may reconfigure drainage networks, thereby leading to large changes in runoff volume and timing. 43 Thermo-erosional gullies and valleys are common in Arctic landscapes, act as important snow accumulation areas, and alter the water, sediment, and organic matter transport from inland permafrost to receiving waters at the local scale. Yet, there are no systematic studies available that classify and quantify thermo-erosional valleys and their impacts at the Arctic scale or that analyze differences in valley network characteristics across different permafrost regions.
As a first step to fill this gap, we provide here an overview of thermo-erosional valleys and their role in the degradation of ice-rich permafrost in the eastern Siberian Arctic. The specific objectives of this study are: (a) to describe the morphometry and spatial distribution of thermo-erosional valleys across three Laptev Sea coastal lowland regions based on remote sensing, geoinformation, and field data; and (b) to relate the identified spatial patterns and morphologies to the topographical and cryolithological settings of the study areas and their evolution through the Holocene. In particular, we address the following questions: How abundant are thermo-erosional landforms in the study region? Which drainage patterns and valley types occur? How do the valleys and their networks differ between the study areas and why? and we propose a conceptual model for valley evolution in the study region during the Holocene.

| REGIONAL SETTING
This study focuses on three sites bordering the Siberian Laptev Sea coast ( Figure 1). The stratigraphy in all three study areas is comparable and consists of late Pleistocene polygenetic Yedoma Ice Complex deposits underlain by middle to late Pleistocene sandy deposits and covered by various Holocene deposits. [44][45][46][47][48][49][50][51][52] The Yedoma Ice Complex is 10-40 m thick and consists of ice-supersaturated silty to sandy sediments and buried cryosols. Large syngenetic ice wedges of up to several tens of meters thickness truncate the sediment column and contribute to a total ground-ice volume of up to 80%. All three sites are situated in the Arctic tundra zone 53 and in the continuous permafrost zone with several hundreds of meters permafrost thickness and mean annual ground temperatures from −9 to < −11 C. 15 The western site around Cape Mamontov Klyk belongs to the Olenyok-Anabar coastal lowland. It is bordered by the Pronchishchev Ridge to the south and is inclined slightly in the NNE direction towards the Laptev Sea. 30 The central site in the Lena River Delta is composed of insular Yedoma remnants and, together with a few locations of bedrock outcrops, they constitute the third geomorphological main terrace of the delta, which accounts for 5.9% of the total delta area of 29,000 km 2 . The second main terrace consists of fluvial sandy deposits and covers 21% of the Lena Delta area. The third and second terraces represent late Pleistocene non-deltaic units, whereas the Holocene deltaic deposits are considered the first geomorphological main terrace of the Lena Delta. 33,34,45,47,54 The eastern site, the Buor Khaya Peninsula, is the westernmost part of the Yana-Indigirka coastal lowland. In contrast to the other sites, the base of the Yedoma Ice Complex deposits is below modern sea level (b.s.l.) here and is currently not exposed. Coring about 800 m off the western shore indicates its location at about 4 m b.s.l. followed by Pleistocene sands that extend down to at least 52 m b.s.l. 55,56 Holocene deposits of 1-2 m thickness cover the Yedoma Ice Complex deposits. Holocene deposits are also found in thermokarst depressions and in river and thermo-erosional valleys. They consist of peat, silty to sandy sediments with high organic matter and ice content, and lacustrine silty sediments. 57 The high ice content and large thickness (10-40 m) of the Yedoma Ice Complex deposits cause substantial terrain subsidence due to thawing and has provided the conditions for widespread surface changes from thermokarst and thermal erosion since the Lateglacial/Holocene transition about 12-10 ka BP. 35 Valley density was then calculated as the sum of rivers, streams, and intermittent streams outside valley floodplains and the valley floodplain centerlines divided by the total area of the study area. Figure S1 provides an illustration of the difference between drainage density and valley density and their calculation. Terrain analyses and the extraction of valley transverse profiles were performed using the ArcticDEM, which is based on stereophotogrammetrically derived elevations from sub-meterresolution satellite imagery. 67 We used mosaic tiles of Release 6, v2.0 with a 5-m resolution that were averaged from multiple 2-mresolution strip DEMs and provide a consistent data source over large areas. The suitability of the ArcticDEM to extract thermo-erosional valley profile characteristics has been confirmed by comparing GNSSmeasured valley profiles with those extracted from the ArcticDEM at the same locations ( Figure S2).

| RESULTS
The summary characteristics of the three study sites are given in Table 2. Cape Mamontov Klyk is extensively covered with valleys and streams ( Figure 2) and therefore has the highest drainage density   In all three study areas, straight, short gullies are often radially located around alasses and thermokarst lakes ( Figure 3f). They start on the Yedoma uplands, cut into the slopes, and end abruptly at the foot of the alas slope or at the lake level.  (Table 2).

| Categorization of thermo-erosional valleys
Based on the results of spatial analyses and field observations ( Figure 5) we categorize the thermo-erosional landforms in the study area into the following types (Table 3): 1. Short, straight gullies around thermokarst lakes and alasses cut down to 10 m deep into the slopes of the thermokarst features and mostly follow the same gradient as the slope, as indicated by the uniform width of the bank slopes along the gullies (Table 3a, Figure 4a).
They are densely covered with more vital and higher growing wet sedge tundra vegetation than the surrounding slope sections.
2. In alasses, small streams often occur as drainage pathways and connect residual and secondary thermokarst lakes on the alas floor with the hydrological system outside the alasses. In cases of catastrophic drainage of large thermokarst lakes, the draining water cuts outflow channels into the unfrozen soft sediments at the lake bottom that persist as valleys on the floor of the newly formed alas and that can reach 20 m depth (Table 3b, Figure 4b).

| Data accuracy
The geo-dataset of thermo-erosional landforms analyzed in this study was compiled using existing datasets that had been mapped based on topographical maps and satellite imagery with differences in spatial resolution and acquisition times between the study areas (Table 1).
Due to this heterogeneity of the dataset, we cannot entirely exclude differences in the degree of mapping completeness and hence the representation of the lateral extent of thermo-erosional landforms between the three study areas. In this case, these differences would have propagated into the results of the analyses of spatial metrics for a particular area, such as total stream length or drainage density (Table 2). Regarding the subsequent comparison between the study areas, in an unfavorable case this could have added up to inconsistencies with the potential to decrease the validity of this cross-area comparison. However, we consider these effects to be minimal for the following reasons.
The best ground resolution of the satellite data used as a basis for the mapping of thermo-erosional landforms was 10 m for the Cape Mamontov Klyk and Lena Delta areas and 6.5 m for the Buor Khaya Peninsula. This ground resolution sufficiently allowed for the detection of all major thermo-erosional landforms analyzed in this study.
Even though landforms of the smaller categories can be <10 m wide

| Conceptual models of valley formation and evolution
In the following, we propose conceptual models of the formation and evolution of thermo-erosional landforms in the study region. In the first part, we focus on the individual valley categories that we distinguished in section 4.3 and which occur in all three study areas. In the second part, we discuss the factors and conditions that led to regional differences in the evolution of the study areas and their valley networks.

| Valley categories
The short, straight gullies on the slopes of thermokarst lakes and alasses ( Figure 3f, Table 3a)  The V-shaped valleys (Table 3d), U-shaped valleys (Table 3e), and valleys of permanent streams (Table 3f)  source in a mountain ridge with exposed bedrock and may therefore have a braided stream origin. Once coarse sediment supply ceases, braiding intensity decreases, and distinct channels form. 75 In the case of abrupt elevation changes, for example at the junctions of streams or valleys of different orders, valley transverse profiles can change from one type to another. Figure 5b shows such an example, where a steep V-shaped ravine is retrograding into a small and shallow U-shaped valley.

| Evolution of the valley networks
The study areas differ greatly in terms of drainage patterns and densities as well as in the distribution of the different valley categories and thermokarst landforms, even though they all represent Yedoma Ice Complex settings. This implies that factors and conditions governing valley evolution and Yedoma Ice Complex degradation must have affected the study areas differently. Here we propose a conceptual model including regional landscape evolution and site-specific development that led to the present-day valley network characteristics ( Figure 6).
In the study region, development of the modern hydrological system began with the Lateglacial/early Holocene transition period, when accumulation of the Yedoma Ice Complex deposits ceased and the regional climate shifted to warmer and wetter conditions, promoting the activation of rivers and widespread permafrost degradation. 31,34,[46][47][48]58 The coastline of the Laptev Sea was located several hundred kilometers further north during that period, 76 and extensive ice-rich permafrost lowlands with low gradients provided favorable conditions for thermokarst (Figure 6a).
Thermokarst processes slowed down substantially with mid-to late Holocene cooling 35,58,77 (Figure 6c). In the Lena Delta, the Yedoma surface was cut off from the mountain ranges to the south as well as from the accumulation plains to the north by large deltaic channels and has been eroded into small disconnected remnants. They are now elevated up to 66 m above deltaic streams and floodplains, with higher relief in the western than in the eastern delta (Figure 6d). 44,47 The Buor Khaya Peninsula was strongly affected by thermokarst processes, which degraded a large proportion of the Yedoma Ice Complex deposits (85-90% of total area), 80 leaving only remnants of undisturbed Yedoma uplands on a low-elevation plain (Figure 6e). We argue that this marked difference between the western and eastern region reflects neotectonic activity (Figure 6b).
The study region is seismotectonically very active, because it is located at the zone of intense deformation between the North American and Eurasian plates, 81  In a warming Arctic with active-layer deepening 94 and precipitation regime shifts, 95 permafrost degradation is expected to increase.
The consequences of this increase include more thermal erosion and thermokarst, intensified organic matter cycling, and release as well as shifts in the seasonality and intensity of the hydrological regime. 10,96 In the landscapes analyzed in this study, this will lead to an expansion