Ground subsidence and polygon development due to thermokarst in the Lena-Aldan interfluve, eastern Siberia, revealed by satellite remote sensing data

Abstract

Thermokarst development is a topographic change in the landscape that is commonly associated with permafrost degradation in ice-rich permafrost regions. The Lena-Aldan interfluvial area in Central Yakutia in eastern Siberia has undergone extensive thermokarst development in the last three decades, particularly in the vicinity of settlements. Despite the negative effects of thermokarst development on the inhabitants of these settlements, no quantitative observation methods have been developed to investigate the surface displacement due to thermokarst development over entire towns. This study utilized interferometric synthetic aperture radar to reveal ground-surface displacement associated with thermokarst near the settlements of selected towns. The findings showed that significant subsidence was detected in disturbed areas (farming and abundant arable land) near the towns. The magnitude of subsidence in the Tyungyulyu and Mayya areas was less than that in Churapcha and Amga. Polygon density in a defined area in each town was examined using high-resolution optical images. The polygon density in Churapcha was considerably lower than that in Mayya, whereas polygonal texture in some areas in Tyungyulyu and Amga was unclear. Spatial frequency analysis using satellite optical images showed clear differences in averaged spectrum models between well-developed and less-developed polygons, which may reflect trough depths and density of vegetation between polygons. Satellite-based subsidence maps and statistics describing polygon development may be useful for evaluating both initial and subsequent stages of thermokarst development.

1 Introduction

Thermokarst development is an irreversible process associated with landscape changes caused by the melting of ground ice (ice wedges) underlain by ice-rich permafrost (Czudek and Demek 1970; van Everdingen 2005). An example of ice-rich permafrost, the Yedoma Ice Complex, is found in eastern Siberia. In particular, the area between the Lena and Aldan Rivers (Fig. 1) is characterized by abundant ground ice with a volumetric ice content above 60% (Schirrmeister et al. 2013; Shestakova et al. 2021). Recent climate changes have promoted permafrost degradation in the interfluve (Iijima et al. 2014), and thermokarst-induced polygonal subsidence patterns, hereafter referred to as polygons, have started to appear in residential areas (Crate et al. 2017). Thermokarst development in the region has had a negative impact on the daily lives of local inhabitants as the changes in the landscape have been associated with extensive destruction of infrastructure (Shiklomanov et al. 2017; Hjort et al. 2018; 2022).

Fig. 1

Map of the Lena-Aldan interfluve in eastern Siberia. Black dots indicate the location of representative towns, and dashed rectangles indicate observation areas of ALOS-2 images

Thermokarst development can be divided into four main stages (Bosikov 1991; Desyatkin et al. 2009). In the initial stage, called byllar, the ground surface starts to subside, and polygonal patterns begin to appear at the surface due to melting of ice wedges. As thermokarst development progresses, the extent of subsidence increases and the bumpy relief become more apparent, often forming small lakes. These first two stages of initiation and enlargement typically occur over a period spanning several decades or up to a few hundred years (Soloviev 1973; Fedorov et al. 2014a). The ice wedge continues to thaw, and the resulting depression is filled with the meltwater from the ice wedge and precipitation, which promote further thawing of the permafrost and the formation of large thermokarst lakes. After water evaporation, stable thermokarst basins, called alases, are formed. Approximately 16,000 alases exist in the lowlands of Central Yakutia (Bosikov 1991), and alases cover approximately 30% of the territory’s surface. Although all the stages of thermokarst development can be observed in the interfluve, thermokarst dynamics are still not fully understood. For example, during the first two stages of thermokarst development, the spatiotemporal evolution of surface subsidence remains unclear. Long-term field observations of thermokarst that have been carried out in Yukechi (Fig. 2a) since the 1990s have revealed inter-annual subsidence (Fedorov and Konstantinov 2008) and thermokarst lake expansion (Ulrich et al. 2017). However, the large scale of the areas that have undergone thermokarst subsidence in the interfluve has meant that investigations are logistically difficult, and conventional remote sensing techniques are not capable of detecting subsidence at a level of several to tens of centimeters per year. As a result, the temporal changes in thermokarst subsidence and prediction of further changes are difficult to assess.

Fig. 2

Study area. a Elevation, towns (stars) and an observation site by Melnikov Permafrost Institute (dot) in the study area close to Lena River; b Map of terraces distinguished by Soloviev (1959); c Elevation along the Amga River around Amga and towns (stars)

