Geomorphological and Climatic Drivers of Thermokarst Lake Area Increase Trend (1999–2018) in the Kolyma Lowland Yedoma Region, North-Eastern Siberia

Abstract

Thermokarst lakes are widespread in Arctic lowlands. Under a warming climate, landscapes with highly ice-rich Yedoma Ice Complex (IC) deposits are particularly vulnerable, and thermokarst lake area dynamics serve as an indicator for their response to climate change. We conducted lake change trend analysis for a 44,500 km² region of the Kolyma Lowland using Landsat imagery in conjunction with TanDEM-X digital elevation model and Quaternary Geology map data. We delineated yedoma–alas relief types with different yedoma fractions, serving as a base for geospatial analysis of lake area dynamics. We quantified lake changes over the 1999–2018 period using machine-learning-based classification of robust trends of multi-spectral indices of Landsat data and object-based long-term lake detection. We analyzed the lake area dynamics separately for 1999–2013 and 1999–2018 periods, including the most recent five years that were characterized by very high precipitation. Comparison of drained lake basin area with thermokarst lake extents reveal the overall limnicity decrease by 80% during the Holocene. Current climate warming and wetting in the region led to a lake area increase by 0.89% for the 1999–2013 period and an increase by 4.15% for the 1999–2018 period. We analyzed geomorphological factors impacting modern lake area changes for both periods such as lake size, elevation, and yedoma–alas relief type. We detected a lake area expansion trend in high yedoma fraction areas indicating ongoing Yedoma IC degradation by lake thermokarst. Our concept of differentiating yedoma–alas relief types helps to characterize landscape-scale lake area changes and could potentially be applied for refined assessments of greenhouse gas emissions in Yedoma regions. Comprehensive geomorphological inventories of Yedoma regions using geospatial data provide a better understanding of the extent of thermokarst processes during the Holocene and the pre-conditioning of modern thermokarst lake area dynamics.

1. Introduction

Thermokarst lakes are widespread in Eurasian and North American permafrost lowlands and their formation occurs in regions formed by ground ice-rich unconsolidated sediments [1]. Olefeldt et al. [2] estimate coverage by thermokarst landscapes of about 20% of the permafrost zone, where the thermokarst lake area is on average about 7% and in places reaches up to 30–40% [3,4]. These lakes are of great importance for the sensitive environmental systems in the northern high latitudes not only due to their large spatial coverage in Arctic coastal lowlands, and their hydrological and habitat functions [5,6,7], but also as a landscape element important for carbon sequestration and an incubator for subsequent methane production from thawing permafrost [2,8,9,10]. Modern climate warming leads to widespread permafrost degradation throughout the Arctic, which is reflected in active layer deepening, and thermokarst and thermoerosion processes activation (e.g., [11,12,13,14,15,16]). Studying thermokarst lakes provides insights into permafrost dynamics in response to climate change. The dynamic behavior of thermokarst lake area is a well-known phenomenon that has been reported from continuous and discontinuous permafrost regions across the Arctic (e.g., [17,18,19,20,21,22,23]).

Contrasting patterns of lake growth, shrinkage, and drainage can be observed at the same time and in the same area. However, both processes can be related to thermal permafrost degradation and erosion of permafrost deposits. The heat capacity of lake water favors thaw along the lake shores and ground subsidence under the lake bottom, leading to lake expansion [18,24]. In surrounding areas, accelerated coastal erosion and thermo-erosional development along deeply incised valleys cause lake tapping and subsequent drainage and lake disappearance [25,26,27]. In the thinner permafrost of the discontinuous permafrost zone, lakes may also drain due to thawing through the permafrost layer [28]. Although both lake growth and drainage effects may level out on average without major effects on the total net lake area change in some region [22,29], a general trend of decreasing lake area is observed in many regions of the discontinuous and continuous permafrost zones [23,30,31]. In the continuous permafrost zone, decreasing thermokarst lake coverage was detected in Subarctic Alaska [22,32,33], in the European part of Russia [34,35], and Western Siberia [21]. Increasing lake areas were observed in the continuous permafrost zone of Northwest Canada on the Tuktoyaktuk Peninsula [17], in the Yedoma region of Central Yakutia [22,36,37,38], and on the Qinghai-Tibetan Plateau [39]. Studies in Arctic Yedoma regions show a general trend of lake area decrease in Alaska [18] and the Yakutian coastal lowlands [21,22,40,41].

Among the regions with the highest abundance of thermokarst lakes are the territories with Yedoma Ice Complex (IC) deposits distribution [2], which are widespread in previously unglaciated lowlands of north-eastern Eurasia, Alaska and Northwest Canada [1]. Yedoma IC deposits are characterized by a high content of buried well-preserved soil organic matter (OM) [42,43,44], which has a potential of fast decomposition upon permafrost thaw and thus the formation of CH₄, CO₂, and other greenhouse gases [45,46,47,48]. Because of the very high ice content ranging from 60 to 90% due to the presence of large syngenetic polygonal ice wedges [27,49,50,51], the Yedoma IC is generally very vulnerable to climate warming, which is reflected in its drastically diminished areal extent since the onset of the late Pleistocene–early Holocene [52,53,54]. To assess whether this tendency is intensifying under modern climate warming, thermokarst lake area dynamics serve as a proxy to quantify the current rate of permafrost degradation processes in Yedoma regions.

Landsat satellite images with 30 m spatial resolution are available for the last decades and allow to estimate long-term landscape changes (e.g., [22,31]). Nitze et al. [22] developed a highly automated and hybrid approach based on Landsat data to analyze lake area changes. This approach has been proven for effectively monitoring in different Arctic regions and is designed to detect long-term lake change trends over periods longer than 10 years. Furthermore, this method reduces the influence of short-term fluctuations and seasonal variations. However, to put modern lake area dynamics into perspective, it is also necessary to incorporate thermokarst process development and its geomorphological conditions into the analysis for regional scales. Research efforts more frequently reflect on the impact of in situ cryolithilogical and relief conditions on thermokarst lake area changes [28,38,55,56,57,58,59,60,61,62,63]. Beside these locally highly detailed investigations, there is a gap of larger regional scale remote sensing-based monitoring efforts that take geological and relief conditions and their influence on lake area dynamics into consideration.

Analyzing hydrometeorological influences along with investigations of lake dynamics are an established standard approach. In some regions, the role of annual air temperature [39,59], precipitation [21,29,64,65], and summer weather conditions [20,34,66] were evaluated. Our study builds on this approach by combining detailed landscape-scale lake area change detection based on moderate and high-resolution remote sensing with available climate data analysis.

