Coastal Retreat Due to Thermodenudation on the Yugorsky Peninsula, Russia during the Last Decade, Update since 2001–2010

: Thermodenudation on the Kara seacoast, the Yugorsky Peninsula, Russia, is studied by analyzing remote-sensing data. Landforms resulting from the thaw of tabular ground ice, referred to as thermocirques, are formed due to polycyclic retrogressive thaw slumps, during the last decade 2010–2020. We calculate the retreat rate of the thermocirque edge using various statistical approaches. We compared thermocirque outlines by the end of each time interval deﬁned by the dates of available very-high-resolution imagery. Six thermocirques within two key sites on the Yugorsky peninsula are monitored. We correlate each of the thermocirque edge’s retreat rates to various climatic parameters obtained at the Amderma weather station to understand the interrelation patterns better. As a result, we ﬁnd a very low correlation between the retreat rate of each thermocirque and summer warmth, rainfall, and wave action. In general, the activity of thermodenudation decreases in time from the previous decade (2001–2010) to 2010–2020, and from 2010 towards 2020, although the summer warmth trend increases dramatically. A single thermocirque or series of thermocirques expand in response to environmental and geological factors in coastal retreat caused by thermodenudation.

Coastal dynamics results from the activation of thermodenudation, the most widespread process in the area with tabular ground ice enclosed in geological sequences. The importance of the massive (tabular) ground ice occurrence as a major factor in coastal retreat is recognized in the literature [7,[21][22][23][24][25].
We marked thermodenudation as the dominating coastal destruction mechanism in the study area [24]. Taken the process's polycyclic pattern, we use different terms for the landslide landform and landsliding as a process instead of one term retrogressive thaw slump (RTS). Term RTS causes complications when studying various stages and mechanisms of coastal retreat and slope mass waste. In Russian literature, we use the term thermodenudation, meaning tabular ground ice thaw and slope mass waste for the process of RTS formation. We use the term earth flow or RTS for a single feature produced by this process. We use the term "thermocirque" (TC) for an extensive landform resulting from a series of multiple-aged RTS occupying various positions within one landform (Figure 1). In this paper, we will use the term RTS only for the single earth flow path. TC at the coasts are step-shaped with two retreating planes and respective edges: a TC backwall with retreating "upper" edge, equal to headwall+headscarp after [13], and a dropwall to the beach with retreating "lower" edge ( Figure 1b). Monitoring by remote-sensing data analysis allows most reliably locating the upper edge of TC and measuring its retreat rate.
The Yugorsky peninsula coast of the Kara Sea near Amderma settlement was studied since 1998. It was established then that warmer summers resulted in an increased rate of TC growth [2,10,24,[26][27][28][29][30][31][32]. Our particular interest in the study of coastal retreat due to thermodenudation at Yugorsky coast is related to a probable reduction in the retreat in the last decade due to warming that could have caused insulation of ice exposures by excessive amounts of thawed material, specifically for low headwalls as suggested in [33]. To resolve this task, we applied processing of very-high-resolution imagery of various time spans for 2010-2020, covering areas with thermodenudation activities monitored in 2001-2010.

Study Area
The Yugorsky Peninsula at the southern coast of the Kara Sea ( Figure 2) is a Pay-Khoy Ridge piedmont, with rolling hills, continuous permafrost zone, layers of tabular ground ice are the main features [34][35][36]. Two key sites are studied: the Pervaya Peschanaya key site (PP), west of the Pervaya Peschanaya River, and the Shpindler key site (Sh), west of the Hubt'Yakha River. Three TCs are identified within each site: the Western (WTC), the Central (CTC), and the Eastern (ETC). In this paper, we will use the term RTS only for the single earth flow path. TC at the coasts are step-shaped with two retreating planes and respective edges: a TC backwall with retreating "upper" edge, equal to headwall+headscarp after [13], and a dropwall to the beach with retreating "lower" edge ( Figure 1b). Monitoring by remote-sensing data analysis allows most reliably locating the upper edge of TC and measuring its retreat rate.
The Yugorsky peninsula coast of the Kara Sea near Amderma settlement was studied since 1998. It was established then that warmer summers resulted in an increased rate of TC growth [2,10,24,[26][27][28][29][30][31][32]. Our particular interest in the study of coastal retreat due to thermodenudation at Yugorsky coast is related to a probable reduction in the retreat in the last decade due to warming that could have caused insulation of ice exposures by excessive amounts of thawed material, specifically for low headwalls as suggested in [33]. To resolve this task, we applied processing of very-high-resolution imagery of various time spans for 2010-2020, covering areas with thermodenudation activities monitored in 2001-2010.

Study Area
The Yugorsky Peninsula at the southern coast of the Kara Sea ( Figure 2) is a Pay-Khoy Ridge piedmont, with rolling hills, continuous permafrost zone, layers of tabular ground ice are the main features [34][35][36]. Two key sites are studied: the Pervaya Peschanaya key site (PP), west of the Pervaya Peschanaya River, and the Shpindler key site (Sh), west of the Hubt'Yakha River. Three TCs are identified within each site: the Western (WTC), the Central (CTC), and the Eastern (ETC).

Pervaya Peschanaya Key Site
PP occupies ancient thermodenudation depressions sloping up from the seashore relatively steeply from 20 to 45 m elevation towards the Pay-Khoy Ridge. Active TCs expose Quaternary sandy and silty deposits, enclosing one layer of tabular ground ice [27,30,35,36]. Layers of tabular ground ice exposed in the backwall of each TC are 2-6-mthick. The edge of the backwall is 50 to 200 m away from the beach, and thus is not directly affected by the wave action ( Figure 3). Before 2010, only two TCs were mapped and monitored at PP. Activation of WTC and ETC occurred in 2000-2001. From 2010-2020, three TCs are observed ( Figure 3). The latest field surveys discovered that the upper limit of the tabular ground ice layer within the TC backwall actively retreating is relatively close to the surface [37,38]. Average retreat rate measured by land-based surveys along transects in 2001-2007 range at 1.6-4.2 m/year with the maximum around 5.8-6.3 m/year.

Pervaya Peschanaya Key Site
PP occupies ancient thermodenudation depressions sloping up from the seashore relatively steeply from 20 to 45 m elevation towards the Pay-Khoy Ridge. Active TCs expose Quaternary sandy and silty deposits, enclosing one layer of tabular ground ice [27,30,35,36]. Layers of tabular ground ice exposed in the backwall of each TC are 2-6-m-thick. The edge of the backwall is 50 to 200 m away from the beach, and thus is not directly affected by the wave action ( Figure 3).

Pervaya Peschanaya Key Site
PP occupies ancient thermodenudation depressions sloping up from the seashore relatively steeply from 20 to 45 m elevation towards the Pay-Khoy Ridge. Active TCs expose Quaternary sandy and silty deposits, enclosing one layer of tabular ground ice [27,30,35,36]. Layers of tabular ground ice exposed in the backwall of each TC are 2-6-mthick. The edge of the backwall is 50 to 200 m away from the beach, and thus is not directly affected by the wave action ( Figure 3). Before 2010, only two TCs were mapped and monitored at PP. Activation of WTC and ETC occurred in 2000-2001. From 2010-2020, three TCs are observed ( Figure 3). The latest field surveys discovered that the upper limit of the tabular ground ice layer within the TC backwall actively retreating is relatively close to the surface [37,38].

Shpindler Key Site
Sh is only 30 km off the PP (Figure 2) yet has a very different appearance. This site occupies a 30-45-m-high terrace composed of sandy and clay deposits enclosing two Before 2010, only two TCs were mapped and monitored at PP. Activation of WTC and ETC occurred in 2000-2001. From 2010-2020, three TCs are observed ( Figure 3). The latest field surveys discovered that the upper limit of the tabular ground ice layer within the TC backwall actively retreating is relatively close to the surface [37,38].

Shpindler Key Site
Sh is only 30 km off the PP (Figure 2) yet has a very different appearance. This site occupies a 30-45-m-high terrace composed of sandy and clay deposits enclosing two tabular ground ice layers. The upper layer is up to 12 m thick, and the lower layer is up to 5 m thick [27,34,35,[39][40][41].
Three TCs at various states of activity were mapped here ( Figure 4). Only ETC and CTC were monitored in 2001-2010, as WTC was stabilized then. After 2010, all the three TCs showed evidence of activity.
Remote Sens. 2021, 13, x FOR PEER REVIEW 4 of 21 tabular ground ice layers. The upper layer is up to 12 m thick, and the lower layer is up to 5 m thick [27,34,35,[39][40][41]. Three TCs at various states of activity were mapped here ( Figure 4). Only ETC and CTC were monitored in 2001-2010, as WTC was stabilized then. After 2010, all the three TCs showed evidence of activity.

Materials and Methods
The study is based on remote-sensing monitoring of TC edge retreat and the relation of retreat rate to several climatic parameters as direct and indirect controls of thermodenudation.

Processing Remote-Sensing Data
We used a series of multi-temporal satellite images of the high and very high spatial resolution of 0.5 to 2.0 m/pixel, dated from 2010 to 2020 (Table 1). To improve the accuracy of orthorectification and imagery spatial alignment in areas with significant changes in surface height (areas of TC growth), we propose a methodology of DEM reconstruction.

Materials and Methods
The study is based on remote-sensing monitoring of TC edge retreat and the relation of retreat rate to several climatic parameters as direct and indirect controls of thermodenudation.

Processing Remote-Sensing Data
We used a series of multi-temporal satellite images of the high and very high spatial resolution of 0.5 to 2.0 m/pixel, dated from 2010 to 2020 (Table 1). To improve the accuracy of orthorectification and imagery spatial alignment in areas with significant changes in surface height (areas of TC growth), we propose a methodology of DEM reconstruction. The spatial alignment accuracy of multi-temporal images is important, as it determines the reliability of retrieved data on TC edge outline and retreat dynamics. High accuracy of imagery spatial alignment, as a rule, is provided by orthorectification with ground control points (GCPs) and a digital elevation model (DEM), as was carried out in [15]. However, our study area is not provided with instrumentally measured GCPs, so we used block adjustment for the spatial alignment of multi-temporal images.
The most optimal is the processing of stereopairs of satellite images for each analyzed date since it allows creating DEM consistent with satellite images, and, consequently, performing accurate orthorectification. However, due to the lack of historical images sufficient for this type of processing, we used the method of reconstructing the available DEMs in the zone of active TC or coastal retreat.
This methodological approach is relevant when images precede the available detailed DEM date. The DEM is "restored", that is, the TC areas are reconstructed in the direction from the edge to the center. Thus, we obtained a retrospective DEM with an initial (assumed) height of all analyzed lines of the TC edges. Such an approach with DEM restoration is advisable for the assessment of the coastal dynamic. This method was successfully tested in our study of the dynamics of the Kolguev Island coasts [42].
The current method is based on morphological dilation, applied to images with discrete brightness values (in our case, altitudes of the DEM matrix). The essence of the method is to assign the maximum value obtained within a structural element belonging to a given neighborhood (local window) to the current pixel of the processed matrix. Mathematically, this can be written as follows: where i, j-the numbers of the row and column of the matrix, g-the reconstructed matrix, f -the structural element, and W-the local window. The size and shape of the structural element both were selected experimentally followed by expert judgment of the results obtained. Hence, we found that the most suitable configuration for the current problem is a local window by 3 × 3 pixels, and the shape is as follows:  Convolution was applied iteratively on a manually defined area of interest, while several iterations controlled the reconstruction width.
We used ArcticDEM Strips [43] as a base for reconstruction. The ArcticDEM Strip, compiled from stereo-pair WorldView-1 of 4 September 2011, was used as a base for PP, and ArcticDEM Strips, compiled from stereo-pairs WorldView-1 and WorldView-3 of 14 June 2014 and 16 September 2015, respectively, were used as a base for Sh. Before performing a photogrammetric task, ArcticDEM Strips processing includes: (1) masking water bodies; (2) replacing single extreme of altitude by median filtering; and (3) expert judgment, masking, and interpolation of incorrect portions of the DEM. To mask water bodies, water identification was performed, and a vector polygonal mask was created, then, the values of the DEM pixels belonging to the mask were assigned the value of the average height along the mask edge. Then, the filtering and interpolation of gaps were undertaken.
Then, reconstruction of topography for retreated parts of the coast was performed. We applied 15 iterations, which correspond to 30-m-width reconstruction (according to ArcticDEM Strips pixel size).
Photogrammetric processing includes two main stages: (1) pairwise spatial alignment of satellite images; and (2) orthorectification. Block bundle adjustment with automatic tie-point measurements and zero-order correction for RFM (WorldView data), as well as rigorous model (Formosat and SPOT data), were all used at an alignment step. Then, image verification was undertaken to prove that TC edges on all images were within the reconstructed area. Reconstructed DEMs were used at an orthorectification step and aligned multi-temporal satellite images with sub-pixel accuracy. Imagery relative spatial alignment calculated as a root mean square error for each image is shown in Table 1. Processing is performed using software SCANEX Image Processor, version 3.5.

Determining Thermocirque Dynamics (TC Edge Retreat Rates)
Coastal TC edges were digitized manually on orthorectified imagery. The TC dynamics was analyzed only in areas with active growth, where the edges were represented by a sharp line that was confidently interpreted on satellite images. Variations in spatial resolution of images and their superposition accuracy resulted in some errors drawing TC edge lines. This error for a consecutive pair of TC edges was estimated as a square root of squares' sum from Table 1. Root mean square error accuracy allowed estimating resulting spatial alignment error within 0.4 to 1.1 m for a consecutive pair of TC edges. Thus, we obtained an accurate mutual spatial alignment of multi-temporal orthorectified imagery, comparable with achieved in [15]. This allowed a detailed assessment of TC edge dynamics.
The TC edge retreat was measured according to the method applied in [15] along the transects placed every 10 m using the digital shoreline analysis system (DSAS) v5.0 add-in to Esri ArcGIS [44]. DSAS allowed the calculation of the dynamics of changes in the position of multi-temporal coastlines along a series of transects, which were automatically built perpendicular to the baseline. This approach provides the best fit to the perpendicular of transects to both TC edge lines between which they appear. We calculated the baseline position for our TCs as Voronoi polygons between TC edge lines 2010 and 2020. Some transect lines were reoriented and reshaped as necessary to reflect the compound form of TC growth. These transects were used to measure the distances between TC edge lines using net shoreline movement (NSM) calculation. Thus, we obtained a series of TC edges retreat values for each TC.
The type of statistical distribution of the obtained values ( Figures A1 and A2) does not obey the Gaussian distribution. This statement was confirmed by additional statistical tests using D'Agostino's K-squared method. The resulting p-value was much less than alpha (5% or 0.05) (Tables A1 and A2). So, this does not allow using the mean arithmetic when evaluating the changes in the rate of TC edge retreat with time.
Statistical measures used in this work are the median, the maximum, the first (Q 1/4), and the third (Q 3/4) quartiles following the method applied in [20]. Results are presented in box-whisker plots. Median characterizes the highest probability of TC edge linear retreat within each time interval under consideration. The interquartile range between Q1/4 and Q3/4 represents the distribution range covering the most probable retreat values, cutting off the rarest minima and maxima.

Processing Climatic Data for 2010-2020
Data were obtained by processing the records from weather station Amderma (WMO# 23022, 69.77°N, 61.68°E). The interannual variability and main trends of air temperature and precipitation are shown on a diagram ( Figure 5). One can see positive trend of air temperature and decline of atmospheric precipitation during two decades of study. However, we had to operate with values placed within the time interval between two adjacent acquisitions, which differ from the standard year period. So, all records were summarized for the period between acquisitions (Table 1). This time interval varied in duration for two key sites because different sets of remote-sensing data were available for each. To eliminate the difference between time intervals for two key sites, we processed However, we had to operate with values placed within the time interval between two adjacent acquisitions, which differ from the standard year period. So, all records were summarized for the period between acquisitions (Table 1). This time interval varied in duration for two key sites because different sets of remote-sensing data were available for each. To eliminate the difference between time intervals for two key sites, we processed climatic records, as presented in Table 2. The parameter for the entire time interval was subdivided (normalized) by the warm period duration (WPD). Parameters, both summarized for the warm period of each time interval and normalized by WPD, are thaw index, degree-days; rainfall, mm; the repeatability of northern wind (wind with a

Results
We have to operate with the available data characterizing periods of different duration between measurements, with a different set of climatic parameters, though obtained from one weather station (Table 2).

Overview of Retreat Pattern within the Yugorsky Coastline
The retreat rate of TC edges at two key sites changed significantly within subdivided time intervals (Figures 6-8). As seen from diagrams in Figure 8, the retreat rate is higher at PP (Figure 8, left column, bottom pane) compared to Sh (Figure 8, right column, bottom pane). The higher range (Q1/4-Q3/4) corresponds to a higher scattering of data, and higher values of Q1/4 indicate less zero retreat transects within each TC at PP compared to Sh (Figure 8, bottom panes). So, at Sh with generally lower retreat rates, practically no stable portions of the edges are observed. In 2010-2020, medians fluctuated with the maximum in 2010-2012 at PP and 2013-2015 at Sh, minima being asynchronous as well ( Figure 2).         Each separate TC ( Figure 6) pattern differed ETC showed a maximum retreat rate in 2016-2017, when other TCs were at their minimum. If we excludes the influence of different observation intervals, there would be no pronounced peaks, and the retreat activity would decrease almost uniformly over time (Figure 8, left column, middle pane, Table 2).

Climatic Controls of Coastal Retreat
Analysis of diagrams for the area and linear retreat for both key sites, both summarized and in each TC separately (Figure 8, middle panes), were compared to climatic parameters ( Figure 8, upper panes) and showed considerable differences in the rates of TC growth versus climatic parameters.
The median and maximum indicators of retreat in PP are at their peak in 2010-2012. However, this time interval is characterized by a relatively cool summer, moderate northern winds, an open water duration, and only high rainfall. The hot summer of 2012 being added to the cool summer of 2010-2011 probably did not add much to a warm peak, yet was sufficient to provide the noticeable activity of the TC edge retreat though added to less active years within the time interval (Figure 8, bottom pane, Table 2). The hot summer of 2016 is included in the period from 17 June 2016 to 27 July 2017, and affects maximum retreat rate indicator (not median and average values, though).
In Sh, peaks of retreat rate indicators are not so consistent in a pattern. Various indicators shift in time relative to each other and between various TCs. The maximum retreat rate indicator is at its peak from 26 June 2012 to 17 July 2013. The hot summer of 2012 is included here in contrast to PP. The median peak occurs from 17 July 2013 to 17 August 2015, which includes two half-summers and one whole summer. However, a medium thaw index is observed, as well as medium precipitation, and the lowest strong northern wind repeatability plus open-water duration (Figure 8, bottom pane, Table 2). Therefore, we explain the median peak only due to summarizing several years of retreat. The other close-in-value peak of the median occurs from 26 June 2012 to 17 July 2013, with one summer composed of two-year portions, including the hot summer of 2012, practically providing the highest relative median.
Detailed analysis of the retreat rates within time intervals showed the weak relation of retreat indicators to climatic parameters or their combination (Figure 9). In general, the best of low correlations belongs to the maximum retreat rate indicator versus the thaw index for PP (R 2 =0.433). For Sh, the best fit corresponds to strong northern winds along with a high open water duration. The correlation is negative with R 2 =0.21-0.48 for the wind repeatability versus all the retreat-rate indicators, and R 2 =0.51 between the open water period and the median retreat rate indicator. It is hard to suggest an explanation to the effect of open water because all backwalls there are 300-500 m off the shoreline, and dropwalls are 7-12-m-high over the beach.
Consistency in the highest values of retreat parameters observed at the PP site (Table 3) compared to the inconsistent distribution of the highest retreat rates at the Sh site can be explained, in part, by the dates of activation. At PP, activation of TCs was observed in 2000-2001, while, at Sh, active TCs were already present on the aerial images of 1947. We suggest that mature TCs of Sh manifest individual rate of backwall retreat to a higher extent compared to "younger" TCs of PP.
A better correlation is obtained when studying the interrelation of climatic parameters and the number of RTS. It was noted [15] that warm summers of 2011 and 2012 are promoted over 200 RTS. However, the same publication states that there is no correlation between individual RTS retreat and warming as well as precipitation. Ref [13] determined that, during 1984-2015, the number of RTS increased 60-fold due to warm summers. Unlike Canada and Russia, in China [14], activation was observed in 2010 and 2016. This activation was mainly attributed to anomalously high summer temperature and abundant precipitation. Intensive rainfall is a reason for significant activation [5]. The increase in TC number in Central Yamal was related to the warm summer of 2012, but retreat rates did not increase in the even warmer summer of 2016 [12,45], which is the same as in the study area on the Yugorsky peninsula.
Other climatic controls of the coastal retreat are also a matter of analysis in several papers. The highest coastal retreat rates are measured on the north-west-facing shorelines, which correspond to strong waves' main direction [46]. Relation of retreat rate to scarp height was found in [47]. They also established that slope orientation does not affect the retreat rate of RTS. In the mountainous regions of the Canadian Arctic, the slope aspect plays an essential role in RTS development. South-and west-facing slopes show higher retreat rates [33]. Ice thickness was mentioned as an essential factor of RTS distribution, while RTS activity and initiation depend on a slope angle [7]. The impact of ice thickness for PP was proved at the previous monitoring stage [30].
All weaknesses in the relationship between retreat rates and climatic parameters should be attributed to other factors, such as relief, environment, ice content, and distribution in the section. Environmental controls of retreat in the study area were suggested [32] and relatively reliably show the retreating edges of TCs in 2010-2020 expansion into less resistant environmental units.

Comparing Retreat Rates in the Yugorsky Peninsula and Other Regions of the World
The long-term retreat rate calculated from the analysis of remote-sensing data at the Yugorsky coast before 2010 [24] was close to that obtained for various arctic areas, approximately 1 m/year for the period 1948-2001. Somewhat similar areas of the Canadian Arctic show a retreat rate of 1.03 m/year in 1970-2000 [23]. Though records averaged for an extended period (53 years for the Yugorsky Peninsula and 30 years for Canada) show relatively low values, annual data display much higher rates.

Regions Farther North from the Yugorsky Peninsula
In the Eureka Sound Lowlands area, the Ellesmere and Axel Heiberg Islands, the average individual TC retreat rates in 2011-2018 comprised 0-26.7 with a maximum of 79 m/year [15]. The maximum rates were noted in 2011-2012 due to high summer temperature. At Yugorsky, in the same period of 2010-2012, an intensive growth of TCs and the first peak of maximum rates (59.1 m for the period) were recorded at the PP. After the peak of 2011-2012, the growth rates in most RTS decreased by 2018 [15]; in the Yugorsky Peninsula, the trend of decreasing thermodenudation rates is similarly observed.

Conclusions
We applied a method for processing multi-temporal satellite images that ensure the accurate alignment of analyzed images. Thus, we calculated the dynamics of TC edges with sub-meter accuracy even without support by land-based surveys.
Measuring the TC edge retreat for several time intervals from 2010-2020 showed an oscillatory character of TC development. In contradiction with our previous studies and some published observations, retreat rates of TC edges do not depend on any climatic drivers at least for the study period. Maximum retreat shows stronger dependence on the thaw index than average values. High medians coinciding with the peak of maximums for a given time interval reinforce the superposition of favorable climatic conditions on the favorable environmental features of the area. The maximum retreat rates, not accompanied by high medians within a time interval, suggest the superposition of high summer warmth and relatively unfavorable environmental conditions along the retreating edge with only local exceptions producing complicated configuration of the retreating TC edge characteristic of the Yugorsky Peninsula.
Different TCs in two key sites develop at their rate, some activate in much warmer but are drier in 2010-2020 compared to 2001-2010; however, in general, they stabilized within a two-decade time interval. Stabilization is most likely caused by increased accumulation of thawed sediments insulating ice exposures, but may also respond to changes in ice layer thickness reduction at different portions of the coastal plain. Data Availability Statement: Publicly available climatic dataset was analyzed in this study. This data can be found here: http://meteo.ru/, accessed on 20 June 2020, weather station #23022, Amderma. Remote-sensing data was obtained from Airbus Defence, Maxar and Scanex due to license agreements, and are available at https://www.airbus.com/space.html, accessed on 20 June 2020; https://www. maxar.com/, accessed on 20 June 2020; and https://www.scanex.ru/en/, accessed on 20 June 2020, with the permission of Airbus Defence, Maxar, and Scanex, respectively.
Acknowledgments: Centre of collective usage «Geoportal», Lomonosov Moscow State University provided access to remote sensing data.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.