Recent advances in remote sensing techniques have enabled us to examine ground-surface changes from space without any equipment in situ. In particular, interferometric synthetic aperture radar (InSAR) well suited for detecting surface displacement (e.g., Hanssen 2001). InSAR time-series analysis using a number of interferograms can be used to derive spatiotemporal changes in surface displacement (e.g., Berardino et al. 2002; Schmidt and Bürgmann 2003; Biggs et al 2007), and several studies on permafrost dynamics using this method have been published to date (e.g., Liu et al. 2015; Daout et al. 2017; Antonova et al. 2018). In eastern Siberia, van Huissteden et al. (2021) used C-band Sentinel-1 InSAR data to derive vertical displacement in an area of the Chokurdakh tundra (70.49°N, 147.29°E) during the thawing season. In the Lena-Aldan interfluve, the average rate of surface subsidence due to permafrost thawing has been examined by L-band InSAR stacking in Mayya and Churapcha (Abe et al. 2020; Iijima et al. 2021; Fig. 2a). Abe et al. (2022) recently examined the inter-annual and seasonal surface displacement of a fire scar in the southern Mayya using L- and C-band InSAR analysis. They reported that L-band InSAR was suitable for monitoring long-term surface displacement in terms of coherence. On the contrary, C-band InSAR with a high temporal resolution is very useful for examining short-term seasonal displacement, but the C-band coherence can sometimes be lost, even in 12 days, particularly in forested taiga areas (Abe et al. 2022). Consequently, L-band InSAR time-series analysis is better suited to deriving long-term surface subsidence due to thermokarst. However, the studied areas were very limited, and the spatiotemporal changes in surface displacement remained unknown.

Polygonal terrain is formed by thermal contraction cracking, characterizing in terrestrial periglacial regions. Especially in ice-rich permafrost regions, high-centered polygons are resulted from the widening and deepening of the polygon-outlining troughs due to ice-wedge degradation (e.g., French 2017). Ice-wedge polygons are the most common type of thermal contraction polygons with 5–20 m in diameter in Arctic regions (Black 1976). Therefore, the size and density of polygons that form in the initial stage of thermokarst reflects the local geomorphological characteristics (Ulrich et al. 2011). However, few studies have attempted a statistical analysis of polygons in the Lena-Aldan interfluve. Saito et al. (2018) used an UAS to obtain high spatial resolution images (i.e., at a spatial resolution of several centimeters), and to reveal the spatial distribution and size of polygons in anthropogenically disturbed areas in Churapcha (Fig. 2a). These statistical analyses require micro-topographic data, images with resolutions of less than tens of centimeters, and field observations together with a UAS and/or aerial photographss (Ulrich et al. 2011). Although applications of UAS in cryosphere studies have been increasing and very useful (Gaffey and Bhardwaj 2020), the coverage and acquisition timing of UAS images are limited. It is thus necessary to examine whether it is possible to identify and evaluate polygon development using high-resolution satellite images, such as those collected by the WorldView-series optical satellites. It is considered that the extraction and statistical analysis of ice wedge polygon characteristics from optical satellite images in different areas will facilitate a deeper understanding of thermokarst dynamics, especially in the initial stages of thermokarst development.

The objective of this study was to characterize thermokarst subsidence in the study area, and to evaluate the extent of permafrost degradation in the Lena-Aldan interfluve. L-band InSAR data obtained by the Advanced Land Observing Satellite-2 (ALOS-2) platform were used to derive the spatiotemporal changes in surface displacement in the study area. WorldView and Pleiades high-resolution optical images were used to identify ice wedge polygons in areas where subsidence has been detected, and to collect data for statistical analysis. A combination of InSAR-based subsidence measurements and estimates of the degree of polygonal texture development were used to produce an index for evaluating the initial stage of thermokarst development in the disturbed areas.

2 Study area

The Yedoma Ice Complex is widely distributed in the Lena-Aldan interfluve, which is terraced and mainly divided into five terraces named after geomorphology (Soloviev 1959). Figure 2 shows map of the terraces and their elevation. Typical well-developed alas depressions and many thermokarst features are identified mainly on the Tyungyulyu and Abalakh terraces, and both terraces are characterized by ice-rich deposits over 25 m-think in some places (Soloviev 1959; Ulrich et al. 2017). On the Tyungyulyu terrace from Borogontsy to the Mayya area (Fig. 2a and b), many deforested areas (farming and arable land) and alases surrounded by forest are found in the vicinity of towns. The town of Mayya is representative of a residential area located in a typical thermokarst landscape. Yukechi, which is 10 km east of Mayya, is on the Abalakh terrace (Fig. 2a and b) and is the site of field observation station managed by the Melnikov Permafrost Institute in Yakutsk. Thermokarst features were also observed in this area (Fedorov and Konstantinov 2008; Fedorov et al. 2014a; Ulrich et al. 2017). Churapcha is located approximately 200 km to the east of Yakutsk and is located on the same Abalakh terrace as Yukechi. Many dry grasslands and abandoned arable lands, which may be at risk of thermokarst development, are widely distributed around the settlement of Churapcha (Iijima et al. 2021). The population of Churapcha has increased in the last 10 years, which has resulted in expansion of the settlement into the dry grasslands and abandoned arable lands. Therefore, this town was selected for evaluating the effect of thermokarst development in anthropogenically disturbed areas.