The Kolyma Lowland tundra is one of the largest Yedoma regions in the Arctic with a high abundance of thermokarst lakes, whose dynamics and influencing factors are not well understood. The objective of our study is to estimate modern thermokarst lake area dynamics of the Kolyma Lowland tundra Yedoma region in North-Eastern Siberia and to analyze the endogenous geological and geomorphological, and exogenous hydroclimatic factors impacting the spatial pattern of lake area changes.

2. Study Region

The study region is in North-Eastern Siberia within the Kolyma Lowland tundra zone with an area of about 44,500 km² (Figure 1). The territory is bordered by the East Siberian Sea to the north and the Khallerchin tundra composed of sandy deposits to the east [67,68]. The western part is bounded by the Suor-Uyata Ridge and the western boundary of the Alazeya river basin. The southern boundary of the tundra zone was delineated manually using Landsat images for the 2001–2010 years based on the spectral differences for forested and non-forested areas.

The region is in an area of steady tectonic subsidence during the Oligocene-Pleistocene, which contributed to the accumulation of loose sediments with a thickness reaching up to 500 m [69,70]. During the Pleistocene, the region remained unglaciated [71]. Pleistocene and early Holocene marine transgressions occurred at river mouths and the lowest coastal areas [70]. The uppermost part of the geological sequence is represented primarily by late Pleistocene Yedoma IC and Holocene Alas complex deposits [71,72].

Quaternary deposits were studied previously on key outcrops in the Kolyma, Alazeya, Bolshaya Chukochya and Kon’kovaya rivers [71,72,73,74,75,76,77], and by drilling boreholes [44,70,78]. The Ice Complex is a specific polygenetic geological formation and represented by syngenetically frozen sediments, composed mainly of silt, and penetrated by large polygonal ice wedges [79,80,81,82,83]. Yedoma IC is the youngest Ice Complex with a thickness up to 40–50 m, which has formed during the late Pleistocene from 60 ka to 12.5 ka BP [42,52,84]. During the late Pleistocene, extensive accumulative plains of Yedoma IC dominated the landscape [71,72]. At the end of the late Pleistocene–early Holocene, due to the climate warming, the Yedoma IC was significantly reworked by thermokarst, thermoerosion and thermodenudation processes [52]. Thawing of highly ice-rich permafrost led to thermokarst lake development and thermokarst depressions (alases) formation after lake drainage. The sediments formed within these topographical depressions are called Alas complex and consist of three layers—taberal, lacustrine, and peat, which reflect three stages of the alas development [52,85,86]. The resulting landscape of uplands and interconnected drained thaw lake basins formed the specific yedoma–alas relief (Figure 2). Nowadays, yedoma–alas landscapes are characterized by a high degree of relief dissection, where the relative elevation of the yedoma surfaces over the alas bottoms is on average 20–25 m [53,87].

The main geomorphological features are represented by Yedoma IC remnants—yedoma uplands, drained thaw lake basins (alases), river valleys, marshes, and deltas (Figure 2). Elevation gradually decreases from 70–90 m a.s.l. in the south of the Kolyma Lowland tundra zone to 5–10 m at the coast of the East Siberian Sea in the north. Radiocarbon dating revealed that significant thawing of IC deposits and formation of thermokarst lakes (a so-called “thermokarst wave”), occurred within a short period around ca. 13–12 ka BP [52,88].

As a result, the yedoma–alas relief had already been widely formed mostly by ca. 10–8.5 ka BP [52]. There have been several distinguished elevation levels characterized by neotectonic activity with different rates and direction that are largely controlled by the degree to which thermokarst and thermal erosion processes reworked late Pleistocene yedoma plains in the Holocene and subsequently limnicity (fraction of lake area) properties [53,87,89,90,91]. The hydrologic system evolved concurrently with the emergence and further expansion of thermokarst lakes, which became its integrated part [86,90].

In the region, continuous permafrost is present with depths of up to 500–600 m [92]. Permafrost temperatures vary from −7 °С on the southern tundra boundary up to −10 °С on the north-eastern part [93]. Permafrost has warmed recently in monitoring sites within the study region (cape Maly Chukochy, Bolshaya Chukochya river middle stream) at 0.1 °С/year [94]. Active layer thickness (ALT) is 20–70 cm, depending on the landscape conditions. At most sites, a trend of increasing ALT is observed [95,96].

The Kolyma Lowland tundra is the eastern part of the most continental Arctic tundra region and exhibits one of the coldest permafrost conditions. Mean annual air temperatures vary from −15 °С in the coastal region to −11 °С at the boundary with the northern taiga [97]. Precipitation decreases from coastal regions from 160–210 mm to 130–190 mm in the southern part of the tundra zone. Currently, interannual air temperature increase is observed along with an increase of precipitation [96,98].

Typical soil types are cryosols, formed on well-drained yedoma uplands, gleysols, and histosols within bogged alas depressions. The region is located within the zones of typical and southern tundra [99], which corresponds to the C, D, E vegetation subzones of the Circum Arctic Vegetation Map [100]. For a typical tundra on drained yedoma uplands, a hummocky graminoid, carex, dwarf willow, dryada, lichen, moss vegetation community is common.

3. Materials and Methods

3.1. Geomorphological Analysis

3.1.1. Yedoma–Alas Relief Types

We analyzed the yedoma upland distribution based on an improved Quaternary deposits map by Veremeeva and Glushkova using Landsat imagery [53], which was created due to the lack of the geological maps at 1:200,000 scale, a digital elevation model (DEM) and the existing geological map of 1:1,000,000 scale. The mapping was done on a work scale of 1:30,000 using ArcGIS (ESRI) and significantly improved the knowledge about the Yedoma IC distribution and its degradation by thermokarst processes. This map (Figure 3) shows that the Yedoma IC covers an area of 6715 km² or 17.3% from all yedoma–alas area (15.2% of the study region), which is 2.5 time less than on geological map of 1:1,000,000 scale [70]. The alases cover 32,062 km² or 82.7% of all yedoma–alas area (72.6% of the study region). Alluvial deposits cover the 4027.5 km² or 9.1% from all area and others—alluvial–marine, marine, and bedrock—1342.5 km² or 3% of study area.

