Geomorphology and InSAR-Tracked Surface Displacements in an Ice-Rich Yedoma Landscape

Ice-ridge Yedoma terrain is susceptible to vertical surface displacements by thaw and refreeze of ground ice, and geomorphological processes of mass wasting, erosion and sedimentation. Here we explore the relation between a 3 year data set of InSAR measurements of vertical surface displacements during the thaw season, and geomorphological features in an area in the Indigirka Lowlands, Northeast Siberia. The geomorphology is presented in a geomorphological map, based on interpretation of high resolution visible spectrum satellite imagery, field surveys and available data from paleo-environmental research. The main landforms comprise overlapping drained thaw lake basins and lakes, erosion remnants of Late Pleistocene Yedoma deposits, and a floodplain of a high-sinuosity anastomosing river with ancient river terrace remnants. The spatial distribution of drained thaw lake basins and Yedoma erosion remnants in the study area and its surroundings is influenced by neotectonic movements. The 3 years of InSAR measurement include 2 years of high snowfall and extreme river flooding (2017–2018) and 1 year of modest snowfall, early spring and warm summer (2019). The magnitude of surface displacements varies among the years, and show considerable spatial variation. Distinct spatial clusters of displacement trajectories can be discerned, which relate to geomorphological processes and ground ice conditions. Strong subsidence occurred in particular in 2019. In the wet year of 2017, marked heave occurred at Yedoma plateau surfaces, likely by ice accumulation at the top of the permafrost driven by excess precipitation. The spatial variability of surface displacements is high. This is explored by statistical analysis, and is attributed to the interaction of various processes. Next to ground ice volume change, also sedimentation (peat, colluvial deposition) and shrinkage or swelling of soils with changing water content may have contributed. Tussock tundra areas covered by the extreme 2017 and 2018 spring floods show high subsidence rates and an increase of midsummer thaw depths. We hypothesize that increased flood heights along Siberian lowland rivers potentially induce deeper thaw and subsidence on floodplain margins, and also lowers the drainage thresholds of thaw lakes. Both mechanisms tend to increase floodplain area. This may increase CH4 emission from floodplains, but also may enhance carbon storage in floodplain sedimentary environments.

Ice-ridge Yedoma terrain is susceptible to vertical surface displacements by thaw and refreeze of ground ice, and geomorphological processes of mass wasting, erosion and sedimentation. Here we explore the relation between a 3 year data set of InSAR measurements of vertical surface displacements during the thaw season, and geomorphological features in an area in the Indigirka Lowlands, Northeast Siberia. The geomorphology is presented in a geomorphological map, based on interpretation of high resolution visible spectrum satellite imagery, field surveys and available data from paleo-environmental research. The main landforms comprise overlapping drained thaw lake basins and lakes, erosion remnants of Late Pleistocene Yedoma deposits, and a floodplain of a high-sinuosity anastomosing river with ancient river terrace remnants. The spatial distribution of drained thaw lake basins and Yedoma erosion remnants in the study area and its surroundings is influenced by neotectonic movements. The 3 years of InSAR measurement include 2 years of high snowfall and extreme river flooding (2017)(2018) and 1 year of modest snowfall, early spring and warm summer (2019). The magnitude of surface displacements varies among the years, and show considerable spatial variation. Distinct spatial clusters of displacement trajectories can be discerned, which relate to geomorphological processes and ground ice conditions. Strong subsidence occurred in particular in 2019. In the wet year of 2017, marked heave occurred at Yedoma plateau surfaces, likely by ice accumulation at the top of the permafrost driven by excess precipitation. The spatial variability of surface displacements is high. This is explored by statistical analysis, and is attributed to the interaction of various processes. Next to ground ice volume change, also sedimentation (peat, colluvial deposition) and shrinkage or swelling of soils with changing water content may have contributed. Tussock tundra areas covered by the extreme 2017 and 2018 spring floods show high subsidence rates and an increase of midsummer thaw depths. We hypothesize that increased flood heights along Siberian lowland rivers potentially induce deeper thaw and subsidence on floodplain margins, and also lowers the drainage thresholds of thaw lakes. Both mechanisms tend to increase floodplain area. This may increase CH 4 emission from floodplains, but also may enhance carbon storage in floodplain sedimentary environments.

INTRODUCTION
Permafrost is increasingly influenced by climate warming. Permafrost thaw is expected to enhance greenhouse gas emissions from permafrost soils and to decrease stability of soils, in particular in ice-rich permafrost areas (e.g. IPCC, 2013;AMAP, 2017). A specific type of ice-rich permafrost formation that occurs over large areas in Siberia and Alaska is Yedoma, also known as Ice Complex: periglacial loess and fluvial silt deposits that were deposited during Quaternary glacials in unglaciated terrain, which obtained large masses of excess ice by the growth of syngenetic ice wedges during sedimentation (Schirrmeister et al., 2013;Murton et al., 2015;Strauss et al., 2017). Yedoma contains a considerable carbon stock (Zimov et al., 2006;Strauss et al., 2017). These deep ice-rich deposits (several tens of meters and more) are very vulnerable to abrupt thaw (i.e. thermokarst) with rapid erosion, or more gradual subsidence, and may release carbon to the atmosphere as CO 2 and CH 4 (Schneider von Deimling et al., 2015;Vonk et al., 2013;Vonk et al., 2015). Surface displacements related to permafrost thaw also affect the stability of human infrastructure in permafrost areas. Structural instability can lead to severe environmental damage when urban, industrial or mining infrastructure is affected (e.g. Shiklomanov et al., 2017;Hjort et al., 2018).
Climate-related changes that may enhance permafrost thaw are higher air temperatures in winter and summer, and changes in precipitation that affect snow cover thickness, soil water saturation, water ponding, river flooding, fire frequency and vegetation change (e.g. reviews by Anisimov and Nelson, 1997;Nelson et al., 2002;Jorgenson et al., 2010;Nauta et al., 2014). Vegetation change also may delay permafrost thaw (Blok et al., 2010;Jorgenson et al., 2010;Kanevskiy et al., 2017). The effects of these climatic factors act unevenly across the landscape, since they are affected by geomorphologically determined terrain characteristics as slope, insolation, drainage and exposure to river flooding and erosion.
Net loss of ice by thaw of ice-rich permafrost over several years is reflected by a trend of surface subsidence, superimposed on the yearly surface displacements due to thaw at the top of the permafrost in summer (subsidence) and regrowth of ice in autumn and winter (frost heave). These surface displacements can be detected by synthetic aperture radar (SAR) interferometry from satellites (InSAR), and can assist in outlining areas that are most vulnerable to soil subsidence from permafrost thaw (Liu et al., 2010;Liu et al., 2012;Liu et al., 2014;Teshebaeva et al., 2020). Ice content is known to be an important parameter in surface subsidence (e.g. Jorgenson et al., 2010;Schneider von Deimling et al., 2015;Turetsky et al., 2019;van Huissteden, 2020), at least on a local scale. The soil ice distribution in Yedoma landscapes is also related to geomorphology and (cryo) stratigraphy of deposits, that determine duration and mode of ice accumulation: growth of wedge ice and small-scale cryostructure ice. This varies from large scale variation, e.g. large ice volumes in deep syngenetic wedges in Yedoma vs. lower ice volumes in deposits of Holocene thaw lake beds and floodplains (Kanevskiy et al., 2013;Morgenstern et al., 2013;Strauss et al., 2017), and variations determined by sedimentary patterns of floodplain and thaw lake sediments (Ulrich et al., 2014). Thus, the spatial variation of vulnerability to surface subsidence is expected to be strongly related to geomorphology.
Here, we present a geomorphological map of an Eastern Siberian Yedoma area in the Indigirka Lowlands, combined with InSAR surface displacement measures using Sentinel-1A/ B data collected over 3 years. This time series is too short to quantify long term surface displacement trends. However, the data result from years with strongly contrasting seasonal weather patterns: 2 years of high winter precipitation resulting in pronounced spring flooding, and a year with low winter precipitation, early spring and warm summer. Therefore these data contain information that permits detection of differences in surface stability between terrain units, differences among years with contrasting weather patterns, and inferences on the processes that determine vertical surface displacements, such as ground ice thaw and accumulation, river flooding, erosion/ sedimentation and soil processes. Also neotectonics influence the landscape on a large scale.
We explore the relation of surface displacements with geomorphological terrain units, and the interaction of these with the contrasting yearly weather patterns, by statistical analysis. Our analysis shows, that patterns of surface displacements vary with geomorphological position, and between years.
The carbon cycle and greenhouse gas emission impacts of geomorphological changes, changes in active layer thickness (ALT), soil subsidence and erosion as detected by InSAR are potentially large (e.g. Liu et al., 2010;Liu at al., 2014;Abbott and Jones, 2015;Natali et al., 2015). We discuss the potential carbon cycle effects of our findings in general terms.