The Yedoma ice complex is also distributed along watersheds and on some terraces in river valleys (Soloviev 1959; Iijima and Fedorov 2019), which may be also at a risk of thermokarst development. Amga is located on the left bank of the Amga River, which is one of the tributaries of the Aldan River (Figs. 1 and 2c). This area is located within the Amga-Aldan region of gently sloping central taiga (Fig. 2c), and there are large ice-wedges up to 20 m thick on the left bank (Lytkin et al. 2021). Lytkin et al. (2021) recently conducted field research and geographical analysis to map the extensive thermokarst development in the Amga town. Their findings showed that areas of byllar are distributed to the northwest and northeast of the settlement. However, the spatiotemporal variation in the surface subsidence is unknown. Since numerous arable areas exist on the left bank of the Amga River, we included this site in our study area.

3 Data and methods

3.1 ALOS-2 data and InSAR analysis

This study used ALOS-2 L-band HH-polarized SAR data with a resolution of 10 m (Stripmap mode 3 data) obtained from 2014 to 2021 (Table 1) to examine inter-annual changes in surface displacement in the Lena-Aldan interfluve, Central Yakutia. The advantage of L-band SAR is that it is less affected by surface vegetation and a small temporal decorrelation (Wang et al. 2017). Thus, L-band SAR data were used to estimate inter-annual displacement in areas surrounded by forested taiga in the interfluve. Our analysis method was essentially the same as that used in previous studies (Strozzi et al. 2018; Abe et al. 2020, 2022; Iijima et al. 2021). We used Gamma software to generate Single Look Complex (SLC) images from level 1.1 ALOS-2 data (Wegmüller and Werner 1997). After co-registration of two SLC images, we performed interferometry with subtraction of topographic effects using the ALOS World Digital Surface Model − 30 m- (Takaku et al. 2020). After adaptive spectral filtering (Goldstein and Werner 1998), a phase unwrapping procedure was employed following the method of Costantini (1998). The spatial resolution after the terrain-corrected data were geocoded and projected onto a UTM map projection was 30 m. Finally, the line of sight (LOS, distance between satellite and ground) change detected by InSAR was converted to a vertical displacement by dividing the LOS by the cosine of the incidence angle because vertical displacement is assumed to be dominant in thermokarst. We generated as many interferograms within a period of less than four years as possible.

Table 1 Detail of ALOS-2 data used in this study

Tropospheric and ionospheric disturbances can sometimes interfere with extracting the extent of surface displacement in interferograms. To address the tropospheric effect, the GACOS product was used to reduce the noise (Yu et al. 2017; 2018). The impacts of ionospheric disturbances typically manifest as a long-wavelength trend. Indeed, several unexpected signals were sometimes observed, particularly in forested areas, which may affect the fitting of the long-wavelength trend. To address this issue, we first clipped the targeted areas in each interferogram to exclude forested areas where possible. We then modeled the trend by fitting a 2D polynomial function and subtracted from original interferogram. A manual quality check of the unwrapped interferograms was also performed to remove any pixels affected by major unwrapping errors.

In InSAR analysis, selecting a phase reference point in interferograms has been demonstrated to be important for deriving surface displacement. However, identifying a stable point that is unlikely to be displaced in permafrost regions is difficult. In a study conducted to the southeast of Churapcha, Iijima et al. (2021) used an alas that was not located in a wetland landscape as a reference point, mainly because alases are representative of the climax geomorphological stage of thermokarst development (Bosikov 1991). In this study, averaged InSAR phases at three reference points were used to correct the phase in each interferogram. The reference points used in this study were also in the bottom of alases and in settlements outside wetlands that were free of any polygonal development. This manipulation was performed to reduce localized tropospheric effects, which sometimes appear and could not be reduced by the GACOS due to their spatial scale. The phase differences between each reference point and the averaged one in each region was examined to evaluate that they did not displace excessively relative to the averaged reference. The averaged deviation of phase differences in areas of Mayya, Churapcha, Tyungyulyu, and Amga was 0.17, 0.11, 0.18 and 0.22 cm, respectively.