The main landforms are yedoma uplands, alases, thermokarst lakes and river valleys (Figure 2 and Figure 3). For further geomorphological analysis of lake dynamics, we only used the yedoma and alas areas, and excluded the large and medium river valleys, deltas, marshes, and bedrocks. Based on mapped yedoma and alas areas, we created a yedoma upland density map related to yedoma–alas area. The yedoma fraction in % was quantified in percentages using a 1 × 1 km-spaced grid: high (32–90%), moderate (15–32%), and low yedoma fraction (0–15%, Figure 3). The intervals were defined on a natural breaks approach preferably used for linear distribution type data, which is characteristic for the yedoma area fraction frequency distribution.

3.1.2. Elevation Analysis Using TanDEM-X DEM Data

We used the high quality TanDEM-X DEM with 12 m spatial resolution for the analysis of the relief height to assess how relief may impact lake area dynamics. The TanDEM-X DEM was derived from Synthetic Aperture Radar (SAR) interferometry based on TerraSAR-X and TanDEM-X radar twin satellites operated by the German Aerospace Center (DLR) since 2010. The vertical accuracy of the DEM with a Root Mean Square Error (RMSE) of ±1.1 m [101] for areas with low or sparse vegetation is sufficient for our mesoscale research at 1:100,000–1:200,000. We quantified the average elevation (m a.s.l.) for each yedoma–alas relief type and the main landforms: yedoma uplands, alases, and thermokarst lakes. The elevation of each lake was extracted from the DEM using the extraction tool in ArcGIS. Elevation variability of lake surfaces was no more than 1 m and, therefore, within the range of the TanDEM-X DEM accuracy for low vegetation areas.

Alas and lake heights were mostly located within the 0–50 m a.s.l. range and accordingly we used relief levels of <10, 10–20, 20–30, 30–40, 40–50 and >50 m a.s.l. for classifying the lake elevation distribution and the subsequent analysis of lake dynamics for the 1999–2018 period.

3.1.3. Thermokarst Lake Distribution

To first estimate the spatial distribution of lakes on an object level, total thermokarst lake coverage was quantified based on Landsat 8 images, acquired during 2013–2014. To ensure consistency across both years regarding surface moisture, only images acquired from August until September have been used. Atmospheric correction of each image was done for radiometric normalization across the dataset. An increase in ground resolution of the 30 m multi-spectral data was achieved through resolution merge (pan-sharpening) with the panchromatic channel to 15 m pixel size. Subsequent mosaicking, classification, and raster to vector conversion was done for the entire Kolyma Lowland tundra. Analysis of the lake area distribution based on relief type map and TanDEM-X DEM data was done using the lakes derived from 2013–2014 Landsat-based dataset.

3.2. Defining Time Frames Based on Hydrometeorological Data

As the Landsat-based trend analysis shows the average values of the lake area changes for the 1999–2018 period, we preliminarily analyzed the hydrometeorological data to assess whether there were periods with significant differences in the weather regime. For data analysis, the Chersky weather station was chosen, as it is the closest station to the study region with the most comprehensive observation records (Figure 1).

The 1999–2018 period shows that after a steady increase, increasing fluctuations of positive air temperature derived thawing degree days (TDD) became common, where minimum and maximum values alternate between 1000 and 2000 during the thawing season. Averaged across all years, an increase of precipitation could be observed, while since 2014 average precipitation of a hydrological year remained at a 1.5 times higher level than in 1999–2013 (350 mm and 220 mm, respectively) (Figure 4). Precipitation either as rain or snow during the warm and cold seasons as well as only for summer months also had constantly higher values during 2014–2018. Therefore, the two periods 1999–2013 and 2014–2018 could be differentiated as periods with moderate versus increased precipitation.

3.3. Thermokarst Lake Area Changes for the 1999–2013 and 1999–2018 Periods

Thermokarst lake area changes were determined for two periods 1999–2013 and 1999–2018. With this differentiation we aim to understand differences in lake area changes based on climatic conditions, particularly precipitation, and the influence of geomorphological factors.

3.3.1. Data and Change Detection Analysis