STUDY AREA
The study area is located at the Kytalyk/Chokurdagh tundra research station (70°49′N, 147°29′E), approximately 30 km WNW of the town of Chokurdagh, and located just north of the Berelegh river, a tributary to the Indigirka river ( Figure 1). From 2003 onwards, data were collected on the carbon balance of the tundra and permafrost soil (Van der Molen et al., 2007;Blok et al., 2010;Parmentier et al., 2011a;Parmentier et al., 2011b;Parmentier et al., 2011c;Nauta et al., 2014;Dean et al., 2020), as well as data on vegetation, heat balance, soils, soil ice content and environmental history (Juszak et al., 2014;Beermann et al., 2015;Juszak et al., 2016;Teltewskoi et al., 2016;Weiss et al., 2016;Juszak et al., 2017;Li et al., 2017;Wang et al., 2018;Magnússon et al., 2020). The site is located on continuous permafrost. It is dominated by three main terrain units: the Berelegh river plain, Yedoma terrace remnants dating from the Last Glacial, and Holocene drained thaw lake basins (DTLB). The area has recently (2007,2017,2018) experienced several high to extremely high spring snowmelt floods (This paper; Tei et al., 2020). The geomorphological map is based on high resolution visible spectrum images (Geo-Eye, Worldview), high resolution digital elevation models (DEMs) derived from radar satellites, and terrain surveys. The permafrost is characterized by a shallow active layer, overlying continuous permafrost of over 300 m thickness. Near the surface, it is underlain by ice-rich Yedoma deposits (Schirrmeister et al., 2011), Holocene fluvial deposits of the Berelegh river and drained thaw lake basin fills (Alas deposits hereafter). The nearest weather station is that of Chokurdagh airport, WMO station code 21,946, 27 km away from the study site. The mean annual air temperature is −13.4°C, with an average July temperature of 10. 3°C (1981-2010). Mean annual precipitation is 196 mm, with 76 mm falling in June to August (1981August ( -2010 (Trouet and Van Oldenborgh, 2013).
The data from the weather station ( Figure 2A, NOAA, 2020) show a conspicuous warming trend, in particular in the autumn months (October, November), winter and late winter (January, March, April). The winter precipitation shows marked variability in recent years, with three consecutive years of excessive snowfall in 2016-2018 ( Figure 2B). The high snowfall winters of 2016-2017 and 2017-2018 are also demonstrated by the ERA5 reanalysis precipitation data in Figure 2C. The years 2017 and 2018 are similar in terms of precipitation and temperature. The year 2019 deviates from the previous 2 years, being considerably warmer and dryer than the previous years, both in summer and winter. This is demonstrated by the freezing and thawing degree days indices and cumulative winter and summer precipitation derived from the reanalysis data ( Figure 2D).
The Berelegh river is an approximately west-east running tributary to the lower Indigirka river. The Berelegh is an anastomosing, multi-channel river with highly sinuous channels and wide marshy floodbasins with lakes. To the north of the floodplain, the landscape is dominated by Late Pleistocene-Holocene aged thaw lakes and drained thaw lake basins (DTLB's), and plateau-like erosion remnants of the Pleistocene Yedoma deposits. In the DTLB's, pingos are common. The DLTB relief is subdued with an elevation in the order of 10-12 m ASL; the DTLB floors are only a few meters higher than the floodplain surface (8-10 m), whilst the Yedoma plateaus rise some 10-25 m above the DTLB floors to elevations up to 35 m ASL. To the south of the Berelegh floodplain, the elevation of the Yedoma surface is higher (40-60 m ASL; Allaikhovsky Highland), and the area of lakes and DTLB is comparatively smaller (Figure 3; Schirrmeister et al., 2011).
The area is largely covered with Pliocene-Pleistocene deposits. At Chokurdagh and on the east side of the Indigirka river Mesozoic basaltic rock and Paleogene-Neogene sandstone crops out. In the Dzhelon-Sisi upland west of the study site, the Eocene Tastakh Formation consisting of clays and sands with coal lenses and gravel is exposed (Schirrmeister et al., 2011;Ren et al., 2013;Akhmetiev, 2015;Drachev, 2016). Below the Yedoma deposits, the Plio-Pleistocene Keremesit and Olyor Formations occur, which consist of fluvial and lacustrine silts and sands, with similar ice-rich characteristics as the Yedoma (Shmelev et al., 2017).
Tectonically, the area lies within the interaction zone of the Eurasian and North American plates (Imaeva et al., 2016a;Imaeva et al., 2016b;Drachev, 2016). The area is situated on the North American plate, consisting here of terranes of various origin. The Cenozoic sedimentary cover is bordered to the south by the Polousnyi mountain ridge, which is an active thrust belt with an WNW-ESE strike. To the SW of the area, the Chersky seismic belt is located (Fujita et al., 2009). Its seismic activity reaches to just south of the study area, with several reports of small (magnitude <4) earthquakes (Imaeva et al., 2016a;Imaeva et al., 2016b). To the north of the area, the intensely deformed Late Jurassic-Early Cretaceous South Anyui suture zone is buried beneath the Cenozoic sedimentary cover (Drachev, 2016).

Data
The SAR dataset used in this study was acquired by Sentinel-1 A/B C-band (∼5.7 cm wavelength) satellite sensor, covering the seasonal time period from June to September for 2017-2019 years. We used 30 Sentinel-1A/B IW SLC SAR images acquired in descending orbits with an incidence angle of 37-39°and a pixel spacing of 2.3 × 14.1 m (range x azimuth).
A Digital Elevation Model (DEM) from TanDEM-X mission is used for InSAR processing, at 90 m spatial resolution. TanDEM-X (TerraSAR-X add-on for Digital Elevation Measurements) is an Earth observation radar mission that consists of a SAR interferometer built by two almost identical satellites flying in close formation. With a typical separation between the satellites of 120-500 m a global Digital Elevation Model (DEM) was generated (Wessel et al., 2018).
In addition a high resolution global DEM from ALOS PRISM (Tadono et al., 2014;Takaku at al., 2016) with spatial resolution of 5 m is used for geomorphology analyses.

InSAR Processing
The InSAR technique exploits multiple SAR images and applies appropriate data processing and analysis procedures to separate the contribution of the phase caused by the deformation, from the other phase components. The technique focuses on the identification of pixels in the SAR image characterized by small noise, which have typically two properties: the radar response is dominated by a strong reflecting object, and remains constant over time (Persistent Scatterer, PS) (Ferretti at al., 2001).
The SAR imagery is processed using ESA's Sentinel Application Platform (SNAP) S1-toolbox (Veci et al., 2014) and Stanford Method for Persistent Scatterers (StaMPS) by the University of Leeds (Hooper, 2008). We used statistical-cost, network-flow phase-unwrapping algorithm (SNAPHU) for unwrapping of the interferometric phase (Chen and Zebker, 2002). The topographic phase is removed using TanDEM-X 90 m DEM (Wessel et al., 2018).
The standard deviation of the mean line-of-sight InSAR velocity is estimated as 4-5 mm/yr. See Supplementary Information for further details on error sources and data processing. A drawback of the InSAR data is the uneven spatial distribution of data points, as is determined by the persistent scatterer algorithm in InSAR processing, which requires stable surface properties over time. This introduces a location bias. Winter data (with variable snow cover thickness) are not included, and the dynamic floodplain (ARF) with varying water cover has a low data point density.
We concentrate in our analysis on the thaw season surface displacement processes. However, summer subsidence may be compensated in winter by heave (Supplementary Figures S1, S2), as indicated by positive surface displacement between the last images of 2017 and 2018 and the first images of the next year. This is attributed to spatially varying ice accumulation in the active layer during the winter.

Soil Ice Content and Active Layer Thickness
Ice content data have been collected by Weiss et al. (2016) in 2012, and in 2013 by the authors. The samples were collected by hammering a steel pipe into the top of permafrost, and carefully extruding the sample in 5-10 cm increments. The samples are weighed and dried, to determine the gravimetric ice content. The ALT at selected sites was measured by probing with a thin steel rod (0.8 cm thickness). Gravelly or stony material is absent in the study area.

Geospatial Analysis
Geospatial data analysis was done with Quantum GIS version 2.18 (QGIS.org, 2016). A landform and surface micro-relief map is produced based on a synthesis of all satellite images, DEM data and field surveys made over the period 2004-2018. The main landform classes are: Yedoma remnant (YR hereafter), drained thaw lake basins (DTLB), active river floodplain (ARF) and river terrace (RT). Table 1 shows further subdivisions in geomorphological details and micro-relief (mainly ice-wedge polygons and their morphology and degree of visibility; for slopes also erosion features such as rills).
Pond density classes were visually estimated. For the DTLB class, an estimate of the relative age was made based on overlapping relations, in combination with ice wedge polygon and surface micro-relief development, which can be used as indicator of age (Hinkel et al., 2003). Based on image time series, flooding frequency classes were derived (permanently flooded, seasonally flooded, extreme floods only and no flooding). ARF flooding frequency was checked by a longer time series of lower resolution Landsat images (Cheung, 2019). The visible spectrum images were also used to derive the normalized difference water index (NDWI, McFeeters, 1996) and normalized difference vegetation index NDVI (near infrared-red)/near infrared+red). NDWI (green + near infrared)/(green−near infrared), indicates the amount of open water within a pixel, NDVI indicates the amount of green biomass. The NDWI data of the 2015 image, with average wetness conditions and completely cloud-free, were combined with the landform data.

Statistical Analysis of InSAR and Terrain Data
The InSAR dataset contains both spatial and temporal variation, and interaction between these sources of variation. The analysis is based on the assumption that spatial variation is related to terrain characteristics expressed by the geomorphological map data, and the temporal variation is caused by the weather patterns throughout the observation years, affecting loss or gain of ice, water or sediment. Interaction between spatial and temporal variation is likely since the weather conditions may influence terrain units differently; e.g. by variation in soil drainage conditions that affect soil water content and thermal properties of the soil. To explore these relations, the InSAR data have been compared with the categorical scale data of the geomorphological map and ratio scale data derived from DEMs and visible spectrum images.
To characterize temporal variability of each data point, we performed a hierarchical cluster analysis (Euclidian distance, Ward's method, divisive algorithm), which groups all data points into dissimilar groups of displacement trajectories, based on the vertical displacements.
Next we applied two pathways of analyzing the relation between terrain characteristics and vertical displacements. The first pathway considers the interval/ratio scale data from the ALOS DEM and visible spectrum Worldview 2 image of 2015 separately. These data are elevation (ELEV), panchromatic greyscale (GS), NDVI and NDWI derived from the Worldview image. ELEV represents drainage conditions and geomorphology, NDVI and GS represent vegetation biomass and type. High values of NDWI indicate areas with much ponded water relative to vegetation, low GS also indicates presence of ponds. First we calculated the Pearson correlation coefficients between these terrain units and respectively the total summer displacement over the 3 years, and the total for each year separately. Next we applied principal component analysis on this dataset (mean and variance standardized) to explore which combination of properties explain most of the variance, and for dimensionality reduction.
The second pathway includes the categorical scale data in the geomorphological map. For relations of surface displacement with the categorical data of the geomorphological map we considered as potential driver variables: • the large landform units (Drained thaw lake basins, river floodplain and terrace, Yedoma) • erosion, sedimentation or ground ice thaw phenomena • pond density • river flooding classes A one-way analysis of variance (ANOVA) for each combination of yearly InSAR data and the categorical classes Landform, EST, Pond density and Flooding) was applied, to detect which of the these categories are linked to surface displacements. A Tukey-Kramer multi-comparison was used to identify classes that differ significantly from others.
Next, the numerical and categorical data selected by PCA and ANOVA are combined into a generalized linear regression model which was fitted stepwise to the surface displacement data of each year separately, and to the summed surface displacement over all years (Σ years). Interaction terms were not included. The regression model gives information on the strength of relationships.

Landforms
The area is dominated by overlapping DTLB (54.5%) and the recent floodplain (ARF, 39.1%). Yedoma remnants (YR, 3.7%) and the river terrace (RT, 2.8%) make up a small fraction of the area. The YR are plateau-like features with a nearly flat upper surface, likely representing an older land surface. An isolated hill consisting of Yedoma deposits testifies former expansion of the Yedoma surface in the northern part of the study area (1 in Figure 4). This old surface slopes gently towards the river (Supplementary Figure S3). The Yedoma deposits have a Late Pleistocene age; calibrated 14 C ages of organic matter and bone at a lake bank exposure are 31,745 and 20,829 years BP, ages of TABLE 1 | Subdivisions of the main geomorphological map units into subunits, and micro-relief classification for these units.

Main unit Geomorphological subunits Micro-relief classes
Yedoma remnant plateau, slope, erosion gully low/high centre polygons, obscured polygons, rill erosion, alluvial fan, colluvial deposition, ponds Drained thaw lake basin basin, shore platform, beach ridge, delta, active lake, pingo, and drainage channel low/high centre polygons, obscured polygons, rill erosion, ponds River floodplain channel, point bar, levee, crevasse splay, flood basin, lake, pingo low/high centre polygons, obscured polygons, ridge and swale topography River terrace terrace plateau, slope low/high centre polygons, obscured polygons, rill erosion, ponds cryoturbated organic matter in the top of the succession are 1,674 and 11,449 years BP (Weiss et al., 2016). The youngest date likely contains young soil material. The YR surface is generally well-drained, but on larger plateau areas low center polygons with intermittent ponded water occur. Tussock tundra vegetation dominates; surface drainage occurs through sedge-lined diffuse water courses, parallel to each other. On steeper slopes these grade into small rills between tussocks, and eventually in gullies where these rills concentrate.
The RT has similar surface characteristics as the YR. The RT may be the river valley equivalent of the YR, although traces of channel systems are not visible. In one location (2B in Figure 4; Supplementary Figure S3B), lateral continuity is suggested by a gradual slope from Yedoma plateau to terrace. Dating information from RT's is not available; a younger (Holocene) age than the Yedoma cannot be excluded. The RT is likely underlain by ice-rich permafrost, as shown by erosion features at RT edges.
The DTLB surface has a 10-20 m elevation difference with the YR surface. It is covered by overlapping drained lake basins. Differences in elevation up to 3 m are common, representing sedimentary features of the lakes, and elevation differences between individual basins. The basins differ in surface microrelief, ranging from well-defined low center polygons in younger basins, to areas where ice wedge polygons are obscured by peat formation in older basins (Teltewskoi et al., 2016;Weiss et al., 2016). In the latter case ice wedges are visible by the formation of linear thaw ponds. A pattern that appears to develop over time on slightly sloping surfaces is a ribbed pattern with parallel, diffuse drainage channels with a dense sedge vegetation, alternating with ±0.5-1 m higher ridges covered by Betula nana heath ( Figure 5). Well-drained areas are dominated by tussock tundra. Features indicating ground ice thaw (thaw ponds) are ubiquitous in DTLB's ( Figures 5, 6). Furthermore, Sphagnum peat growth is widespread, indicated on panchromatic images by light tones.
Individual basins are lined with a low ridge marking their former banks, interpreted as a beach ridge or ice shove ridge built up by material bulldozed ashore by a floating lake ice cover ( Figure 6). The elevation difference between individual basins is small, 1 m or less. Figure 6 also shows a fluvial delta, built into a lake by a channel connecting the lake with the river floodplain. At present, drainage through this channel has reversed, with flow towards the river, likely caused by river downcutting. Other relief elements are slightly elevated shore platforms that line the lake banks on the lake side, and a central bulge of the lake floor, shown as a better drained area in the center of the lake basin.
The overlapping relations of the lake basins were reconstructed, based on polygon maturity (Hinkel et al., 2003), surface elevation, and intersection of beach ridges (Supplementary Figure S4). Only from the older lake bed depicted in Figures 5, 6 14 C dates are available (Schirrmeister et al., 2011;Weiss et al., 2016). Weiss et al. report a calibrated age of 3,251 years BP from a buried peat bed, and ages of 2,884 and 8,016 years BP for cryoturbated peat. Schirrmeister et al. (2013) report uncalibrated ages of 1,632 ± 32 and 2,144 ± 33 from peat. The age of 8,016 years BP indicates that the lake was drained in the Early Holocene. Renewed clastic sedimentation led to the burial of a peat bed after 3,251 years BP. Teltewskoi et al. (2016) report a calibrated basal 14 C age of 3,630 ± 40 for the surface of this lake bed.
Within the DTLB unit, lakes are common. Some of these are actively expanding thaw lakes (e.g. at 5 in Figure 4). Other lakes are remnants after partial lake drainage, or infilling of lakes with peat or clastic sediment (e.g. at 6 in Figure 4). These remnant lakes also may start to expand again, which occurred at 7 in Figure 4. New lakes are created by surface subsidence and merging of polygon ponds and thaw ponds; these are generally smaller and shallower lakes or ponds and have an irregular shape.
Pingos (bulgannyakhs) of all sizes occur abundantly on the drained lake basins (Supplementary Figure S4). Some consist of irregular shaped complexes. The pingos are of the hydrostatic type (Mackay, 1998). A field survey of pingos in the central section of the study area did not show any water seepage activity associated with these pingos, suggesting low or absent activity (Mackay, 1998).
Actively eroding slopes between YR and DTLB occur only near expanding lakes (5 in Figure 4); erosion is dominated by thaw slumps (e.g. at 1 and 5 in Figure 4). When lake bank erosion stops, slump activity continues for some time; at 1 in Figure 4 lobes of thaw slump debris overriding the adjacent DTLB are found. Older slopes have lost their slump morphology by rill erosion and colluviation. Erosion occurs by rills between vegetation hummocks, along ice wedges and occasionally gully erosion. On lower slopes, colluviation and deposition of small alluvial fans likely dominates.
The ARF is separated by clear terrace scarps from the RT and often also from the adjacent DTLB surface; those separating ARF and RT are highest. This indicates two phases of downcutting of the river, after formation of the RT surface and after formation of the oldest DTLB's.
The ARF has a large proportion of lakes. By contrast to the lakes in DTLB, the water level and size of most of these lakes varies with the flood level of the river (e.g. Figure 7, Cheung, 2019). Many lakes likely have an origin as thaw lakes which have become connected with the river (1, 3 in Figure 7), or as thaw lakes eroding the RT unit (e.g. 2 in Figure 7). Several deeper lakes maintain a stable water level at low river stands. These lakes are not connected with outlet channels to the river channel, but nevertheless receive river water during high flood stages.
The river channel pattern is characterized as a high sinuosity river with multiple channels, but dominated by one main channel (Figure 8). At low stage, the non-lake floodplain surface in FIGURE 6 | Border between two lake basins with different age and surface micro-relief. The elevation difference is less than 1 m. The border of the younger lake is marked by a low ridge. In the center a former river delta. Inset: beach ridge of the same lake seen from point 1 in Figure 4; in the lower left corner the beach ridge is overridden by a thaw slump debris lobe, photo J. vanHuissteden. Location of the image frame: see Figure 4 at nr 4. WorldView © 2019 MAXAR.
between the channels consists of mostly vegetated, non-flooded flood basins, dominated by grasses and sedges, and channel levees with willow brushes. A part of these flood basins consist of pointbars with ridge and swale topography created by river bend expansion. The vegetation is denser and more productive than on similar flat areas outside the floodplain (DTLB, YR plateaus). The lowest floodplain levels (channel banks, drained lakes) are unvegetated. Ice wedge polygons are ubiquitous. At low stage in winter, the channel ice surface is some 2-3 m below the floodplain surface, exposing steep channel banks. These banks are sources of aeolian dust, as shown by streaks of dust over the winter snow surface.
Active sedimentation occurs on floodplain surfaces on levees and point bars, but also by suspended sediment and organic debris deposition in flood basins. Sediment cores in flood basins show successions of vegetation horizons, alternating with layers of fine-grained clastic material. In a core in flood basin sediments taken on a crevasse splay, a sedimentation rate of 0.48 ± 0.10 cm yr −1 was determined by radiocarbon dating of a vegetation horizon (185 ± 30 years, Groningen lab. nr. GrA59784, at depth of 90-95 cm). This sedimentation rate is corrected for a 15 vol. % excess ice content below the 0.5 m thick active layer.

Ice Content and Active Layer Depth
The ice content at the top of the permafrost varies strongly over short distances. The average ice content for DTLB is 42.60 ± 20.32% (n 34), for YR 38.5 ± 18.46% (n 12) and for ARF 26.45 ± 19.45% (n 12). The YR and DTLB have approximately equal topsoil ice content. The ice content of ARF is lower than that of YR and DTLB; it differs significantly from that of DTLB (Tukey-Kramer test, p 0.04). Fresh YR exposures show large ground ice masses in ice wedges. A 7 m deep borehole in DTLB sediments near the research station showed a decrease of ground ice after the top 3 m (G. Iwahana, pers. comm.).
On 19-20 July 2018, three transects of midsummer active layer thickness (ALT) were measured perpendicular to the maximum flood lines in 2017 and 2018 (Figure 9), and compared with ALT measurements in grids outside the flooded area. Transect 1-3 are located entirely on usually mesic tussock tundra, transect 4 in a wet area dominated by hummocky Sphagnum and Salix vegetation and a thicker surface organic horizon. The ALT for transect 1 to 3 show a significant increase at the flood limits (t-test, p < 0.01); this thickness increase is absent in the fourth transect. The grids a (100 × 100 m) and b (10 × 10 m) are active layer monitoring sites; grid a is located outside the flood limits (vegetation: tussock tundra and Betula nana heath), grid b (Sphagnum and Salix) inside the flood limits. The 100 point average of these grids, measured in the same week as the transects, are shown alongside the transect data. The average ALT of grid a is lower than that of transect 1-3 outside the flood limit; that of grid b is similar to transect 4.

Statistical Analysis of Velocity Data
The cumulative surface displacement over the thaw seasons of 2017-2019 are shown in Figure 10. There are clear regions where data points show similar displacement trends, either positive or negative. These displacement patterns vary among the years (Supplementary Figure S5). A hierarchical cluster analysis (Euclidian distance, Ward's method) using the displacements detected from the successive SAR images is used to cluster the data points into 11 groups of trajectories. In Figure 11, most clusters also represent fairly tight spatial clusters.
To evaluate relations of surface displacement with quantitative parameters of geomorphology, drainage and vegetation we have compared the total surface displacement over 3 years and the summed surface displacements over each year, with surface cover FIGURE 7 | Open water on the river floodplain at different river stages. Backdrop: geomorphological map units, see Figure 4 for legend. Top: low stage, end of summer; mid: normal, ± bankfull river stage, below: extremely high flood, exceeding the normal floodplain limits. 1: river-connected thaw lakes; 2: thaw lake developed on floodplain in river terrace; 3: locations where connections between thaw lake and river is developing during high flood. WorldView © 2019 MAXAR. data derived from the Worldview image of 2015 (NDWI, Panchromatic greyscale (GS), NDVI) and the ALOS DEM. GS is related to both vegetation type and water ponding, NDVI reflects vegetation, NDWI water ponding. Table 2 shows the correlations of surface displacement with the quantitative image data and DEM. NDWI, GS and NDVI correlate with surface displacement, but not equally among years and with varying signs. The surface displacement in 2017 correlates weakly negative with NDVI and NDWI; for 2018 there are no significant correlations; for 2019, positive correlations are found. Because of their relation with ponded water, GS -NDWI, and NDVI -NDWI correlate negatively.
Elevation has a strong positive correlation with surface displacement in 2017 (higher areas -less subsidence or more heave), but negative in 2019 (higher areas-more subsidence). Although there are significant correlations of displacement with vegetation-related parameters NDVI, greyscale and NDWI, these are weak in comparison to elevation. It indicates poor correlation with vegetation-related parameters and surface displacement, although this may have been affected to some extent by differences in spatial resolution between the InSAR and NDVI/NDWI data (see Data and Methods). Correlations of surface displacement between years have different signs; 2018 deviates, with negative correlation between the previous year and the next. These spatially different patterns of displacement in successive years is also visible in the trajectory clusters of Figure 11.  A principal component analysis on the standardized data shows that three components explain 95% of the total variance. The first principal component explains 60% of the variance, with a strong loading on the total surface displacement of 2019 (Supplementary Figure S6). The second principal component contains 27% of the variance, with the highest loading on the 2017 data. From the other variables, only elevation determines a small part of the total variance, in the third principal component (8% of total). This indicates that variability of the surface displacement is dominated by variability between years, related to elevation, but not to vegetation or drainage.
The categorical classes of the geomorphological tested in the further analysis are: The ANOVA to detect the effect of the categorical classes of the geomorphological map (Landform, EST (erosion, sedimentation or ground ice thaw), POND and FLOOD) yields significant differences at the p < 0.01 level ( Table 3, Supplementary Figure S7). A Tukey-Kramer multicomparison was used to identify classes that differ significantly from others. The multi-comparison tests of class differences indicate significant differences (p < 0.01) for most of the classes, except for the EST and POND, where several classes do not differ significantly. In Table 3, only those classes have been listed that differ significantly from other classes within the same variable, and depart more than the uncertainty range (3 mm) from 0 mm displacement.
The categorical variables from the geomorphological map explain a significant amount of the variability of the InSAR surface displacement. However, the behavior of each class is not consistent between years. For the landform classes, the floodplain and river terrace (ARF, RT) shows the most consistent behavior with predominantly subsidence throughout all years. For EST, classes that indicate net sedimentation in nonflooded areas (colluvial deposition and peat growth) show clear, and among the years consistent, net positive surface displacement. However, sedimentary environments on the    river floodplain or DTLB areas show net subsidence. Areas with inferred ground ice thaw (indicated by ice wedge thaw ponds) do not show subsidence, and even positive displacement in 2019. Pond density is the least predictive, only the "very low" density class (1 in Table 3 A generalized linear regression model ( Table 4) was fitted stepwise to the surface displacement data of each year and the summed surface displacement over all years (Σ years), with Elevation (ELEV), and the categorical geomorphological map classes listed in Table 4 as predictor variables. Elevation was included because of its contribution to the variance in the principal component analysis. Interaction terms were not included. All models capture a significant amount of the variance in the data (p < 0.01, F test).
For 2017 and 2019, all variables are retained by Akaike's (Akaike, 1974) criterion; for 2019, ELEV was removed from the model ( Table 2). For Σyears, the model is reduced to LANDFORM, EST and FLOODING variables; for 2017, the positive coefficient for ELEV is significant.
Within the Landform group, the DTLB class was taken as the reference variable. For 2017, the coefficient for LANDFORM-YR is positive, in agreement with the positive coefficient of ELEV in that year and the higher elevation of YR. In 2019 however, the YR term has a strongly negative coefficient. Of the other LANDFORM groups, the ARF and RT have negative coefficients, indicating subsidence relative to DTLB throughout all years.
The EST variable includes processes that can result either positive or negative surface displacements. The widely occurring class with evidence of ground ice thaw (mostly ice wedge thaw) was taken as reference here. Colluvial deposition on slopes and peat growth show pronounced positive coefficients for 2017 and 2018, in line with expected positive displacement for these classes.

DISCUSSION
Geomorphology and Late Quaternary Development Figure 12 summarizes the landscape development of our study area that we deduced from our geomorphological and geological data. Drained thaw lake basins (55%) and the present river floodplain of the Berelegh river (39%) dominate the study area. We assume that the ice-rich Yedoma deposits have had a much wider extension but have been reworked by ground ice thaw and erosion. This is based on an erosion remnant (at 1 in Figure 4), and presence of Yedoma hills north of the study area. The original Yedoma surface had an elevation of 30-38 m above present sea level in the study area, sloping gently towards the river (Supplementary Figure S3). The erosion remnants at larger distance of 2 km do not indicate higher levels than 38 m, suggesting an essentially flat surface at distance from the river. The DTLB surfaces have an elevation between ±10 and ±14 m, indicating the removal of atleast 16 m of Yedoma deposits and ground ice over a large area, and subsequent redistribution in fluvial and lacustrine sediment. Remarkably, most of the Yedoma remnants occur at short distance from the river. A possible explanation that emerged from landscape-scale thaw lake modeling, is that expansion of thaw lakes close to rivers is restricted since the probability of early lake drainage is higher near the river, protecting Yedoma to large scale thaw lake erosion (Van Huissteden et al., 2011).
We assume that the RT unit is of the same age as the Yedoma surface. Its surface is higher by several meters (up to 20 m) than the highest DTLB surfaces (up to 14 m). Former fluvial channels are invisible on the terrace remnants, but may be obliterated by surface processes (sedimentation of loess, overbank river sediments, peat, or erosion).
Thaw lake formation started between at least 20,829 years BP (youngest 14 C date on Yedoma deposits) and 8,016 years BP (oldest 14 C on drained thaw lake basin fill) in the study area. In a recent analysis of Arctic-wide lake sediment basal dates (Brosius et al., 2021) two peaks in northern lake formation, including thaw lakes, were detected at 13,200 and 10,400 years ago, both following rapid increases in North Atlantic air temperature. This agrees well with the radiocarbon dates in our study area. The lake basin development in East Siberia is assumed to last until the end of the Boreal period (9.0-7.5 kyr BP, Romanovskii et al., 2004;Kaplina, 2009;Morgenstern et al., 2013). Figures 4, 8 indicate the formation of large, elliptical to eggshaped lake basin with smooth shores. This indicates the operation of wind-driven longshore currents during lake formation , and a tendency to elongation in southerly direction. The overlapping relations of lake basins (Supplementary Figure S4) show that lake basin formation was well underway before the large, prominent lake basins in the center of Figure 4 have been formed. The lake basin at 3 in Figure 4 was drained before 8,016 years BP, but overlaps older lake bottoms (Supplementary Figure S4, right side of the basin). This suggests an early, Late Glacial start of widespread lake formation, similar to the East Siberian area described by Morgenstern et al. (2013). After drainage of this lake, thaw lake activity shifted to the north with the formation of a new large lake basin at a slightly lower level (at 4 in Figure 4), indicating that the development of lakes continued after ca 8,000 years BP. This lake received river sediment, as shown by the presence of a delta (Figure 7). The presence of an elevated central bulge in many lake basins is attributed to a longer duration of sediment accumulation in the centre of expanding lakes, compared to near-shore areas where subsidence and sedimentation began later, as demonstrated by modeling by Kessler et al. (2012).
Present actively expanding thaw lakes are small compared to these older drained lake basins. They are restricted to the older basins floors, but may expand rapidly into remaining Yedoma (at 5 in Figure 4). At 7 in Figure 4, a remnant lake, left after lake drainage, appears to have been rejuvenated and is expanding again. Thaw lake formation and expansion needs large amounts of ground ice in the subsurface (Jorgenson and Shur, 2007;Morgenstern et al., 2008;Morgenstern et al., 2013), and therefore old thaw lake basin floors are less susceptible to rapid lake expansion. In the older lake basin floors significant accumulation of ice wedges and ice-rich organic deposits has occurred (Teltewskoi et al., 2016;Weiss et al., 2016) but most of the ground ice volume is restricted to the upper meters as shown by deeper boreholes (See Results-Ice Content and Active Layer Depth).
After major thaw lake expansion, a small 1-2 m high terrace scarp was created between ARF and DTLB, likely by fluvial downcutting. The downcutting occurred after the formation of the DTLB at 4 in Figure 4. Before drainage of this basin, a delta existed in this lake, with sediment supply from the river towards the lake through a feeder channel ( Figure 5). Presently, the drainage is reversed towards the river. Fluvial downcutting should have facilitated this reversal, and subsequent lake drainage and the northward extension of the channel into the drained basin. The age of this river downcutting is assumed to be younger than Mid Holocene age, based on the age of the DTLB's.
On the ARF, several lakes have a morphology that suggest ground ice thaw as mechanism for their expansion: steep banks and rounded shapes, in particular bordering the RT unit. At 8 in Figure 4, thaw lakes developed within the river terrace. Some of these lakes are drained by the river channel during low flows (Figure 7, top), but others maintain a drainage threshold with the channel. Within the DTLB unit, lakes occur which appear to be thaw lakes that have been connected with the river quite recently, and drain during low stands (e.g. at 2A and 2B in Figure 4; see also Figure 7 at 1 and 2). During extreme floods, also lakes at larger distance from the floodplain have become connected to the river (at 3 in Figure 7 bottom).
On non-permafrost floodplains, lateral erosion is the primary mechanism for floodplain expansion (e.g. Howard, 1996). However, in our study area the capture of thaw lakes by subsidence and erosion of drainage thresholds may be induced by ground ice thaw during floods, representing a mechanism of floodplain expansion (see discussion of InSAR data). Capture of thaw lakes with subsequent avulsion of river channels through the former lakes may also be the cause of the irregular, patchy distribution of the river terrace (RT) unit on the floodplain, rather than the more common location of river terraces on valley sides. The height of the threshold between the lakes and the river system strongly affects river water input to the lakes, and as a consequence sediment and nutrient supply to the lake water (Marsh et al., 1999;Emmerton et al., 2007). Once connected to the floodplain, sedimentary infilling of these lakes may accelerate.
The compressional tectonic stress in the area has resulted in fault movements. In the DEM of Figure 3, two clear W-E trending lineaments can be discerned: the steep and nearly straight southern limit of the Berelegh floodplain, and a further terrain step of about 10 m to the south of it. Furthermore, there is a marked elevation difference between the Yedoma plateaus in the Allaikhovsky Highland (±50 m) to the south of the Berelegh river and the lower elevation of Yedoma plateau surfaces (30-38 m) north of the river. This suggests neotectonic movement along these lineaments, with a downward warping of the floodplain. On a North-South profile, the InSAR data indicate an average velocity of −5.3 ± 3.1 mm yr −1 north of the lineament bordering the Berelegh floodplain, south of this lineament −0.17 ± 3.2 mm yr −1 , indicating a relative displacement of −3.6 mm yr −1 . The elevation differences of 12-20 m, likely developed since the Late Glacial, results in an estimate of the tectonic subsidence rate of 1-2 mm yr −1 . Further lineaments, with a less clear terrain expression, occur north of the Berelegh valley (Figure 3), but apparently did not result in a large amount of displacement.
Tectonic subsidence may have contributed to the extensive lake formation in our study area by impeding drainage and creating more flood-prone areas. Jorgenson and Shur (2007) invoke a large scale water balance change at the Last Glacial termination to explain the rapid expansion of large thaw lakes. In a poorly drained subsiding basin this should have had larger impact than in an area without subsidence. Tectonic movements resulting in base level changes, as well as climate change resulting in the change of the drainage basin water balance and sediment supply are likely causes of river downcutting, and may have acted in concert. For the prominent downcutting phase that formed the RT river terrace water balance changes at the Last Glacial Termination are a likely cause, next to tectonism. Jorgenson and Shur (2007) cite evidence for water balance changes at the LGT that also initiated thaw lake formation in Alaska. Another mechanism of floodplain erosion is catastrophic drainage of large thaw lakes, which could create peak flows larger than normal spring floods (Marsh and Neumann, 2001;Marsh et al., 2009;Arp et al., 2020). This likely occurred repeatedly in the drainage basin of the Berelegh river given the large DLTB areas north of the river.

Surface Displacements Detected by InSAR
We did not consider a tectonic component in the InSAR data analysis over the study area of Figure 4, which is entirely north of the most active lineaments. The principal component analysis shows that variation among years explains 87% of the variance in the first two principal components (Supplementary Figure S6). The large contrast between the weather conditions in 2019 and 2017-2018 ( Figure 2) causes most of the variability. Terrain variable Elevation appears in the third principal component, which explains 8% of the variance while the contribution of vegetation (NDVI) or its ponded water component (NDWI) is negligible.
This is likely the result of variation in weather conditions. The year 2019 was exceptionally warm and dry (Figure 2), with a warm winter, an early spring and warm summer. The year 2017 and 2018 had high precipitation during the preceding winter. The net summer precipitation in 2017-2018 showed a small precipitation deficit, without long periods of warming in summer. By contrast, 2019 had a large summer precipitation deficit, and prolonged periods of hot weather. The dominance of the surface displacements in 2019 in the first principal component (60% variance), indicate that in particular warm and dry weather conditions induce surface subsidence on a large scale.
Although the yearly weather conditions dominate, the trajectory clusters in Figure 11 show rather consistent spatial clustering, indicating that the terrain conditions contribute significantly to the spatial variation in surface displacements. This is also confirmed by the ANOVA and regression analysis. Weather conditions interact with terrain conditions. Several groups of weather and climate related processes could contribute to surface displacements: 1. Thaw or aggradation of ground ice (Liu et al., 2010;Liu et al., 2012;Liu et al., 2014), induced by changes in heat transport into the soil, by summer heat, a thick snow cover, change of surface vegetation, forest fire, flooding; e.g. Nauta et al., 2014;Liu et al., 2014;Jones et al., 2015; references in van Huissteden, 2020). 2. Shrinkage or swelling of soils on decrease or increase of soil water content. In particular fine-grained soils and peat are sensitive (e.g. Camporese et al., 2006;Succow and Joosten, 2012); this shrinkage/swelling can be detected by InSAR (Cuenca, 2008). 3. Surface displacements resulting from peat decomposition have been detected by InSAR (Cuenca, 2008). Decomposition of peat is affected by soil moisture and temperature (e.g. Succow and Joosten, 2012). The reverse process, peat accumulation, increases surface elevation. Sphagnum peat growth occurs widely in the study area, in particular within the DTLB unit (Magnússon et al., 2020). 4. Erosion and sedimentation are triggered by increased runoff and flooding; in permafrost areas the release of water by ground ice thaw on slopes triggers thaw slumps and active layer slides (Lewkowicz, 2007;Lewkowicz and Way, 2019).
The summer warmth and dryness in 2019 has contributed to significant surface subsidence compared to the previous years, for several of the clusters of data points in Figure 11, notably clusters. 1, 2 and 7. The regression analysis (Table 3) shows that in 2019 most categorical terrain variables have negative coefficients, indicating dominant subsidence. Most terrain classes show smaller positive coefficients, more negative coefficients, or reversals from positive to negative, notably the Yedoma landform class. A notable exception is Pond density class 4, of which the regression coefficient turns significantly positive in 2019, from negative in 2017-2018.
Both loss of ground ice and shrinkage of peat soils could have contributed to the pronounced surface subsidence in 2019. Trajectory clusters 1, 2 and to a lesser extent 7 all show pronounced subsidence in 2019 compared to 2017 and 2018. Clusters 1 and 2 are located in the large DLTB on the center-north of the study area (Figure 11), 7 in a neighboring area immediately to the south in an older DTLB. Both are areas with peat growth. Cluster 1 and 2 show low center polygons, cluster 7 represents an older DTLB area where polygons have been overgrown with peat ( Figure 6). Shrinkage of peat during the dry summer therefore could have contributed to subsidence. However, small perennial frost mounds (without peat growth) in these clusters also show subsidence, indicating that ground ice thaw has occurred also.
Positive surface displacements in 2019 occur in cluster 6 and 10. Most of these points are areas with intermediate -high pond density and ice wedge polygons (cluster 6), or areas on the lower slopes of Yedoma remnants and higher ground in DTLB's (cluster 10). The lower slopes of the YR may have been subject to colluvial deposition or downslope creep of soil material by increased water release from thawing permafrost in 2019. In 2017 and 2018, areas with potential colluvial deposition and peat growth also showed aggradation (Supplementary Figure S5 and Table 3). For areas with ponds in 2019, vegetation and Sphagnum peat growth may have been stimulated by decrease of water depth in ponds (Magnússon et al., 2020).
Positive surface displacements in 2017 correlate significantly with elevation (Table 1), the YR landform (Supplementary Figure  S6 and Table 3) and trajectory cluster 4 ( Figure 11). This is less clear in 2018, and in 2019 elevation correlates negatively with surface displacement. Pronounced positive surface displacement data points in 2017 are located on higher, well-drained parts of the Yedoma surface and DTLB surfaces. These areas have dry soils in summer with a thick active layer, consisting of silt and organic matter. The release of large amounts of meltwater from the thick snowpack after the winter of 2017, followed by a summer with relatively high amounts of precipitation ( Figure 2) likely increased soil moisture. Downward movement of this additional water towards the transition layer at the top of the permafrost contributes to ice lens formation and frost heave in summer (Mackay, 1983). In addition, swelling of the soil by high moisture content could have contributed to the surface heave.
Seasonal flooding may contribute to sedimentation in floodplain environments, reflected by the positive regression coefficients for seasonal flooding in the regression data for all years in Table 4. However, this is probably a spurious effect; the number of points on seasonally flooded floodplain areas is very small.
The areas outside the floodplain that are reached by extreme floods only, show significant subsidence in 2017 and 2019 in the ANOVA results (Table 3), but the regression coefficients for extreme flooding in Table 4 remain just below the p 0.01 significance level. The extreme flooding areas are represented by a low number of points, in clusters 4, 5, 9 and 10. Nevertheless, from Figure 9 a pronounced increase of summer thaw depth of dry soils that were flooded in 2017 and 2018 (transect 1, 2 and 3 in Figure 9) can be inferred at flooded locations. It is unlikely that the thaw depth recovers later in summer. Ponded water and wet soil conditions after flooding result in increased heat flux into the soil by changing the heat balance at the surface and increasing the soil thermal conductivity (e.g. Westermann et al., 2016), creating a thaw depth legacy of the flooding after flood recession. Nearby floodplain areas that are seasonally flooded have a thicker late summer ALT than non-flooded areas (van Huissteden et al., 2005;van Huissteden et al., 2009). Therefore, flooding is expected to increase ALT throughout the summer and fall with loss of ground ice (Figure 11), by heat transfer from flood water into the soil. However, this depends on soil conditions; wet soil with a Sphagnum peat cover was not affected (transect 4 in Figure 9).
Also for the RT landform unit, subsidence is pronounced ( Figures  10, 11). During extreme floods, the RT are small isolated islands on the river floodplain surrounded by water (Figure 7). Prolonged flooding in 2017, 2018 likely adds to heat transfer in this unit.
Given the absence of river gaging data, it is difficult to specify a flood return period for the high floods of 2017 and 2018. However, the low terrace scarp bordering the floodplain, which limits most floods, was exceeded over a large distance in 2017 and 2018. Downstream in the Indigirka river, the return period of the 2017 flood level was estimated at ≥49 years (Tei et al., 2020). A study of Landsat images of the study area during the spring flood period dating back to 1999 did not show any higher flood limits than the 2017 flood (Cheung, 2019).
A higher frequency of extreme floods likely induces subsidence at ice-rich floodplain borders, making these areas in turn more floodprone. Additionally, subsidence of drainage thresholds could facilitate the connection of thaw lakes and DTLB areas to the floodplain. The sites numbered 3 in Figure 7 are critical areas in that respect. For the westernmost point, the drainage threshold was just not exceeded in 2017, for the easternmost point, the flood invaded two lakes in a large DTLB, and a channel connecting these lakes towards the floodplain is forming. We hypothesize that ground ice thaw by an increase of flood extremes will expand floodplains in ice-rich permafrost. A future increase of flood height and flooding frequency therefore may have a large impact on floodplain development in ice-rich permafrost floodplains (Shkolknik et al., 2017).
So far, we have only considered thaw season surface displacement processes. However, as Supplementary Figures  S1, S2 show, surface subsidence is compensated by spatially varying heave in winter. This is attributed to ice lens growth, determined by local water availability in the active layer. Points with strong subsidence in one or more thaw seasons may have gained in height overall years when winter heave is included. Overall, net subsidence dominates in the study area over 2017-2019, but there are marked areas where net heave occurred.
For instance, in the relatively young drained thaw lake basin in the north-central area of the study area, a cluster of points occurs with pronounced subsidence in 2019, resulting in considerable cumulative subsidence over the thaw season ( Figure 11 cluster 1, 2, Figure 10), but the net effect of winter heave and summer subsidence over 3 years is substantial heave. (Supplementary Figure S2). This is an area of active Sphagnum peat growth, with a high soil water content. Also the Yedoma surfaces show net heave over the 3 years, which is strongest in 2017 (cluster 3, 4 and 8 in Figure 11). This is mainly driven by soil water accumulation in 2017 as discussed above.
Pingos that were visited in the study area do not show any recent growth activity. Compared to normal DTLB surfaces, pingos tend to subside a few mm over the years in the thaw season (t-test, p < 0.01) and also on average over the 3 years, which is like the result of slope processes and loss of ice.

Carbon Cycle Implications
The InSAR data presented here do not allow estimates on the carbon loss or greenhouse gas emissions from thawing permafrost, since this would require quantitative data on soil moisture conditions, active layer depth and carbon content of the soil. However, general trends can be inferred: 1. Increased flooding enhances ground ice thaw and active layer depth along floodplain limits and may expand floodplains. This enhances CH 4 emission (van Huissteden et al., 2005;van Huissteden et al., 2009), by exposing old permafrost carbon to anaerobic decomposition, and potentially by creating a more productive wetland vegetation through nutrient supply by flood water. Carbon emission in floodplain environments is fueled by higher ecosystem productivity of floodplains; most of the released fluvial-derived carbon is of recent origin (Dean et al., 2020). However, floodplain subsidence and capture of thaw lakes also creates storage space for burial and sequestration of organic-rich sediments (Emmerton et al., 2007;Van Huissteden et al., 2013;Vonk et al., 2015). In this location a sedimentation rate of organic-rich silt of 0.5 cm year −1 was recorded in a borehole. However, tectonic subsidence also contributes to the sedimentation rate. 2. Outside floodplain areas, warm or wet summer conditions contribute to thaw of ice-rich permafrost, subsidence and enhanced exposure of permafrost organic matter to aerobic or anaerobic decomposition (e.g. Dorrepaal et al., 2009;Natali et al., 2015). However, the net effect on carbon loss and the ratio between CO 2 and CH 4 emission is not straightforward, since multiple interactions on various time scales between vegetation primary production, nutrient liberation by permafrost carbon decomposition and soil water content occur (Natali et al., 2015). Subsidence therefore cannot be used as a proxy for carbon loss. 3. Wet soil conditions enhance subsidence, as shown by increased subsidence in 2017-2018 in areas with intermediate-high pond density. On the other hand, during the warm/dry year the well-drained Yedoma was more prone to subsidence, and areas with high pond density showed positive surface displacement, probably linked to enhanced vegetation growth/peat accumulation (see at 3 below). Thaw ponds in the DTLB unit in this area release mostly contemporary CH 4 and CO 2 carbon from recent vegetation and soil (Nauta et al., 2014;Dean et al., 2020). This indicates that release of old carbon is minor. However, the loss of recent vegetation and soil carbon by pond formation may be substantial, and the release of CH 4 decreases the greenhouse gas sink or increases the net greenhouse gas emission of the area. Greenhouse gas emission of ponds is partly compensated by recolonization of ponds by vegetation (Magnússon et al., 2020); see also previous section). 4. Some of the areas where peat formation by Sphagnum mosses occurs, showed pronounced subsidence in dry and warm 2019. This suggest dryer conditions by shrinkage and/or ground ice thaw. Dry conditions may limit peat growth and its carbon sequestration. 5. Both warm and wet years enhance mass wasting and colluviation on slopes of Yedoma landforms. Mass wasting by thaw slumps, active layer slides and gully erosion exposes permafrost and vegetation carbon to decomposition; sedimentation of mass wasted material and colluviation buries soil and vegetation carbon (Grosse et al., 2011;Abbott and Jones, 2015).

CONCLUSION
In our study area, erosion remnants of a paleo-surface of the Glacial Yedoma landscape testify of the former morphology of this landscape. The Yedoma surface extended with gentle slopes into the present river floodplain, where an old terrace level testifies of the floodplain surface at that time. The fluvial downcutting that created this terrace could relate to tectonic movements and/or change of the water balance of the drainage basin. The formation of thaw lakes during the Late Glacial and Early Holocene eroded most of the Yedoma deposits, leaving only minor erosion remnants near the river floodplain. Probably formation of large thaw lakes was enhanced by tectonic subsidence. The Yedoma remnants were preserved preferentially near the river floodplain, because drainage by fluvial erosion limited the growth of thaw lakes. Catastrophic lake drainage could have contributed to floodplain erosion.
The potential for short-term geomorphological changes is shown by the InSAR surface displacement data. These data, and midsummer active layer thickness data collected after the floods of 2017 and 2018, show that heat transfer into the permafrost by flooding during high flood stages results in subsidence along floodplain margins, lowering drainage thresholds of thaw lakes near the river. We hypothesize that flood-induced subsidence is an important mechanism for floodplain expansion if future climate change increases spring flood heights.
Variability between the observation years contributes most to the variance of the InSAR data, but interacts with the local terrain characteristics. Surface displacement patterns differ considerably between the wet years 2017 and 2018, and the warm and dry year 2019. In general, during the wet years well-drained, elevated terrain appears less sensitive to subsidence by ground ice thaw and even may show heave. During the warm and drier year 2019, overall subsidence is more prominent.
Subsidence or heave are generally caused by changes in the ice volume at the top of the permafrost. Increased soil water availability in wet years contributes to gain of ice at the top of the permafrost, in particular at relatively dry sites. However, variations in soil moisture content also may contribute to surface displacement by shrinking or swelling of peat and fine-grained Yedoma soils. Peat growth and colluvial deposition on lower slopes result in positive surface displacement. However, detailed observations on changes in ice content and sedimentation are necessary to constrain contribution of each of these processes.
The relation between subsidence and the carbon cycle is not simple because of the many soil -vegetation interactions on various time scales. Combined with water ponding, subsidence increases CH 4 emission. Carbon cycle implications of floodplain expansion are the creation of a larger storage space for organic-rich fluvial sedimentation, but also expansion of areas with a high CH 4 emission.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusion of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
JH did the geomorphological interpretation of the satellite images and field data, the statistical analysis and wrote most of the text. KT analyzed the radar images and computed the InSAR data. YC analyzed river flood levels from older Landsat image data, the data were provided by NLR and her work was supervised by HN. RM participated in fieldwork, discussions and provided a part of high resolution satellite image data. SK participated in fieldwork organization and data collection. TM and AD supervised fieldwork organization and contributed to discussion of the research results.