After generating and correcting interferograms following these procedures, we performed InSAR time-series analysis as described in previous studies (Schmitt and Bürgmann 2003; Biggs et al. 2007, Liu et al. 2015; Yanagiya and Furuya 2020; Abe et al. 2022). Selected InSAR pairs in each site are shown in Figure S1 in Additional file 1. In the time-series analysis, pixels affected by noise were eliminated using a coherence filter (i.e., coherence above 0.35 in 50% of the interferograms, Rouyet et al. 2019). To estimate the propagation errors in the time-series analysis, we assumed that each interferogram contains an error of 0.2 cm (Yanagiya and Furuya 2020; Abe et al. 2022) and created an InSAR data covariance matrix following the method of Biggs et al. (2007). The 2σ propagation error reached up to approximately 2 cm, which depends on the number of epochs in the time-series analysis.

3.2 High-resolution satellite optical images and extraction of polygon statistics

WorldView and Pleiades images (panchromatic: 0.5 m resolution, multispectral: 2 m resolution) were used to identify polygons caused by thermokarst. Details of the datasets used in this study are shown in Table 2. We first performed visual comparisons between an UAS orthorectified image taken in 2017 by Saito et al. (2018) and a WorldView-2 pansharpened image taken in 2021 to confirm whether number of polygons were identical in both image types (Fig. 3a and 3b). Next, we counted the number of polygons in a defined area (white inset) and calculated the density of polygons (Fig. 3b). Note that the selected area should contain only polygons and exclude features such as forests, lakes, and buildings, as much as possible.

Table 2 Detail of high-resolution satellite optical images used in this study

Fig. 3

Visual comparison of a an orthorectified UAS image, and b a WorldView pansharpened image. Each red symbol indicates an identified polygon, and the white rectangle indicates the area in which the density of polygons was calculated. c, d, and e Procedure of spatial frequency analysis. An example of c a clipped panchromatic image, d an amplitude spectrum in 2D spatial domain, and e a 1D amplitude spectrum

In addition to calculating the number of polygons, we attempted to estimate the degree of polygon development by spatial frequency analysis. Polygons were generally circular in shape with a diameter of approximately 10 m, and the borders are represented by the shading in the image (Fig. 3). When the polygons are well-developed, their borders can be easily identified, because we consider the amplitude spectrum in spatial frequencies corresponding to the diameter of polygons becomes large. The spatial frequency analysis was conducted as follows: 1) clipped panchromatic images containing only polygons in an area measuring 100 × 100 m² were prepared (Fig. 3c); 2) discrete Fourier transformation was applied to the clipped panchromatic images to obtain the relationship between spatial frequency and amplitude spectrum in the two-dimensional (2D) spatial domain (Fig. 3d); and 3) the 2D amplitude spectrum was converted to a one-dimensional (1D) amplitude spectrum for fitting to an exponential function (except for spatial frequencies below 0.05 m⁻¹) to obtain a model representing the degree of polygonal texture development (Fig. 3e). The average of the models between well-developed and less-developed polygons was compared to discuss the degree of polygon development.

4 Result

4.1 Ground displacement derived from ALOS-2 InSAR time-series analysis

Figure 4a shows the WorldView-2 image of Mayya town. The settlement is in the top center of the image. Mayya is situated on almost flat ground at an elevation of approximately 140 m above sea level (a. s. l.); the southern part of the town is higher and has an elevation of up to 200 m a. s. l. (Fig. 4b). Alas depressions can be identified by their slightly lower elevation compared to the surrounding inter-alas meadows (Fig. 4b). Figure 4c shows the spatial distribution of the cumulative vertical displacement around Mayya town, with positive and negative signals indicating uplift and subsidence, respectively. Significant subsidence of up to 10 cm was observed in the deforested inter-alas areas. Temporal changes in the displacement at the four selected subsidence sites (M1–M4, Fig. 4d) were almost linear, except near the end of the observations when they appear to have accelerated somewhat.

Fig. 4

(a) WorldView-2 image captured in 2019 over Mayya; (b) Elevation around Mayya; (c) Spatial distribution of cumulative surface displacement. M1–M4 denote sites used for extracting time-series displacement data shown in (d). Stars indicate the reference points used to calibrate InSAR phases; (d) Time-series plots of the four selected sites shown in (c)

Figure 5a shows WorldView-2 images of Churapcha town. The settlement of Churapcha is located to the west and the center of the large lakes in the center of the satellite image. In the south, dry grasslands and localized areas of abundant arable land were distributed. The elevation of Churapcha ranges from 170 to 220 m a. s. l. (Fig. 5b), with areas to the north, west, and southeast having a slightly higher elevation, while areas in the center are lower than the surrounding areas. Figure 5c shows the spatial distribution of the cumulative surface subsidence, and that maximum subsidence of up to 16 cm was observed in an area of abundant arable land in the north. Other areas of arable land and grassland subsided by more than 10 cm. Interestingly, temporal changes in the four selected subsidence sites (C1–C4, Fig. 5d) showed little inter-annual variations and the rates of subsidence were almost linear.