Thermokarst lake area changes were studied based on the Landsat data archive for the 1999–2013 and 1999–2018 periods. The processing follows the lake change processing workflow from Nitze et al. [22,23]. We used Landsat Collection-1 data from TM, ETM+ and OLI sensors with a spatial resolution of 30 m and six overlapping spectral bands for the peak summer months July and August, and a cloud cover of less than 70%. Images were preprocessed to surface reflectance by the ESPA processing interface of the United States Geological Survey (https://espa.cr.usgs.gov). In total, 1303 scenes were processed with an automated workflow, undergoing several necessary processing steps, such as masking, data distribution and calculation of multi-spectral indices (MSI: Landsat Tasseled Cap—Brightness, Greenness, Wetness, NDVI, NDWI, NDMI)). MSI were calculated for each unobstructed (cloud-, shadow- and snow-free) observation between 1999 and 2018. A robust linear trend analysis (Theil-Sen) was applied to each pixel for the temporal representation of changes of different land surface properties over both observation periods.

3.3.2. Thermokarst Lake Trend Analysis and Area Change Calculation

For calculating lake area changes, we used the machine-learning classification (MLC) of the pre-calculated spatio-temporal trend information and object-based image analysis (OBIA). Supervised MLC of robust trends of the six MSI based on Random Forest method was used to classify land-water interactions with four classes; stable classes “land” and “water” and dynamic classes “land to water” and “water to land”. Training areas were distributed across different Arctic and Subarctic regions (see Nitze et al. [23]). Based on the classification results, the lake objects were detected and differentiated into stable and dynamic zones. Stable zones were assigned a class weight of 1, where the entire pixel was covered by water. In the dynamic margin zones, subpixel probabilities of the MLC were used as weights for fractional pixel area of each class. Lakes with a maximum extent smaller than 1 ha were automatically removed from the analysis. A more detailed description of the lake area dynamics analysis is provided in Nitze et al. [22,23]. The Landsat-based trend approach is designed to detect long-term lake dynamics to minimize the influence of short-term fluctuations (e.g., due to the high precipitation events, or seasonality), and requires periods longer than 10 years to safely calculate meaningful results. Therefore, it was not possible to calculate lake area dynamics separately for the periods 1999–2013 and 2014–2018. However, period intercomparisons allow the estimation of commonalities and differences including the impact of precipitation changes.

3.3.3. Detection of Lake Expansion Along Shores Composed of Yedoma IC

Individual increases in lake area are caused by thermal erosion of lake shores. It is important to quantify the Yedoma IC lake shore thermoerosion, as the eroded ice- and organic-rich Yedoma IC deposits deliver poorly decomposed OM into lakes which accelerates the greenhouse gas (GHG) emissions and especially methane emission from lakes [46,102]. We used the Quaternary deposits map to detect those lakes with shores formed by Yedoma IC deposits. Yedoma lakeshores are often characterized by active baydzherakh slopes. Baydzherakhs are residual thermokarst mounds with heights of up to 4–6 m and that formed due to the polygonal ice wedges thawing of Yedoma IC (Figure 5a).

Baydzherakhs are typical for the yedoma upland slopes and could be inactive covered by vegetation and exhibit developed soil profiles, or active with vegetation-free mounds and pioneer plant species in depressions between mounds (Figure 5a). Active baydzherakh slopes are well detectable on very high-resolution satellite images (Figure 5b) as well as on medium resolution Sentinel 2 data (Figure 5d). We created a Sentinel 2 mosaic using July–September images of 2018 to detect active baydzherakh slopes formed by Yedoma IC for the study region (Figure 5c). Mapping expanding lakes with thermal erosion of Yedoma IC at the lake shores was done manually in ArcGIS (ESRI) using the Sentinel 2 mosaic, the yedoma distribution map, and the Landsat-based trend data of the lake area dynamics.

3.4. Database of Thermokarst Lake Changes

The Quaternary deposits map was used to distinguish lakes with primarily thermokarst origin and located within the yedoma–alas domain from other lakes in floodplains, marshes, etc. with heterogeneous or non-thermokarst origin. For instance, river valleys cover about 10% of the research area and floodplain lakes exhibit significant inter- and intra-annual area changes. For our focus on thermokarst lake dynamics we therefore excluded such lakes of oxbow–thermokarst and other origin from the further analysis. Based on the yedoma–alas relief types map, we derived lakes for each relief type using overlay operations. Trend analysis data of lake area changes present each lake as one object with information about changes. After manual inspection and correction of relief type boundaries, intersecting lakes were assigned to a specific relief type.

Based on the Landsat trend analysis and change calculation, each lake object has associated data of lake area at start and end, net change, and decrease and increase values for 1999–2013 and 1999–2018 periods which were used for the further analysis. Thermokarst lakes were divided by area size into four groups: small (1–10 ha), medium (10–100 ha), large (100–1000 ha) and very large (>1000 ha). The size groups correspond to other classification efforts (e.g., Plug et al. [17]). Using these size classification categories will allow the inter-regional comparison with other studies. Finally, each lake object carries its geospatial information, change data, relief types, absolute heights, and elevation level.

4. Results

4.1. Yedoma–Alas Relief Type Characteristic

Areas with a low yedoma fraction have the highest area of alases covering 96.4% of the class area with a lake area of 18.5% from alas area (

Table 1. Yedoma, alas, and lake area for yedoma–alas relief types.

Yedoma Fraction	Yedoma–Alas Area, km2/% from All Yedoma–Alas Area, km2	Yedoma Area, km2/% from Class Area	Alas Area, km2/% from Class Area	Thermokarst Lake Area (2013–2014)
				km2	% from Class Area	% from Alas Area
Low	17,242/45	627/3.6	16,614/96.4	3283	18.5	19.8
Moderate	13,736/35	2804/20.4	10,932/79.6	1468	10.7	13.4
High	7798/20	3283/42.1	4515/57.9	5609	11	19
All yedoma–alas areas	38,776/100	6714/17.3	32,588/83	4221	14.5	17.2

). Moderate yedoma fraction areas are characterized by 79.6% of alas area, where lakes cover 9.5%. In areas with a high yedoma fraction, yedoma uplands cover 42%, while alases occupy 58% with a limnicity of 10% in the alas areas (Table 1).

In general, thermokarst lakes cover 14.5% of the entire yedoma–alas area and 17.2% of the alas areas. Most lakes are located within alas depressions and generally represent residual lakes after the initial primary thermokarst lake drained partially or fully. Only very few lakes can be considered to be primary thermokarst lakes surrounded by Yedoma IC deposits. Low yedoma fraction areas are characterized by the highest limnicity (19% of the area), while moderate and high yedoma fraction areas have a lower limnicity (about 11% of the area).

The ratio between total alas area and modern thermokarst lake area within alases reflects general limnicity changes during the Holocene until present. The decrease of thermokarst lake area is similar among sites with different alas fractions and is about 80–86% (Table 1). The slightly lower lake area in regions with moderate alas fraction likely corresponds to better drainage conditions due to the intermediate elevation position of this relief type (Figure 6 and Figure 7).

The yedoma–alas relief of the study region has an average elevation of 25.4 m a.s.l. The average height of yedoma uplands is 38.5, of alases 22.6 and of thermokarst lakes 20.7 m a.s.l. However, there is a difference between the yedoma–alas relief type height levels (Figure 6 and Figure 7). The low yedoma fraction areas are typically found within low-lying relief positions with average yedoma, alases and lakes heights of 30, 18, and 15 m a.s.l., respectively. The moderate yedoma fraction areas are characterized by yedoma average heights of 35 and 25 m a.s.l. for alases and thermokarst lakes. For high yedoma fraction areas the average yedoma heights are 43 m, alases and lakes heights are 34 m a.s.l. (Figure 6).

Therefore, yedoma–alas relief types correspond to different elevation levels but have similar relative elevation differences between average yedoma and alas heights across all relief types, which are about 10 m (Figure 6). The average heights of alases and thermokarst lake surfaces are the same for moderate and high yedoma fraction areas, but in low yedoma fraction areas lake level heights are on average 2 m lower than average alas height.

4.2. Geomorphological Analysis of Thermokarst Lake Distribution

Analysis of the lake size area distribution in relation to elevation shows that within each yedoma–alas relief type, average lake level heights of different lake size classes are grouped roughly within the same elevation range with a slightly larger difference of 10 m for low yedoma fraction areas compared to 5 m for moderate and high yedoma fraction areas (Figure 8). Average heights of small lakes are on the same level with average heights of large and very large lakes for low and moderate yedoma fraction areas. This indicates that small lakes are mostly residual lakes that remain after drainage of initial large thermokarst lakes. The interquartile range of the upper and lower elevation value distribution is the narrowest within the low yedoma fraction (from 8 to 23 m). For moderate yedoma fraction relief the elevation distribution ranges from 12 to 40 m, and for high yedoma fraction areas between 20–45 m. Generally, large and very large lakes are characterized by less variation and are mostly located in low-lying relief parts.

The largest lake area has been revealed for the low yedoma fraction regions with a significant contribution of large and very large lakes (Figure A1). In areas with a medium and high yedoma fraction, the largest contribution derives from large and medium lakes. Small lakes have a similar share in all relief types.

4.3. Lake Area Changes of 1999–2013 and 1999–2018 and Their Dependence on Geomorphological Conditions

4.3.1. Lake Area Changes within Yedoma–Alas Relief Types

A general increase in thermokarst lake area was revealed for both periods 1999–2013 and 1999–2018 (Figure 9,

Table 2. Lake area change for the periods 1999–2013 and 1999–2018.

Yedoma Fraction	Limnicity Changes in Relation to 1999		Increased Lake Area in Relation to 1999		Decreased Lake Area in Relation to 1999
	km2	%	km2	%	km2	%
1999–2013
Low	19.75	+0.62	119.72	+3.73	99.97	−3.12
Moderate	19.51	+1.29	55.68	+3.69	36.18	−2.4
High	9.95	+1.24	27.7	+3.45	17.76	−2.21
All yedoma–alas areas	49.2	+ 0.89	203.11	+3.68	153.91	−2.79
1999–2018
Low	155.97	+4.98	190.14	+6.07	34.18	−1.09
Moderate	45.87	+3.12	67.27	+4.58	21.4	−1.46
High	21.3	+2.76	30.04	+3.9	8.74	−1.13
All yedoma–alas areas	223.13	+4.15	287.45	+5.35	64.32	−1.2

). However, there is a difference in spatial distribution of lake area dynamics. In 1999–2013, lake shrinkage was observed in all areas, while in the northern part expanding lakes prevailed. When considering the full period 1999–2018, increasing lake areas prevailed in all areas with most significant increases in the region with predominantly large lakes on the north-eastern coastal part (Figure 9).

The most widespread hydrologic processes, detected by the tasseled cap trend index of the Landsat imagery for the period 1999–2013 were drainage events of thermokarst lakes, increasing surface area of the largest lakes and inundation along low-lying streams and rivers (Figure 9 and Figure A2). In contrast, the full period 1999–2018 was dominated by river inundation and thermokarst lake expansion, while only some thermokarst lakes experienced partial or almost complete drainage.

The overall limnicity of the yedoma–alas relief excluding alluvial, deltaic, and marine deposits areas, had increased slightly by 0.89% during the early period 1999–2013 and strongly by 4.15% when looking at the full period 1999–2018 (Table 2, Figure 9). During 1999–2013, the minimum of increased and maximum of decreased lake area was found for the low yedoma fraction areas (+0.62% and −3.12% respectively). In the moderate and high yedoma fraction areas stronger lake area increase (+1.29 and +1.24 respectively) and less pronounced decrease (−2.4% and −2.1% respectively) were observed (Table 2, Figure 9).

In contrast, during the 1999–2018 period, maximum lake area increase was detected for the low yedoma fraction region (+4.98%) and minimum for high yedoma fraction areas (+2.76%). Lake area decrease was most pronounced in moderate yedoma fraction areas (−1.46%). Comparing only increasing lake areas, in 1999–2013 the values for all relief types were similar, while in 1999–2018 increase values consistently rise from high- to low yedoma fraction areas (Table 2).

The highest contribution in the lake area change trend for both periods is due to lake area dynamics in low yedoma fraction areas with highest lake area percentage changes and largest water body area additions from large and very large lakes (Figure 10). Lake area change differences between 1999–2013 and 1999–2018 periods were most expressed for low yedoma fraction areas, where a substantial lake area increase was outpacing generally lower decrease events in the 1999–2018 period in comparison with the much smaller difference between increase and decrease values during 1999–2013 (Figure 10). The same pattern for moderate and high yedoma fraction areas was observed but with smaller absolute area values.

Analysis of the lake dynamics dependent on lake size shows that larger water area increase for all lake size groups in different relief types occurred in 1999–2018 in comparison with 1999–2013 (Figure A3). The strongest lake area increase during 1999–2018 was due to expanding large and very large lakes. In contrast, during 1999–2013 an increasing trend prevailed for small and medium lakes within all relief types, while large and very large lakes tended towards a decreasing trend. If considering lake area gross increase and gross decrease separately, in 1999–2013 the decreasing trend for all lake size groups of all relief types was more expressed.

To estimate the contribution of large and very large lake area dynamics we analyzed separately lakes with a significant area decrease of more than 10% and an increase of more than 5% (Figure 11 and Figure A4;

Table A1. Large and very large lakes with largest area loss in 1999–2013 and 1999–2018.

Lake Area Decrease Group	Relief Types, Yedoma Fraction	Area Changes, ha (Decreased Lake Number)
		1999–2013	1999–2018
Large lakes (area decrease 10–100%)	All Low Moderate High	−4844.3 (34) −2603.3 (16) −1548.4 (11) −692.7 (7)	−2153.8 (14) −1013.8 (6) −988.8 (5) −151.2 (3)
Large lakes (area decrease 50–100%)	All Low Moderate High	−3972.4 (15) −2090.1 (7) −1345.3 (5) −537.0 (3)	−1810.4 (7) −792.7 (3) −917.2 (3) −100.6 (1)
Large lakes (area decrease 25–50%)	All Low Moderate High	−356.4 (5) −237.0 (2) −119.4 (3) -	−133.9 −96.2 - 37.6 (1)
Large lakes (area decrease 10–25%)	All Low Moderate High	−515.5 (14) −276.1 (7) −83.7 (3) −155.6 (4)	−209.5 (5) −124.8 (2) −71.6 (2) −13.0 (1)
Very large lakes (area decrease 10–25%)	All Low Moderate High	−2124.0 (3) −2124.0 (3) - -	−110.3 (1) - −110.3 (1) -

and

Table A2. Large and very large lakes with largest area gain in 1999–2013 and 1999–2018.

Lake Area Increase Group	Relief Types, Yedoma Fraction	Area Changes, ha (Decreased Lake Number)
		1999–2013	1999–2018
Large lakes (area increase > 5%)	All Low Moderate High	1120.4 (74) 790.8 (50) 287.8 (20) 41.8 (4)	4466.8 (136) 3332.4 (98) 1050.6 (33) 83.8 (5)
Large lakes (area increase >20 %)	All Low Moderate High	35.6 (1) - 35.6 (1) -	2166.4 (6) 1576.7 (4) 589.7 (2) -
Large lakes (area increase 10–20%)	All Low Moderate High	200.5 (9) 122.8 (6) 77.6 (3) -	368.9 (14) 355.4 (13) 13.4 (1) -
Large lakes (area increase 5–10%)	All Low Moderate High	884.4 (64) 667.9 (44) 174.6 (16) 41.8 (4)	1931.6 (116) 1400.2 (81) 447.5 (30) 83.8 (5)
Very large lakes (area increase >5%)	All Low Moderate High	1087.9 (3) 1087.9 (3) - -	5749.3 (7) 5749.3 (7) - -
Very large lakes (area increase >10–20%)	All Low Moderate High	812.5 (1) 812.5 (1)	4361.9 (3) 4361.9 (3)
Very large lakes (area increase >5–10%)	All Low Moderate High	275.5 (2) 275.5 (2)	1387.4 (4) 1387.4 (4)

). In comparison with overall changes, substantial contributions of decreasing large and very large lakes in 1999–2013 was found for low and moderate yedoma fraction areas. During 1999–2018, despite relatively high values of increased lake area in comparison with overall changes, their contribution was not significant. This suggests that overall decreased large and very large lakes contributed most to overall lake area changes during 1999–2013, while the increasing trend during the extended observation period 1999–2018 was due to area gains from many lakes of different size, but not due to the increase of several very large lakes.

The distribution, area, and numbers of large and very large lakes that experienced decreasing areas by more than 10% and an increase of more than 5% separated by relief types are shown in Table A1 and Table A2. Decreased lakes number and area was about two times more in 1999–2013 in comparison with 1999–2018. The largest number and area of growing lakes were observed during 1999–2018. Only three large lakes drained substantially with 10–25% of lake area loss (Figure 11).

4.3.2. Lake Area Changes for 1999–2013 and 1999–2018 Analysis with Respect to Lake Size and Elevation Levels within the Yedoma–Alas Relief Types

Analysis of the lake area dynamics depending on absolute heights shows a difference in distribution patterns of lake changes for different relief types that in turn reflect different elevation levels of alases (Figure 12). However, a general increase trend for mostly all lake size groups at different elevation levels was observed within all relief types during 1999–2018 in comparison with 1999–2013.

In low yedoma fraction areas most lakes are located at elevation levels from <10 to 20 m a.s.l. with an average elevation of alases of 18 m a.s.l. Due to the presence of very large lakes in this relief type compared to medium and high yedoma fraction areas, the difference in area changes is most expressed here. The decreasing trend in 1999–2013 prevailed for large and very large lakes at the 10–20 m a.s.l. elevation level and very large lakes at <10 m. When considering only increasing lakes, small- and medium-sized lakes had similar values for both periods. A significant increase in 1999–2018 in comparison with 1999–2013 was detected for large and very large lakes at the lowest elevations and large lakes at the 10–20 m a.s.l. level, while other lake size groups at different elevations showed similar values.

In moderate yedoma fraction areas, the most significant increases in 1999–2018 in comparison with 1999–2013 were observed for medium and large lakes at <10 m elevation level and large lakes at 10–20 m heights. In 1999–2013, a decrease trend was observed for large lakes at <10 m, 10–20 m and 40–50 m levels, and in 1999–2018 a decrease trend prevailed only for large lakes at 40–50 m.

However, medium-sized lakes at 40–50 m and >50 m a.s.l. experienced lake area increase during the 1999–2018 period. Since the elevation levels higher than 40 m corresponds to the average heights of yedoma, these lakes have the highest potential for thermokarst development. However, the picture of decreasing large and increasing medium-sized lakes suggests that once lakes reach a certain size on top of yedoma uplands, their growth is limited by topographical constraints.

High yedoma fraction areas are characterized by a high abundance of medium lakes at different elevation levels for both periods (Figure 12 and Figure A1), with more pronounced lake area increase in 1999–2018 in comparison with 1999–2013. The maximum lake area increase trend in 1999–2018 was found at 30–40 m a.s.l. Largest decrease values were detected for large lakes in 1999–2013 at 30–40 m and 40–50 m and to a lesser degree for large lakes in 1999–2018 at 30–40 m elevations. The average yedoma elevation of high yedoma fraction areas is 43 m, which means that the lakes of 40–50 and > 50 m are at the same level of yedoma uplands and therefore have the maximum potential for the thermokarst development. On this elevation level the increase trend for both periods and all lake size groups was observed with the exception of a prevailing decrease trend for large lakes on the 40–50 m elevation level. In 1999–2018 the values for increasing lakes were higher than in 1999–2013 with the largest difference for medium lakes on elevation level >50 m.

Analysis of the lake change values distribution dependent on elevation levels shows that the highest variations are found for the low yedoma fraction area and the smallest for the high yedoma fraction areas (Figure 13). In general, in the 1999–2013 period variations were smaller than in 1999–2018, while variation increased from small to large lakes and from high to low elevations. The variations in the increasing trend were significantly larger than for the decreasing trend for all relief types. In high yedoma fraction areas the value variations on elevations higher than 40 m corresponded to the level of yedoma average heights were highest for large lakes in both periods. The smallest lakes have similar lake change variations for both periods. Medium and large lakes have higher value variations in 1999–2018 in comparison to 1999–2013 with the most pronounced changes for large lakes on the 40–50 elevation level in 1999–2018 compared to 1999–2013. Variations in lake area decrease were higher in 1999–2013 for all elevation levels and highest for 40–50 elevation level.

4.3.3. Detection of Thermal Erosion of Lake Shores Formed by Yedoma IC

We identified 342 expanding lakes with active baydzherakh slopes and expressed thermoerosion due to Yedoma IC lakeshore degradation (Figure 14 and Figure A5). They were most commonly located within the high yedoma fraction areas (Figure 14 and Figure A5), where lake area mostly increased in both periods except for some lakes that exhibited a decreasing trend in 1999–2013 which was compensated by an increasing trend during the extended observation period 1999–2018. Thermodenudation processes on active baydzherakh slopes surrounding lake depressions were detected for both periods, but were more expressed in 1999–2013. For some lakes, the Sentinel 2 images exhibit the contribution of eroded Yedoma IC sediments into lake water (Figure 14, the example 2; Figure A5, examples 3,7).

5. Discussion

5.1. Yedoma–Alas Relief Types Formation

We revealed a similar elevation difference about 10 m between average yedoma and alases surface heights within each of the different relief types (Figure 6), which could be explained by the cover character of the initial yedoma surface and suggests a similar thickness of Yedoma IC across the study region. From our point of view, this supports an eolian-dominated sedimentation and genesis of the Yedoma IC deposits on the Kolyma Lowland region. Although in some recent studies the eolian genesis of Yedoma IC is supported by the metagenom, palynological, and other analyses [103,104,105], there is also evidence of polygenetic Yedoma IC origins in other regions including alluvial, fluvial, and niveo-aeolian processes [106]. Despite the ongoing Yedoma IC genesis discussion, the concept of an initially extensive Yedoma IC cover deposits blanketing the landscape is crucial for understanding the spatial development of thermokarst during the Holocene and resulting yedoma–alas relief types. It also determines the modern and future changes of these landscapes.

Our results also support preliminary estimates of average incision depth of 25 m by thermokarst by Veremeeva and Glushkova [53], in line with a thickness of about 20–25 m detected by drilling and outcrop studies in the region [71,74,76,78] and proving the validity of a geomorphometric approach. However, the elevation difference of about 10 m between yedoma upland average heights and alas floors measured using a TanDEM-X DEM in comparison to 25 m from the relief dissection map is explained by using the minimum and maximum heights in the relief dissection map and averaged elevations across yedoma and alas areas with the DEM. This mainly reflects the remnant character of yedoma uplands with extensive gentle slopes around their margins that were shaped by thermodenudation and surface subsidence [53]. The drainage conditions are more favorable on yedoma uplands and slopes due to the higher ALT [107]. However, Günther et al. [27] found stronger land surface subsidence rates where ALT is low, and these were mostly located on the most elevated parts of a yedoma upland. This pattern has been repeatedly confirmed by Chen et al. [108] and can be attributed to pure ground ice close to the surface on yedoma tops in some cases.

An important question and largely understudied aspect is the thickness of buried Yedoma IC under alas deposits as well as under the lake taliks. Drilling data [44,109] within the Kolyma Lowland showed a 2–10 m thickness, which means that Yedoma IC thawing in some cases has not been complete during the Holocene thermokarst wave and suggests that lake drainage events play an important regulatory role in the landscape development. Therefore, thermokarst lake area loss within classes of different alas fraction can be considered to be the result of a simultaneous development of thermokarst lakes and the fluvial drainage network during the Holocene [90,110]. Thermokarst lakes are closed systems only during their initial formation stage while actively thawing Yedoma IC. Subsequent arising runoff from lakes either quickly creates new water flow pathways [111] or connects them with existing drainage system and ultimately with the absolute erosion basis, establishing a dynamic equilibrium of the hydrological network [110]. Similar limnicity changes within the different settings of yedoma–alas relief types confirm that the hydrological network connections down to the absolute erosion basis level are an important factor for lake area dynamics and drainage. This conclusion is supported by a study of the Arctic Coastal Plain in Alaska where the lake coverage is fairly constant (~20%) while thermokarst basin coverage varied from 19 to 47% depending on geological and geomorphological conditions [56]. Our limnicity quantification of 14% revealed that the Kolyma Lowland tundra is in line with other Yedoma IC thermokarst regions. Thermokarst lakes cover about 14% of the Bykovsky Peninsula [112], 9–15% of the area in the northern part of the Anabar–Olenek interfluve [113,114], 13.3% of the third terrace in the Lena river delta [115], 15% of the Shirokostan Peninsula [40], and 7% of the Seward Peninsula, Alaska [18].

5.2. Comparison with Previous Data from Yedoma Regions

Recent studies report about the general trend of decreasing thermokarst lake area in different Arctic regions of the continuous permafrost zone [23,30,31,33]. On the Seward Peninsula in Alaska, the lake area decreased by 14.9% from 1950 to 2007 [18]. In the Yana-Indigirka Lowland tundra zone and the Kolyma Lowland taiga zone over the 1970–2000 period, lakes increased by 0.9%, while most lakes decreased at 2.2% [21]. On Kurungnakh Island, a yedoma remnant in the southern Lena river delta, over the period 1964–2006 a general trend of lake area decrease (3.5%) has been revealed, while persistent lakes increased about 2% [116]. In contrast, some works report about lake area increasing trends. Walter et al. [117] revealed a lake area increase of 14.4% for a period 1974–2000 in a small area on the right bank of the Kolyma river near Chersky. A prominent example of lake area increase is Central Yakutia, located within the high continental middle taiga zone [22,36,37,38]. For the 1999–2014 period a lake area increase of 48.8% was calculated by Nitze et al. [22] in this particular region. Over an extended historical observation period from 1944, the lake area increase occurs in alases as well as formation of new lakes on yedoma uplands, which is related to an increase of air temperature, winter precipitation and ALT [38].

In a previous study of the eastern part of Kolyma Lowland, Nitze et al. [22] detected a general trend of limnicity decrease of 0.51% over the 1999–2014 period. Although this is not entirely consistent with the increasing trend of 0.89% for the period 1999–2013 in our study, the difference derives from the focus on lakes with thermokarst genesis only in this study. On top of it, the decreasing trend in the previous study has been strongly influenced by lake area changes within the Khallerchin tundra, which is characterized by ice-poor sandy deposits and that has not been a subject in our research.

Lake dynamics in the Yedoma region of northwestern Alaska revealed catastrophic thermokarst lake drainage events in 2018 [33]. There, 192 lakes were completely or partially drained following the winter in 2017/2018, which was the warmest and wettest winter on record and that exceeded the average drainage rate by a factor of ~10 and doubled the rates of the previous extreme lake drainage years of 2005 and 2006.

Based on these reports from different regions in the Arctic except for Central Yakutia, the most recent increase trend of thermokarst lakes in the Kolyma Lowland was a quite unexpected result and highlights the importance of differentiating between lake genesis types. The combination of different factors such as geographical settings, geological and geomorphological conditions, and regional climate peculiarities led to an increasing trend in Kolyma Lowland tundra.

5.3. Factors of Modern Thermokarst Lake Dynamics and Potential for Future Thermokarst Development

5.3.1. Climate Factors

NE Siberia is one of the Arctic regions where the strongest climate warming is observed since 1960 during both cold and warm seasons. In the region, the air temperature increase for the period 1961–2014 was 2–4 °C for both seasons [118] and a wetting of climate was also observed. The increase of air temperature and precipitation is likely connected to the reduced sea ice cover in the Arctic Ocean and the East Siberian Sea [119].

Annual air temperature and precipitation at the Chersky weather station show an increasing trend for both parameters since 1940 (Figure 15). Daily air temperature analysis for June-September period of the Chersky weather station since 2005 shows no increasing trend, while daily precipitation data reveals that since 2014 precipitation did increase (Figure A6). It is important to note that more heavy rain events were observed in August and July during the 2014–2017 period, while in 2018 the days with maximum precipitation were observed earlier at the end of June, when a heavy rainfall with 50 mm fell for three days, which is twice the amount of the average precipitation sum for the summer months (Figure A6). Because at the end of June the ALT is still quite shallow, heavy rainfall and poor drainage conditions likely played a role in rising water level in lakes and rivers in 2018.

Monitoring of ALT changes since 1996 for Kolyma Lowland tundra region CALM (Circumpolar Active Layer Monitoring program network) sites reveal the general trend of increasing ALT on the yedoma upland sites while there is no or negative trend of the alas sites ALT (Figure 16).

However, as alases dominate the area, but ALT in alases did not increase, while at the same time ALT on yedoma uplands increased, this together with stronger precipitation could have supported the lake area increasing trend due to the poor drainage conditions in low yedoma fraction areas and Yedoma IC thawing in areas with moderate and high yedoma fraction.

The increasing trend of lake area for both periods is likely the result of a warming climate and wetting, which also led to permafrost thawing. Therefore, the lake water levels could rise due to the deepening of the ALT on yedoma uplands together with stable ALT on alas depressions, activation of thermokarst and thermoerosion, as well as due to the increased precipitation as snow in cold and rain in warm time.

5.3.2. Geomorphologic Factors of the Lake Area Changes

Yedoma–alas relief types reflect the intensity of the thermokarst processes in the past, dependent on elevations of the initial surface, and the geological conditions, formed the modern relief and spatial patterns of lake distribution as well as geomorphological conditions for thermokarst lake dynamics. Despite the generally decreasing lake area trend and drainage of 80% of existing thermokarst lakes during the Holocene, our study shows an increasing trend, especially expressed in the most recent years. However, the spatial analysis of lake area changes during both periods, with moderate and high precipitation, revealed a difference in dynamics dependent on relief type, elevation, and lake size. In low yedoma fraction areas, the strongest dynamics of large and very large lakes were detected at the lowest elevation levels: decreasing in 1999–2013 and increasing lake area for both periods, but in 1999–2018 the increase was about three times higher. The areas with moderate yedoma fraction play a transitional role with largest lake area changes at the 10–20 m elevation level and significant changes for elevations of 20–50 m. High yedoma fraction areas have the lowest absolute changes in area. Highest lake area changes were detected for medium lakes on 40–50 m. We revealed an increasing lake area trend for both periods, but the increase was most pronounced in 1999–2018. Important to note is that increase trends during both periods were observed despite different weather conditions with moderate precipitation during 1999–2013 and increasingly higher precipitation starting 2014 during the 1999–2018 period. As mentioned above, the Yedoma IC could be present under lake bottoms. In high yedoma fraction regions this could be an important factor for lake deepening. In the region of Cape Maly Chukochy, where lakes are surrounded by Yedoma IC, we detected baydzherakhs on the lake bottom using remote sensing data, indicating ongoing downward thermokarst activity underneath lakes (Figure A7). Our results indicate the continued Yedoma IC degradation due to thermokarst and thermoerosion in areas with high yedoma fraction.

5.4. Potential for GHG Emissions Assessment

The observed thermokarst lake area increase likely acts to accelerate GHG emissions, especially methane, which was shown for other Arctic regions [9,10,45,120]. For the Yedoma region of the Yana-Indigirka Lowland, Dean et al. [48] show that >80% of total inland water carbon was contemporary in age, but pre-age Holocene and late Pleistocene carbon contributed >50% at sites strongly affected by permafrost thaw. Brosius et al. [121] conclude that radiocarbon age profiles from thermokarst lake in Yedoma region (Alaska) showed oldest methane-carbon origins at lake margins and younger profiles towards the center, indicating degradation of freshly thawed old permafrost C associated with lake expansion. In an ongoing warming and wetting climate, the lake area increasing trend likely will continue. Geomorphological analysis and classification of yedoma–alas relief types allow to consider lake area changes separately for regions with low, medium, and high yedoma fraction, which could prove useful in GHG upscaling efforts, especially methane emission from Yedoma regions.

6. Conclusions

We investigated lake area dynamics for the 1999–2018 period of the Kolyma Lowland tundra Yedoma region based on an approach that includes remotely sensed thermokarst lake changes, geomorphology and DEM analysis, and Quaternary deposit mapping. We distinguish yedoma–alas relief types with low, moderate, and high yedoma fraction area (0–15, 15–32 and 32–90% respectively), which reflects thermokarst process intensity depended on geological and geomorphological conditions during the Holocene. We found that yedoma–alas relief types have different elevation levels but similar height differences of about 10 m between average yedoma and alases heights for all relief types, suggesting the initial cover deposit character of Yedoma IC. We revealed a limnicity decrease of the 80–86% that is similar for all yedoma–alas relief types during the Holocene, highlighting the importance of the connection of thermokarst lakes with the hydrological network.

Despite the general lake decreasing trend over the Holocene, modern climate warming led to an increasing trend of thermokarst lake coverage in the region. We considered two periods 1999–2013 and 1999–2018 dependent on meteorological conditions for Landsat-based trend analysis and change calculation. During the 1999–2013 period with moderate precipitation, lakes increased in size by 0.89%. In 1999–2018, this figure increased to 4.15% of lake area expansion due to the prevailed increasing trend for the 2014–2018 period with higher precipitation.

We analyzed the influence of geomorphological factors on lake area changes for both periods depending on yedoma–alas relief type, lake size, and elevation using DEM data. Analysis of the lake area dynamics for both periods depending on their absolute height shows a difference in the distribution pattern of lake changes. However, the vital increase trend during 1999–2018 of mostly all lake size groups was observed at all elevation levels and for all relief types in contrast to the 1999–2013 period, when this pattern was more heterogeneous with shrinking and expanding lakes. Large and very large lakes substantially decreased and contributed most to the overall lake area changes in 1999–2013, while in 1999–2018 the increasing trend was due to lake area gains from a broader base of many lakes and not due to the increase of several large lakes. We find the largest increase of lake area in low yedoma fraction areas. In high yedoma fraction areas, the increasing trend was most expressed for medium-sized lakes located on yedoma uplands. This indicates further Yedoma IC degradation due to thermokarst and thermoerosion processes in areas with high yedoma fraction.

Our results contribute to an enhanced understanding of thermokarst lake dynamics and support the development of impact assessments for landscape and habitat characteristics as well as local land-use management. Geomorphological analysis and differentiation of yedoma–alas relief types allow consideration of lake area changes separately for regions with varying yedoma fractions, which could be used in refined assessments of GHG and especially methane emission in Yedoma regions. Because hydrometeorological conditions and patterns may change quickly over short periods, a continuous lake dynamics monitoring that builds on geomorphological analysis is necessary.

We showed that geological and geomorphological conditions played an important role of thermokarst development during the Holocene and are necessary to analyze spatial patterns of modern thermokarst lake dynamics.