Fig. 5

(a) WorldView-2 image captured in 2011 over Churapcha; (b) Elevation around Churapcha; (c) Spatial distribution of cumulative surface displacement. C1–C4 denote sites used for extracting time-series displacement data shown in (d). Stars indicate the reference points used to calibrate InSAR phases; (d) Time-series plots of the four selected sites shown in (c)

Figure 6a shows the spatial distribution of cumulative surface displacement in the Tyungyulyu area from Borogontsy (north) to Tyungyulyu (south) (see also Fig. 2a), showing that areas where large subsidence was observed (representative areas are indicated by white rectangles) close to local town settlements. Time-series plots of subsidence at three of the selected sites (T1–T3) are shown in Fig. 6b; the plots show that while there were some inter-annual variations, the overall trends between 2014 and 2021 were generally linear. The magnitude of subsidence was up to 8 cm. Figure 6c–e shows enlarged views of the subsidence map for the three areas with satellite optical images (Fig. 6f–h) for comparison. In Syrdakh, areas of subsidence are distributed to the south of the settlement, as well as to the west and east (Fig. 6c and f). The magnitude of subsidence in the area to the south reached up to 8 cm (T1). In Oltyokhsky, areas of subsidence are distributed to the northeast of the settlement (Fig. 6d and g), which corresponds to the distribution of arable land (Fig. 6g). The magnitude of subsidence in these areas was up to 6 cm (T2). In Balyktakh, areas of subsidence are distributed to the northeast of the settlement (Fig. 6e and h), and the magnitude of subsidence was up to 7 cm (T3).

Fig. 6

a Spatial distribution of inter-annual vertical displacement in the Tyungyulyu area from Borogontsy to Tyungyulyu. T1–T3 denote sites selected for extracting the time-series displacement data shown in b. Stars indicate the reference points used to calibrate InSAR phases; b Time-series plots for the three selected sites shown in a; c, d, and e Enlarged view of inter-annual vertical displacement in c Syrdakh, d Oltyokhsky and e Balyktakh, respectively; f, g, and h Pleiades-1 images captured f in 2013 over Syrdakh, g in 2013 over Oltyokhsky, and h in 2019 over Balyktakh, respectively

Figure 7a shows the spatial distribution of cumulative surface displacement in the Amga area from Amga (northeast) to Bolugur (southwest) along the Amga River. Four areas (A1–A4 in white rectangles) showing significant surface subsidence were examined in detail and time-series plots for these areas are shown in Fig. 7b. The changes in subsidence rates were almost linear between 2015 and 2019, followed by a slight increase in the last two years. Enlarged views of the subsidence maps in these areas (Fig. 7c–e) are compared to optical satellite images (Fig. 7f–h). The areas of subsidence in Amga (Fig. 7c) are distributed to the northwest and northeast of the settlement (Fig. 7f), and subsidence of up to 10 cm was detected. In Pokrovka, a large area of subsidence (~ 0.6 km²) was observed to the north of the settlement, and subsidence of 16 cm was detected. Around Bolugur, clear evidence of subsidence was observed around the settlement (Fig. 7a) and significant subsidence was detected in two areas of arable land located to the southwest of the settlement (Fig. 7e and h) in the upstream reaches of the Amga River (Figs. 2c and 7a). The magnitude of subsidence was over 25 cm, which was largest observed in this study.

Fig. 7

a Spatial distribution of inter-annual vertical displacement in the Amga area along the Amga River. A1–A4 denote sites selected for extracting time-series displacement data shown in b. Stars indicate the reference points used to calibrate InSAR phases; b Time-series plots of the four selected sites shown in a; c, d, and e Enlarged view of inter-annual vertical displacement in c Amga, d Pokrovka and e the western part of Bolugur, respectively; f WorldView-2 image captured in 2021 over Amga; g Pleiades-1 image captured in 2013 over Pokrovka; h WorldView-2 image captured in 2017 over the western part of Bolugur, respectively

4.2 Density of polygons and comparison with subsidence and elevation

The top panel of Fig. 8 shows the locations where the polygon density was counted in the four towns; the box plots shown in the bottom panel show the number of polygons in each town. The median polygon density in Mayya was 292 ha⁻¹, which was considerably higher than that in other towns (165 ha⁻¹ in Syrdakh, 76 ha⁻¹ in Churapcha, and 63 ha⁻¹ in Amga). These results imply that the diameter of polygons in Mayya is considerably smaller than those at the other three sites. Examples of polygonal terrain were shown in Figure S2 in Additional file 1.

Fig. 8

(Top panel) Location of sampling points in the selected towns (Syrdakh, Churapcha, Mayya and Amga) for calculating density of ice wedge polygons. White-circle dots indicate the location of displacement time series shown in Figs. 4, 5, 6 and 7, and the rectangle dots indicate the others; (Bottom panel) Box plots of polygon density in the four towns. Median value of polygon density (M) and number of samples (N) are also shown in the figure

Figure 9 shows comparisons between subsidence and elevation, subsidence and polygon density, and elevation and polygon density in areas with well-developed polygons shown in Fig. 8. Results of linear regressions with all data (dashed line) and all data excluding Amga (solid line) were also presented. The R² values for the InSAR-based subsidence and elevation were 0.32 and 0.31 for all data and all-except-Amga data, respectively (Fig. 9a). The R² values for the subsidence and polygon density for all data and all-except-Amga data were much small, 0.07 and 0.17, respectively (Fig. 9b). The R² value for elevation and polygon density for all data was 0.18, but that for all-except-Amga data was 0.45 (p < 0.001, Fig. 9c), which was the highest correlation of all combinations using the three factors (subsidence, elevation, and polygon density).

Fig. 9

Relationship between a subsidence and elevation, b subsidence and polygon density, and c elevation and polygon density. Colored symbols indicate the different towns. Dashed and solid lines indicate linear regressions using all of the markers and all of the markers excluding those from Amga, respectively

4.3 Normalized amplitude spectra in area with polygons

Figure 10a and b shows examples of pan-sharpened images in the well-developed and less-developed polygons in C4 and A1, respectively. Polygons representing circular convexities were clearly identified in C4, but not in A1. We selected each seven areas as well-developed and less-developed polygons like C4 and A1 (see also Additional file 1: Fig. S3) and used the spatial frequency analysis to derive an averaged model that could be fitted to a 1-D amplitude spectrum (Fig. 3e, see also Additional file 1: Fig. S4). Figure 10c shows a comparison of averaged models in 1-D normalized amplitude spectra between clearly- and unclearly- identified polygons. The difference between the two averaged models was pronounced at spatial frequencies in the range 0.05 to 0.2 m⁻¹, which corresponds to a spatial wavelength of 5–20 m. As a comparison with other land covers, two other amplitude spectrum models of water surface and forested areas was also presented (Fig. 10c). The amplitude spectrum at forested areas was much stronger at any spatial frequencies than those polygonal areas, and the frequency dependence at water surface was quite different from those polygonal areas.

Fig. 10

Examples of WorldView pansharpened images indicating a well-developed and b less-developed polygons; c Averaged 1-D spectrum models of (red) well-developed, (blue) less-developed polygons, (yellow green) forested areas, and (cyan) water surface with the standard deviation (shaded areas)

5 Discussion

5.1 Impacts of land use changes on thermokarst development

Inter-annual subsidence maps of the four areas (Figs. 4, 5, 6 and 7) showed that significant subsidence occurred in farming areas and/or in areas of arable land, compared with other land use areas. These open areas were originally covered by forest taiga before the collapse of the Soviet Union in the 1990s; a period of large-scale changes in land use (Fedorov et al. 2014b). The combination of climate change and increased precipitation in the 2000s resulted in warming of permafrost, deepening of the active layer, and melting of ground ice (Iijima et al. 2014), all of which markedly increased the progress of thermokarst development in deforested areas (Iijima et al. 2021). Changes in land use history and subsidence have been investigated previously in Mayya and Churapcha (Abe et al. 2020; Iijima et al. 2021), and the other studied areas (Tyungyulyu and Amga) are also in the vicinity of towns, implying that they are almost in the same situations. These drastic changes have also impacted cattle and horse farming, as well as residential areas that are highly dependent on the local ecosystem (Crate et al. 2017; Takakura et al. 2021). On the contrary to the deforested areas mentioned above, ground in forest areas is considered to be still stable (Fedorov et al. 2019). However, further forest cutting, and climate changes could increase the area affected by thermokarst more than ever before. Therefore, there need to live in a way that does not disturb the ground surface to protect permafrost.

5.2 Relationship among subsidence, elevation and polygon density and its implications

Relationship between InSAR-based subsidence and elevation in areas with well-developed polygons showed that the subsidence seems to be greater at higher elevations (Fig. 9a). However, only data in Churapcha are located at higher elevation than others and they are scattered ranging in a few to several tens of centimeters. Therefore, we consider little correlation between subsidence and elevation. Relationship between subsidence and polygon density (Fig. 9b) showed that the polygon density in Churapcha were scattered regardless of magnitude of subsidence and that the subsidence in other towns were also scattered regardless of polygon density. These suggest that subsidence occurred regardless of polygon density. However, Saito et al. (2018) derived averaged subsidence rates of two sites with different polygon densities in Churapcha and mentioned that the difference of subsidence might be explained by that of volumetric ice contents. We speculate that the difference of relationship between subsidence and polygon density might be caused by that of spatial resolution and analysis method between satellite and the UAS images. The UAS-based “cumulative” subsidence were calculated on the assumption of straight original surface before the thermokarst initiation. The UAS images have high spatial resolutions with a few centimeters, which could detect the bottom of narrow troughs between polygons in both sites. However, the InSAR-based subsidence is a magnitude of that “over a given period of time” and the spatial resolution is 30 m, which might not detect the subsidence distribution in troughs. In initial stage of thermokarst, ground surface gradually subsides particularly in troughs (Bosikov 1991; Desyatkin et al. 2009) and the spatial heterogeneity of subsidence distribution and the difference stages of thermokarst might make the relationship between observed subsidence and polygon density more complicated. The subsidence rate itself is also time-variant due to the stages and meteorological factors, which may be associated with the relationship (Fedorov et al. 2014a).

There was no strong correlation between elevation and polygon density, but with the exception of the Amga data, the correlation was the strongest (Fig. 9c). The elevation and polygon density are time-invariant at least within our observation period, different from subsidence, and they may reflect ground ice structure in the interfluve. Based on those results, we speculate that the higher elevation has older topographic landform and would remain older ice wedges. The studied area is terraced by the Lena and Aldan River, with newer terrace faces closer to the rivers. Therefore, it is likely that newer and smaller ice wedges formed on the Tyungyulyu terrace, which is lower and closer to the rivers, compared with on higher terraces. On the other hand, ice wedges in Churapcha are likely to be older, and have more grown up as times goes by. In view of polygon density, polygons with an average diameter of more than 10 m are distributed in Churapcha. Note that some polygonal patterns with smaller diameters are also distributed within Churapcha, which is thought to be in localized wetter areas (Saito et al. 2018). This situation may suggest that the lower terraces (i.e., Tyungyulyu terrace), the much water it converges, which may form newer and smaller ice wedges. When such ice wedges melt due to recent temperature increases, smaller polygons will be formed on the lower terraces, as observed in Mayya. Thermokarst development is not restricted to sediments with massive ice-wedges; it can occur in areas with smaller ice structures, such as ice lenses in alluvial deposits, and in biogenic deposits (Shestakova et al. 2021). Thus, these differences of polygon density might be affected by the period and soil structure at the time ice wedges developed. Satellite data can reveal ground surface conditions, not in ground, and cryostratigraphic observations such as sampling and analyzing ice wedges are desirable as a future research to better understand differences of polygon sizes and timing of their formation as well as soil conditions at that time (e.g., Opel et al. 2019).

Depositional environment of ice wedges in Amga differs from that at other locations, and that may reflect as the outliers of data in Fig. 9c. This is because Amga is a watershed area that is in close proximity to the Amga River. The subsidence observed in Bolugur along the Amga River was the greatest of all areas studied, and the density of polygons in Amga was small for its elevation (Fig. 9c). Our observed large subsidence may suggest that Amga was underlain by largely developed ice wedges and those melting induced entire subsidence before developing polygonal surface. The presence of high ground ice content along the Amga River was reported by Shestakova et al. (2021) and the extremely wet area along the Amga River is likely to enhance larger ice wedges due to its high ground ice content and/or ground structure. In addition, Amga has experienced occasional spring floods (Tananaev et al. 2021), which may provide excess water supply to the frozen ground and to have grown ground ice and ice lens up possibly that can cause significant subsidence.

5.3 Evaluation of polygonal pattern development based on spatial frequency analysis

There were some areas where the polygonal patterns were unclear, but where clear subsidence was detected (e.g., Bolugur in Fig. 7e). These findings indicate that such polygons are in developmental stage (i.e., byllar) and that the troughs between polygons may be shallow. These conditions prevented us from counting number of polygons everywhere using satellite optical images. Micro-topographic data and super-high-resolution images derived from UAVs or airplanes may facilitate a more detailed statistical analysis of these polygons (Ulrich et al. 2011; Saito et al. 2018; Lytkin et al. 2021). In recent years, analytical techniques have been also improved; for example, a study employed machine learning to identify polygons from the terrain (e.g., Zhang et al. 2020). However, since obtaining training and/or validation data for such approaches can be difficult in our study area, it would be good to have some indices of polygon development derived solely from satellite images.

Our study performed spatial frequency analysis using WorldView and Pleiades panchromatic images to reveal differences in the mean values obtained from models at spatial frequencies and how these correspond to polygon development. We consider that the differences of amplitude spectra between well-developed and less-developed areas are relevant to those of trough depths, which may be also related to the density of vegetation. The water content tends to be high, especially in the regions between troughs, and vegetation may grow as troughs get deeper (C4 in Fig. 10a), which may affect the amplitude spectrum. The spectrum at forested areas was stronger than that at deforested areas regardless of spatial frequencies (Fig. 10c). This means that forest has structures of all size, from small to large such as leaves, branches, trunks and trees, all of which reflect strongly. Conversely. the frequency dependence at water surface was quite different from those polygonal areas (Fig. 10c), where the spectrum was weak regardless of frequency. These results may suggest the possibility of evaluating polygon development could be performed using satellite optical images. However, our study only presented a qualitative discussion of the degree of polygon development. For a quantitative discussion, a comparison with microtopographic data captured by UAVs is needed; such an approach will be the subject of future works. The physics of polygon development and model-based studies (e.g., Abolt et al. 2020) will also help to develop the evaluation method.

5.4 Effectiveness of InSAR-based subsidence and polygonal texture degree estimation for evaluating thermokarst development

InSAR can detect ground subsidence with a spatial resolution of tens of meters, which is useful for investigations focusing on regional-scale thermokarst subsidence. Our study succeeded in clarifying the spatial extent of subsidence in areas of arable land and farmland at a scale of several-hundred meters. On the other hand, the subsidence associated with thermokarst development can be strongly heterogeneous at the scale of several meters, which can make it difficult to compare InSAR-based data with field measurement data. Abe et al. (2020) compared InSAR-based subsidence data with field-leveling data in Mayya in the Lena-Aldan interfluve. They found that localized thermokarst development covered in a-few-tens meters, and they concluded that the discrepancy between the two techniques was due to differences in spatial resolution between InSAR and leveling. To estimate the distribution of thermokarst subsidence at a higher spatial resolution, a combination of InSAR data and UAS data is considered preferable. Iijima et al. (2021) compared InSAR-based subsidence with that obtained using UAS-based digital surface models (DSMs) and found that the two methods were roughly in agreement. However, due to the low vertical accuracy of DSMs and the UAS data spanned only 1 year, producing an accurate subsidence map from UAS data was not possible. We consider that a more detailed study of the subsidence distribution could be produced using InSAR and UAS-DSM data that have been collected over a period of more than 5 years, as the expected subsidence will exceed the level of uncertainty associated with UAS-DSM data (− 10 cm).

To better understand the role of freeze–thaw cycles in ground displacement, simultaneous extraction of seasonal and inter-annual changes is necessary. In this study, we could not estimate seasonal displacement directly due to the limited number of observation data by ALOS-2 over the study area. Future L-band SAR missions, such as the ALOS-4 by JAXA and NISAR by NASA and ISRO, are expected to markedly increase the frequency of observations over the Lena-Aldan interfluve, which will help to reveal temporal changes in the ground displacement more detail. In addition, geophysical models for explaining InSAR-based displacements in permafrost regions are also required in order to evaluate ground ice conditions (Abe et al. 2022). If we can use these models to estimate the amount of ground ice melting, then we will be able to evaluate permafrost degradation more quantitatively and improve our capability for future predictions and adaptation.

Our study shows that the use of satellite optical images to extract statistical characteristics of polygons and estimate their degree of development has the potential to optical satellite images is useful for evaluating thermokarst conditions in the interfluve. Increasing the number of optical satellite images in conjunction with validation data will surely advance the visualization and statistical analysis of polygons and better illustrate the extent of degradation in permafrost regions.

6 Conclusions

Thermokarst development has been investigated intensively in the Lena-Aldan interfluve in eastern Siberia, which is underlain by continuous ice-rich permafrost. However, there were no methods for evaluating thermokarst development over entire towns so far. This study used satellite radar images to characterize the extent of subsidence near settlements in the interfluve. Our results clearly showed that subsidence of up to 25 cm in 2014–2021 had occurred in region, but the magnitude of subsidence varied among towns. The analysis of satellite optical images that were used to estimate the densities of ice wedge polygons in those areas revealed significant differences at each area, possibly due to differences in ground structure and geomorphologic characteristics. Spatial frequency analysis of the optical images showed clear differences in averaged spectrum models describing well-developed and less-developed polygons, which may reflect trough depths and the density of vegetation in the troughs. Our findings showed that InSAR-based subsidence maps combined with the optical image analysis of polygon development may be useful for evaluating thermokarst development in permafrost regions.