The Role of Thermal Denudation in Erosion of Ice-Rich Permafrost Coasts in an Enclosed Bay (Gulf of Kruzenstern, Western Yamal, Russia)

Baranskaya, Alisa; Novikova, Anna; Shabanova, Natalya; Belova, Nataliya; Maznev, Stepan; Ogorodov, Stanislav; Jones, Benjamin M.

doi:10.3389/feart.2020.566227

ORIGINAL RESEARCH article

Front. Earth Sci., 19 January 2021

Sec. Cryospheric Sciences

Volume 8 - 2020 | https://doi.org/10.3389/feart.2020.566227

This article is part of the Research Topic Observations, Interactions, and Implications of Increasingly Dynamic Permafrost Coastal Systems View all 15 articles

The Role of Thermal Denudation in Erosion of Ice-Rich Permafrost Coasts in an Enclosed Bay (Gulf of Kruzenstern, Western Yamal, Russia)

Alisa Baranskaya^1,2*

Nataliya Belova^1,3

¹Laboratory of Geoecology of the North, Faculty of Geography, Lomonosov Moscow State University, Moscow, Russia
²Laboratory of Environmental Paleoarchives, Institute of Geography, Russian Academy of Sciences, Moscow, Russia
³Earth Cryosphere Institute, Tyumen Scientific Centre, SB RAS, Russia
⁴Institute of Northern Engineering, University of Alaska Fairbanks, Fairbanks, AK, United States

Coastal erosion in the Arctic has numerous internal and external environmental drivers. Internal drivers include sediment composition, permafrost properties and exposure which contribute to its spatial variability, while changing hydrometeorological conditions act as external drivers and determine the temporal evolution of shoreline retreat. To reveal the relative role of these factors, we investigated patterns of coastal dynamics in an enclosed bay in the southwestern Kara Sea, Russia, namely the Gulf of Kruzenstern, which is protected from open-sea waves by the Sharapovy Koshki Islands. Using multitemporal satellite imagery, we calculated decadal-scale retreat rates for erosional segments of the coastal plain from 1964 to 2019. In the field, we studied and described Quaternary sediments and massive ground-ice beds outcropping in the coastal bluffs. Using data from regional hydrometeorological stations and climate reanalysis (ERA), we estimated changes in the air thawing index, sea ice-free period duration, wind-wave energy and total hydrometeorological stress for the Gulf of Kruzenstern, and compared it to Kharasavey and Marre-Sale open-sea segments north and south of the gulf to understand how the hydrometeorological forcing changes in an enclosed bay. The calculated average shoreline retreat rates along the Gulf in 1964–2010 were 0.5 ± 0.2 m yr⁻¹; the highest erosion of up to 1.7 ± 0.2 m yr⁻¹ was typical for segments containing outcrops of massive ground-ice beds and facing to the northwest. These retreat rates, driven by intensive thermal denudation, are comparable to long-term rates measured along open-sea sites known from literature. As a result of recent air temperature and sea ice-free period increases, average erosion rates rose to 0.9 ± 0.7 m yr⁻¹ in 2010–2019, with extremes of up to 2.4 ± 0.7 m yr⁻¹. The increased mean decadal-scale erosion rates were also associated with higher spatial variability in erosion patterns. Analysis of the air thawing index, wave energy potential and their total effect showed that inside the Gulf of Kruzenstern, 85% of coastal erosion is attributable to thermal denudation associated with the air thawing index, if we suppose that at open-sea locations, the input of wave energy and air thawing index is equal. Our findings highlight the importance of permafrost degradation and thermal denudation on increases in ice-rich permafrost bluff erosion in the Arctic.

Introduction

Coasts in permafrost regions are extremely sensitive to environmental changes; coastal erosion is among the major destructive geomorphic processes in the Arctic, with circum-Arctic average shoreline retreat rate reaching 0.5 m yr⁻¹ (Lantuit et al., 2013). While coastal erosion in low latitudes mainly depends on wind and wave energy and sediment composition of the shores (e.g., Luijendijk et al., 2018), as well as sea level rise (Meyssignac and Cazenave, 2012) and human impacts, drivers of the Arctic coastal dynamics are unique with respect to the added influence of permafrost and ground-ice (internal drivers) and sea-ice dynamics (external drivers). Together, they comprehend such highly variable processes as thermal abrasion, or thermal and mechanical destruction of permafrost coasts by waves, and thermal denudation, or thawing of the frozen ground at the bluffs and slumping of material due to gravity (Aré, 1988; Razumov, 2001; Leontiev, 2003).

The temporal evolution of coastal erosion rates in permafrost areas mainly depends on changing hydrometeorologic factors, such as air temperature, wind-wave energy, ice-free period, etc. (Vasiliev et al., 2005; Vasiliev et al., 2006; Shabanova et al., 2018). As the frozen cliffs thaw in contact with warmer air and water, shoreline retreat rates are sensitive to air temperature changes (Jones et al., 2009), increasing along with the modern warming in the Arctic (Overland et al., 2019), which exceeds the average world trends (IPCC, 2014; Savo et al., 2016) and trends for the Northern Hemisphere (Serreze et al., 2008). The thawed sediments are further removed by waves, the energy of which depends on wind speed and wind direction. Sea ice is a considerable limiting factor as well, protecting the coasts in winter, and controlling the length of the wave fetch in the sea ice-free period (Barnhart et al., 2014; Ogorodov et al., 2016). With the sea ice reduction in the Arctic since the early 2000s (Kay et al., 2011; Overeem et al., 2011; Stroeve et al., 2014; Barber et al., 2017), the sea ice-free period has become longer (Barnhart et al., 2016), and the length of the wave fetch has increased. As a result, stronger waves affect the coast for a longer period of time each year (Ogorodov et al., 2013; Nielsen et al., 2020), while at the same time there has been an increase in storm frequency (Atkinson, 2005) and positive wave height anomalies have been forecasted (Dobrynin et al., 2015). All of these hydrometeorologic drivers influence shoreline retreat; however, changes in one parameter do not always result in changes in other parameters. In order for a chain reaction to occur they have to act simultaneously, amplifying each other (Baranskaya et al., 2017). For example, an increase in the sea ice-free period duration will only result in active coastal erosion if there are considerable storms coinciding with it. Such complex patterns make correlations of hydrometeorological parameters with temporal evolution of coastal erosion a challenging task (Jones et al., 2018; Shabanova et al., 2018).

The spatial variability of shoreline retreat in the Arctic is determined by coastal morphology, lithology, exposure, ground ice and permafrost properties (Jones et al., 2009; Konopczak et al., 2014; Maslakov and Kraev, 2016; Obu et al., 2016). In the West Siberian Arctic, one of the typical permafrost features is the presence of massive ground-ice beds associated with Late Pleistocene deposits (Streletskaya et al., 2001; Belova, 2015; Vasil’chuk, 2012), large layers of ground-ice inside frozen sediments. Contrary to the East Siberian Arctic, where deposits of the Ice Complex (ice-rich permafrost with numerous ice wedges) compose quickly eroding coastal plains (Lantuit et al., 2012; Günther et al., 2015), massive ice beds are a typical feature of the Yamal and Gydan peninsulas. They are often exposed in coastal bluffs, and their rapid thawing increases rates of retreat of ice-rich segments in comparison with the adjacent areas with little or no ground ice (Belova, 2015). Because of this, they often factor into the development of retrogressive thaw slumps (or thermocirques) (Leibman and Kizyakov, 2007; Lantuit et al., 2012; Khomutovm and Leibman, 2016; Belova et al., 2018) or promote subsidence of coastal plain surface (Lim et al., 2020). The spatial variability of coastal erosion within segments with different morphological and permafrost properties imposed on the temporal evolution of climate parameters creates the rapidly changing and highly variable spatial patterns of coastal dynamics in the Arctic.

The southwestern Kara Sea is among the best-studied polar regions in terms of coastal dynamics. Field observations of coastal erosion began in the 1980s along the coasts of the Baydaratskaya Bay (Kamalov et al., 2006; Novikova et al., 2018), near Kharasavey (Yuriev, 2009; Belova et al., 2017) and Marre-Sale settlements (Vasiliev et al., 2006; Vasiliev et al., 2011); additional remote sensing based studies followed (Kritsuk et al., 2014; Belova et al., 2017; Belova et al., 2018; Novikova et al., 2018; Isaev et al., 2019). They showed average long-term shoreline retreat rates ranging from 0.3 to 3.0 m yr⁻¹; these values are generally lower compared to Eastern Siberia, with its Ice Complex coastal bluffs retreating with rates from 0.9 to 9.0 m yr⁻¹ with extremes exceeding 15 m yr⁻¹ (Günther et al., 2013, 2015; Pizhankova, 2016; Grigoriev, 2019; Nielsen et al., 2020). However, these shoreline retreat rates are comparable to some sites along the Beaufort Sea with average shoreline retreat rates from the 1970s to the early 2000s ranging from 0.5 to 2 m yr⁻¹ (Brown et al., 2003; Radosavljevic et al., 2016; Irrgang et al., 2018), with extremes up to 9.4 m (Jones et al., 2018), increasing in the 2000s–2010s. In this light, shoreline retreat rates of Western Siberia, often exceeding 2 m yr⁻¹ are considerable taking into account that they are long-term average values, representing sites of several kilometers (Kritsuk et al., 2014; Novikova et al., 2018). We use the term “sites” identifying areas of coastal dynamics monitoring of several kilometers’ length with varying coastal morphology and other properties (e.g., Gulf of Kruzenstern, Kharasavey, etc.). By “segments” we mean short patches of coastline within these sites with relatively uniform morphology, exposure, permafrost properties, etc. All of the previously studied sites face toward the open sea and are directly exposed to waves of different directions. Wave energy is therefore one of the most important drivers of the erosion for these segments, as rates of coastal erosion showed an increase in years with longer sea ice-free period and stronger storms, an example being the period between 2005 and 2012–2013 for two stretches of coast of the Baydaratskaya Bay (Novikova et al., 2018; Isaev et al., 2019). However, so far there have been no attempts to quantify the relative role of abrasion driven by the wave factor, and thermal denudation, resulting from air temperature increase and permafrost thawing. One possible way to estimate the role of thermal denudation in eroding the Kara Sea coasts is comparing the results of monitoring at known open-sea sites with shoreline retreat rates of coasts in an enclosed bay with less wave action. In addition, observing spatial variability of coastal erosion in such a bay is a tool to link it to coastal morphology, sediment composition and permafrost properties and assess their influence on shoreline retreat.

In this study, we chose the shallow Gulf of Kruzenstern, Western Yamal, protected by the Sharapovy Koshki Islands, as a test site to compare to erosion studies being conducted along bluffs facing toward the open sea in the same region (Kritsuk et al., 2014; Belova et al., 2017; Novikova et al., 2018). Here, we address the questions related to the relative contribution of factors driving spatial and temporal differences in shoreline retreat rates across sections of coast that are exposed to direct wave action to different degrees. Using multitemporal satellite imagery, we estimate decadal-scale shoreline retreat rates since the 1960s for the Gulf of Kruzenstern, and compare the patterns of erosion in the Gulf to the spatial distribution of outcrops of massive ice beds observed in the field. We also analyze the temporal evolution of the hydrometeorological parameters since the 1970s, and compare them, along with shoreline retreat rates, to the areas of long-term coastal erosion monitoring situated from 70 to 150 km away from the gulf in open-sea settings in the Western Yamal study region (Kritsuk et al., 2014; Belova et al., 2017; Novikova et al., 2018). Our aim is to study the naturally isolated effect of thermal denudation and make the first attempt to quantitatively compare it to areas where its influence is combined with the effect of wave abrasion.

Study Area and Regional Setting

The Gulf of Kruzenstern is situated in the western part of the Kara Sea, near the western Yamal Peninsula (Figure 1). It is protected from the open sea from the west and northwest by a chain of low accumulative sandy barrier islands named Sharapovy Koshki Islands. The size of enclosed part of the lagoon between Yamal peninsula and the Sharapovy Koshki Islands is about 50 × 30 km; the deepest part is near its southern entrance, where depths reach slightly more than 10 m. Most of the bathymetry is very shallow, with average depths of 2–3 m. Areas of less than 1vm depth stretch 200–400 m into the Gulf. The eastern coast of the Gulf consists of several levels of terraces with heights of 10–15 and 25–40 m, with steep bluffs, and low floodplains and deltaic laidas (salt marshes) of not more than 1–2 m elevation. The floodplains stretch for 10–15 km.

FIGURE 1

FIGURE 1. Position of the study area and other key sites of coastal dynamics’ monitoring mentioned in the text. 1—Gulf of Kruzenstern, Eastern Coast; 2—Kharasavey, 3—Marre-Sale, 4—Ural coast of the Baydaratskaya Bay, 5—Yamal coast of the Baydaratskaya Bay. Marresalya hydrometeorological station at Marre-Sale is shown with an asterisk; Kharasavey hydrometeorological station is shown with a diamond. Abbreviations: GK—Gulf of Kruzenstern, CY—Cape Yasalya.

The climate of the region is cold, with a long winter lasting up to 9 months, and a cool, short summer. Mean annual air temperatures in the second half of the 20th century varied from −8 to −10°C (RIHMI-WDC, 2020). Mean annual accumulated positive air temperatures (thawing degree days) crucial for permafrost thawing and, hence, coastal dynamics, rose from 485 to 700°C-days in 1977–2013 at Kharasavey station, 70 km north of the Gulf (Belova et al., 2017), and from about 600 to 780°C-days in 1979–2013 at Marre-Sale, 90 km south of the Gulf (Shabanova et al., 2018). These values are smaller compared to the Barents Sea coasts, where they increased from ∼800°C-days to ∼1,000°C-days since the 1950s at Kolguev Island, and higher compared to stations at high latitudes both in the Kara and Laptev seas such as Dikson, (from ∼300 to ∼500°C-days since the 1970s), Terpiay-Tumsa or Shalaurova, (from ∼200 to ∼300°C-days since the 1960s) (Ogorodov et al., 2020). Average yearly precipitation varied from 300 to 350 mm in the 1950s–1970s (Trofimov et al., 1975) with a recent trend to increase to ∼370 mm in Amderma (Aleksandrov et al., 2005). The wind conditions of the Kara Sea are driven by the monsoon effect: in winter, the Siberian winter high pressure center determines southerly winds (Wu and Jia, 2002), while in summer, northerly and northwesterly winds prevail (Trofimov et al., 1975).

The thickness of Quaternary sediments in the western part of Yamal can reach 200 m (German et al., 1963). Terrestrial permafrost is continuous, with temperatures at the depth of zero annual amplitude varying from −5 to −10°C. The thickness of permafrost varies from 50 m at river mouths to 150 m in interfluvial areas. Depths of the active layer on the coastal plains varies from 0.6 to 1.5 m, depending on sediment composition and slope exposure (Melnikov et al., 2004). One of the typical features of permafrost is the presence of abundant massive ground-ice beds in the form of large horizontal ice bodies. The ground-ice outcrops in coastal exposures on the Ural coast of the Baydaratskaya Bay (Belova et al., 2008; Belova, 2015; point 4, Figure 1B), near Marre-Sale station (Tarnogradskiy, 1982; Astakhov et al., 1996; Forman et al., 2002; Kanevskiy et al., 2005; Kritsuk, 2010; Shpolyanskaya, 2015; Streletskaya et al., 2018, etc.; point 3, Figure 1B) and Cape Kharasavey (Belova et al., 2017; point 2, Figure 1B). The presence of similar massive ice beds at the Bovanenkovo gas site, documented by numerous boreholes and outcrops over a large area (Dubikov and Koreisha, 1964; Tarasov, 1990; Solomatin and Koniakhin, 1993; Baulin et al., 1996; Kondakov et al., 2001; Romanenko et al., 2001; Dubikov, 2002; Vasil’chuk et al., 2009; Vasil’chuk, 2012; Fotiev, 2012; Vasil’chuk et al., 2014) provides evidence that massive ice beds are typical for Western Yamal and are present under most of its plains, occasionally outcropping in coastal bluffs or thermocirques. The origin of these massive ice beds has a broad range of interpretations, from intrasedimental (Dubikov and Koreisha, 1964; Streletskaya et al., 2001; Streletskaya et al., 2009; Slagoda et al., 2010; Slagoda et al., 2012; Vasil’chuk, 2012; Streletskaya et al., 2013; Streletskaya et al., 2018) to buried glacial ice (Astakhov et al., 1996; Astakhov and Nazarov, 2010).

On the southwestern Kara Sea coasts, long-term, high-resolution monitoring of coastal dynamics has been conducted by several research teams. Instrumental measurements have been made since 1988 on the Ural and Yamal coasts of the underwater pipeline crossing of the Baydaratskaya Bay (Ogorodov et al., 2016; Kamalov et al., 2006; points 4 and 5, Figure 1B). Since 1981, they have been made near Kharasavey settlement (Belova et al., 2017; point 2, Figure 1B) and since 1978 near Marre-Sale (Vasiliev et al., 2011; point 3, Figure 1B). Application of multitemporal satellite imagery allowed for the determination of shoreline retreat rates for contiguous segments of the study regions and to extend the period of investigations to since 1969 for Marre-Sale (Kritsuk et al., 2014), since 1964 for the Yamal coast of the Baydaratskaya Bay and since 1968 for the Ural coast of the Baydaratskaya Bay (Novikova et al., 2018) because of the use of Corona satellite imagery from the 1960s.

Methods

Fieldwork

Field investigations on the morphology and sediment composition of the coastal bluffs of the eastern coast of the Gulf of Kruzenstern were carried out in July 2016. The topography and geomorphic processes were described at key points, distributed along the coast every 400 m; the total length of the explored shoreline was about 50 km. Geomorphic profiles were made; in separate outcrops, exposures of Quaternary sediments were cleaned and described. Samples for grain size analysis, diatom analysis and radiocarbon dating were collected. Sections up to 2 m wide were described. Cryogenic structure and permafrost properties were documented; outcrops of massive ice beds were described; their position was marked with a handheld GPS unit Garmin Etrex Vista Hcx with a precision of up to 5 m. The position of thermocirques was also marked along the coast with the same GPS unit. High rates of thawing resulting from extreme summer temperatures of 2016 (up to +16°C) and destruction of the thermoabrasional niches did not allow detailed sampling of the massive ice beds, as the overhanging 10–15 m of sediments threatened to collapse during sample collection and excavation at the outcrops.

Conventional radiocarbon dating was made in the Laboratory of Palaeogeography and geomorphology of Polar Countries and the World Ocean named after V. Köppen, Institute of Earth Sciences, Saint-Petersburg State University. Grain size analysis was done at the Institute of Geography RAS, Moscow; diatom analysis was done in VNIIOkeangeologiya, Saint-Petersburg.

Satellite Imagery

To estimate changes in the position of the shoreline since the second half of the 20th century, multitemporal imagery was used for two coastal segments of ∼20 and ∼17 km length. We investigated erosional segments only, as the shoreline position on low accumulative segments highly depends on tides and surges and is therefore challenging to detect on multitemporal imagery. The accumulative coasts excluded from the analysis made about 30% of the study site. The decision whether the coast was accumulative or erosional was made based on coastal morphology determined during fieldwork in 2016 and on the images of 2010 and 2019. We analyzed recent changes in the position of the shoreline using a WorldView-1 (WV-1) image of 2010 with a spatial resolution of 0.5 m and a WorldView-2 (WV-2) image of 2019 with a spatial resolution of 0.6 m (DigitalGlobe, 2020). The WorldView-1 2010 image initially had Standard 2A level of processing: it was radiometrically corrected, sensor corrected, projected to a plane using the WGS 84 datum, and had a coarse DEM applied to it by the provider, to normalize for topographic relief with respect to the reference ellipsoid (Digital Globe Core Imagery Product Guide, 2016). The WorldView-2 2019 image was provided as Basic 1B level of processing: it was radiometrically corrected and sensor corrected, but not projected to a plane using a map projection or datum. Both images have a 5 m geolocation accuracy. The WV-2 2019 image was georeferenced to the WV-1 2010 image using 30 ground control points and the spline transformation in ArcGIS 10.6.1™. The accuracy of the georegistration of the 2019 imagery was assessed with 50 additional ground check points located within 100 m of the shoreline. The mean uncertainty of georeferencing of the image from 2019 to the image from 2010 imagery was 0.94 m, with a minimum and maximum error of 0.26 and 1.60 m, respectively.

To determine the position of the shoreline in 1964, Corona KH-4 satellite images were used. Corona is declassified military satellite imagery from 1960 to 1972 with global coverage, distributed by the U.S. Geological Survey (EarthExplorer, 2020) as scanned film strips 70 mm × 75.6 cm in size in four image tiles with a 7 μm (3,600 dpi) scan resolution. We used images from the 9th of August, 1964, with a spatial resolution of 2.7 m. The parts of the Corona images covering the area of the study were referenced to the 2010 WorldView-1 image with a spline transformation with 52 referencing points. For more precise georeferencing, additional points referenced by handheld GPS in field were applied. The accuracy of the referencing of Corona imagery (δ_r), estimated by manual comparison of 50 points across the study area, is equal to 4.2 m.

The WV-2 image of 2019 and the Corona image of 1964 that initially had no projection were projected using the respective UTM zone (42N) on a WGS-84 datum.

We digitized the shorelines in 1964 and 2010 manually in the ArcGIS 10.2 (ESRI) software at a sub-pixel accuracy equivalent to the scale from 1:300 to 1:2000. The shoreline was similarly digitized in the 2019 scenes using ArcGIS 10.6 (ESRI). We traced the edge of the cliff at erosional segments only, as the accumulative segments are difficult to interpret because of unstable vegetation position and considerable shoreline shifts depending on the tide.

To calculate the shoreline position changes, we used the Digital Shoreline Analysis System (DSAS) of Thieler et al. (2017), available as an extension to ArcGIS. The program automatically creates transects perpendicular to a previously drawn baseline. The position of the intersection points between the transects and the shorelines are consequently being used for calculating shoreline change statistics. For this study, a transect spacing of 500 m was chosen which resulted in 74 transects. The DSAS program measures the shoreline retreat rates as an offset from the baseline L (t_i) for imagery acquisition dates t_i, and calculates statistics on the overall shoreline retreat rates and their uncertainty.

We calculated shoreline retreat rates for the time span between imagery acquisitions in m yr⁻¹:

r_{L} (t_{i - 1}, t_{i}) = Δ L / Δ τ = (L (t_{i} - L (t_{i - 1}))) / (t_{i} - t_{i - 1}), (1)

where r is the shoreline retreat rate, L is the distance of retreat, and t is time.

The volume of eroded deposits was assessed based on the obtained values of shoreline retreat (r_V) and altitudes of transects extracted from the ArcticDEM (Porter et al., 2018) (h_t) in cubic meters from 1 m of shoreline:

r_{V} = h_{t} \times r_{L} . (2)

Calculating the uncertainty of shoreline position on the WorldView-1 imagery, we took into consideration the topography-induced horizontal displacement (Toutin, 2004), which is

δ_{t} = t g α \times H, (3)

where α is the tilt angle of the spacecraft (14.1° in our case), and H is the maximum relative height at the territory (25 m), resulting in 6.25 m.

The uncertainty of shoreline position is calculated as:

δ_{x} = \sqrt{δ_{t}^{2} + δ_{s}^{2}} (4)

for the WorldView-1 base image. For the referenced WorldView-2 2019 and KH-4 Corona 1964 images, the uncertainty of shoreline position was calculated as:

δ_{x} = \sqrt{δ_{r}^{2} + δ_{s}^{2}}, (5)

where δt is the topography-induced uncertainty, δs is the uncertainty of shoreline digitization, which is one half of the imagery spatial resolution, and δr is the uncertainty of the referencing of Corona and WorldView-2 images to WorldView-1. The total uncertainty of the shoreline position results in 6.25 m for the WorldView-1 of 2010, 0.99 m for the WorldView-2 of 2019 and 4.46 m for the Corona imagery of 1964.

The uncertainty of the rate of linear shoreline retreat from x₁ to x₂ during the period from t₁ to t₂ was calculated as:

δ_{L} = \frac{\sqrt{δ x_{1}^{2} + δ x_{2}^{2}}}{t_{1} - t_{2}}, (6)

and results in 0.2 m yr⁻¹ for 1964–2010, 0.7 m yr⁻¹ for 2010–2019 and 0.1 m yr⁻¹ for 1964–2019.

The volumetric approach utilizes heights of eroded surface (h) as well as rates of planar retreat (r), so the total uncertainty of volumetric rates δ_V comprises:

δ_{V} = \sqrt{h^{2} δ_{L}^{2} + r^{2} δ_{h}^{2}}, (7)

where δ_h is the uncertainty of estimation of the eroded surface height, depending on the surface roughness and equal to the height variance of the eroded surface at the current cliff.

Hydrometeorological Forcing

To estimate the changing hydrometeorological impact, the evolution of the air temperature, sea ice-free period duration and wind-wave energy was estimated. The potential of thermal denudation was assessed by the air thawing index (I_t) showing the number of positive/negative °C∙days per year (Andersland and Ladanyi, 2004):

I_{t} = \sum_{i = 1,}^{N} t_{i}, t > 0, (8)

where $t_{i}$ is the daily mean temperature, and $N$ is the number of warm (t > 0°C) days per year. A similar parameter called positive degree-day (PDD) sums was used in Günther et al. (2015). This index is an evaluation of the annual amount of heat added to the ground and permafrost during warm periods. For its calculation, we used both observational data from Marre-Sale (Marreslaya hydrometeorological station) and Kharasavey (Kharasavey hydrometeorological station, Figure 1C) and ERA Interim reanalysis data on air temperature at a height of 2 m (Berrisford et al., 2009).

The Kharasavey hydrometeorological station was established in 1972; observations were made until 1986. We used the whole obsevation period. Data from Marre-Sale has gaps in 1973–1976, as well as in the 2000s. For the Gulf of Kruzenstern, the nearest cell of ERA Interim grid (1979–2018) with a grid resolution of 0.75° was used. Calculation of the air thawing index (It) at Marre-Sale and Kharasavey was done according to the procedure described in Shabanova et al. (2018) and Belova et al. (2017). To reconstruct thermal conditions in the Gulf of Kruzenstern, I_t from Marre-Sale and Kharasavey were weighted depending on their distance to the Gulf of Kruzenstern. The weights were calculated as the inversed linear distances:

t_{k} = \frac{t_{m} \cdot \frac{1}{d_{k m}} + t_{h} \cdot \frac{1}{d_{h k}}}{\frac{1}{d_{k m}} + \frac{1}{d_{h k}}} = t_{m} \cdot \frac{d_{h k}}{d} + t_{h} \cdot \frac{d_{k m}}{d}, (9)

where $t_{h}$ , $t_{m}$ and $t_{k}$ are temperatures in Kharasavey, Marre-Sale and Kruzenstern respectively, d is the distance between Marre-Sale and Kharasavey ( $d = d_{h k} + d_{k m}$ ), where $d_{h k}$ is the distance from Kharasavey to Kruzenstern, and $d_{k m}$ is the distance from Kruzenstern to Marre-Sale. As $d_{k m}$ is equal to ∼100 km, and $d_{h k} is equal to \sim$ 70 km, the weights $\frac{d_{h k}}{d}$ and $\frac{d_{k m}}{d}$ make 0.4 for Kharasavey and 0.6 for Marre-Sale, respectively.

We chose eleven years in which data from both stations were available (1972 and 1977–1986). The obtained series were used to build the regression model with the ERA Interim data in the nearest “land” cell of the ERA grid, the size of the grid being equal to approximately 0.75°. The calculated regression coefficient was 0.82 and the determination coefficient was 0.94, showing sufficient agreement between the reanalysis and observational data.

The duration of the sea ice-free period was determined based on satellite imagery data. We used the sea ice concentration products of the EUMETSAT Ocean and Sea Ice Satellite Application Facility (OSI SAF) OSI-450 (EUMETSAT, 2017) for 1979—2015 and OSI-430-b for 2016–2018 (EUMETSAT, 2019). The products are provided by the Danish Meteorological Institute; their resolution is 25 km. Previous works based on OSI SAF products showed that the accuracy of the sea ice-free period start and end dates determination is comparable to NSIDC (Shabanov and Shabanova, 2019).

We determine the start and end dates of the sea ice-free period for the cells of OSI SAF data nearest to Marre-Sale, Kharasavey, the Gulf of Kruzenstern eastern coast and to the outer coast of the Sharapovy Koshki islands outside the Gulf of Kruzenstern. To detect the sea ice-free period start and end dates, the original rolling-window method (rolling-window approach or RWA) is used. The method is based on sea ice concentration (SIC) annual variation analysis and detects the robust SIC jumps, which are usually linked to water clearing-up (the last date of the downward jump) or freezing-up (the first date of the upward jump). The details of the RWA are published in Shabanov and Shabanova (2019). This method was introduced to avoid the disadvantages of the common 15%-method which is sometimes inapplicable in coastal zones due to SIC data contamination from land which results in high (20–40%) SIC values during the open water season.

The wave energy potential (WE), called “wave-energy flux” in earlier studies (Popov and Sovershaev, 1982), was calculated according to the Popov–Sovershaev method (Popov and Sovershaev, 1982; Ogorodov, 2002; Shabanova et al., 2018). The method is based on the wave processes theory and applies correlations between wind speed and parameters of wind-induced waves.

For deep-water conditions, when the sea floor does not affect wave formation, the wave energy potential, coming from 1 m of the wave front per second during a storm with the wind speed $V_{j i}$ with a wind speed class j from a chosen wave direction i at the outer coastal zone boundary is proportional to the wind speed to the power of three and to the wave fetch:

W E_{j i} = 3 \times 10^{- 6} V_{j i}^{3} x_{i}, (10)

where $V_{i j}$ is the real wind speed of a chosen direction measured by anemometer at 10 m above sea level (m/s), $x_{i} .$ is the wave fetch (km) along the current wind direction. The $3 \times 10^{- 6}$ coefficient is derived from the wind speed and wave parameters’ empiric ratios and corresponds to the dimension of $ρ / g$ , where $ρ$ is density (kg/m³), and $g$ is gravitational acceleration (m/s²). Thus, $W E_{j i}$ has the dimension of kilograms per second:

\frac{k g}{m^{3}} \cdot \frac{s^{2}}{m} \cdot \frac{m^{3}}{s^{3}} \cdot m = \frac{k g}{s} . (11)

The total wave energy potential, accumulated during the ice-free season, is calculated as:

W E = \sum_{i = 1}^{M} \sum_{j = 1}^{N} 3 \times 10^{- 6} V_{j i}^{3} x_{i} \cdot p_{j i} \cdot n, (12)

where $p_{j i}$ is the frequency of storm with the wave direction i and wind speed class j (calculated for the ice-free season only), n is the sea ice-free period duration in seconds, M is the number of wave-dangerous directions, and N is the number of wind speed classes. The wave energy potential is then expressed in kilograms (kg per season or kg yr⁻¹) and is the measure of the wave energy coming to the coastal zone, expressed as mass of incoming water.

Landward wave-dangerous directions of winds, wave fetches and depths were obtained from ETOPO1 digital elevation model, which also includes bathymetry (NOAA, 2017). The vertical resolution of ETOPO1 did not allow to take into account Sharapovy Koshki islands, because of their shallow topography being lower than the resolution of the DEM. For this reason, the wave energy potential evolution was calculated for the outer coast of the Gulf of Kruzenstern. The wave energy potential inside the Gulf was calculated separately using topographic maps. The frequency of wave-generating winds was estimated for the sea ice-free period only; winter months were not taken into account. For wind data, ERA Interim reanalysis was used (Berrisford et al., 2009); we took data on the wind speed and direction at a height of 10 m. The complete results of calculations of the air thawing index, ice-free period duration and wind-wave energy for the Gulf of Kruzenstern, its outer border at Sharapovy Koshki Islands, Kharasavey and Marre-Sale can be found in the Supplementary material (Table S1)

To estimate the relative role of wave destruction and thermal denudation in the total hydrometeorological forcing, the distribution of the values of the wave energy potential and air thawing index for four sites: Marre-Sale, Kharasavey, the Gulf of Kruzenstern eastern coast and to the outer coast of the Sharapovy Koshki islands outside the Gulf of Kruzenstern, was analyzed. The 0.975 quantile for the wave energy potential of all the four sites together made from 0 to 1,865 × 10⁶ kg per year. The 0.975 quantile for the air thawing index made from 0 to 1,065 degree-days. All the values of WE and I_t were divided by these two values, respectively. Therefore all the transformed values had the range from 0 to 1 and no dimension. The transformed WE (WE_t) and I_t (It_t) values were further summed to calculate total hydrometeorological effect (stress or action) TE acting on the coasts:

T E = W E_{t} + I t_{t} . (13)

Based on these values, the relative role of the wave energy potential and air thawing index in the total hydrometeorological effect for each site could be quantitatively estimated. This estimation is relative, as it depends on the choice of sites for analysis, and shows the relation of the two factors between these sites only. It therefore helps to compare the conditions of neighboring sites only.

Results

Quaternary Sediments and Massive Ice Beds

The sections exposed in the bluffs of the Gulf of Kruzenstern generally consist of three main units: Unit 1) brown silts at the top; Unit 2) poorly sorted silty thin sands in the middle and Unit 3) dark gray silts and clays with massive ground-ice beds at the bottom (Figure 2).

FIGURE 2

FIGURE 2. General section of Quaternary sediments outcropping in the coastal bluffs of the Gulf of Kruzenstern. 1—sands, 2—silts, 3—peat, 4—ice, 5—geological interfaces. (A) Brown silts of the upper layer, underlain by peat; poorly sorted thin sands of the top of the middle layer. (B) Erosional contact between sands of the middle unit and silts and clays of the lower unit; deformations in the upper part of the lower unit. (C) Dark gray silts and clays of the bottom unit. The red box shows the position of material taken for radiocarbon dating; the radiocarbon age, error and sample number are shown in red color.

The upper unit (Figure 2, Unit 1) consists of brown silts which are non-laminated and deformed by cryoturbation, with the thickness of 1–2 m; they cover all the high plains of 10–35 m elevation. Their lower contact with Unit 2 is straight and horizontal, sometimes marked by thin (up to 5 cm) layers of peat debris. Separate benthic freshwater alkaphile diatom algae Staurosira venter were encountered, along with the extinct re-deposited diatoms Paralia grunowii, Aulacoseira praegranulata, fragments of Grunowiella gemmata and non-identifiable fragments of marine centric diatoms. Thin (2–3 cm) layers of peat at the lower boundary (Figure 2A), radiocarbon dated to 8,160 ± 80 years (LU-8457; 9,130 ± 120 cal. Years), give evidence that the upper brown silts formed in the Holocene; they are cover sediments accumulated in terrestrial conditions, derived from the presence of single freshwater diatoms and cryogenic deformations by multiple thawing and re-freezing events.

The middle unit (Figure 2, Unit 2) consists of poorly sorted silty thin sands with uneven lamination and traces of thawing and re-freezing. Separate freshwater boggy diatoms Eunotia parallela were documented, as well as single valves of the extinct Palaeogene Paralia grunowii, fragments of non-identifiable marine centric diatoms and sponge spicula. The presence of single freshwater diatoms, absence of distinct lamination implying sedimentation from water, and traces of thawing and re-freezing suggest sub-aerial or alluvial accumulation. This unit overlies the lower silts and clays of Unit 3 with an unconformity; its lower contact goes up and down, sometimes with a steep dip up to 60^o (Figure 2B).

The lower unit (Figure 2, Unit 3) of dark gray and bluish gray silts and clays is exposed at elevations of up to 5–6 m a.s.l.; its bottom lies below sea level. The silts and clays are either non-laminated or have thin horizontal layers of gray sandy loam with deformed initial lamination; they have a splintered texture (Figure 2C) and a reticulate cryogenic structure, with up to 2 cm thick ice lenses. The upper part of the silts and clays contains layers of brownish sandy loams; this part is folded. The sediments contain few microfossils, including separate valves of the extinct diatom algae Paralia grunowii, non-identifiable fragments of marine centric diatoms, sponge spicula and separate cysts of freshwater Chrysophyte algae.

Massive ice beds, incorporated into the dark gray silts and clays, were noted in several thermoabrasional niches at elevations of up to 2–3 m a.s.l. The niches are steep, with an overhang of clays and silts above the ice. In the warm summer of 2016, the niches were very unstable: every few minutes, large volumes of thawed silts collapsed over the quickly melting ice, suggesting a very short lifetime of every massive ice exposure of not more than a few days. The shape and position of the ice outcrops probably changes every year, with some niches collapsing, and other ones becoming exposed.

The length of the ice exposures along the coast varies from 3 to 10 m; the visible thickness of ice is up to 2.5 m; its bottom is below sea-level. The ice is clean and transparent, with no bubbles and large (up to 2 cm) crystals. The massive ice beds have distinct layers of 0.3–0.5 m thickness; they form micro-steps in the outcrop. The upper contact of the ice with clays cuts these layers (Figure 3). From the contact, ice veins of up to 2 cm thick go up and sideways, creating the reticulate cryogenic structure in the dark gray silts and clays. Around the ice, they are non-laminated; layers of brown sandy loam are not seen. At 0.2–0.3 m above the massive ice beds with transparent clean ice, lenses of dirty turbid ice of up to 0.3 m thickness occur.

FIGURE 3

FIGURE 3. Coastal exposure of massive ice beds near Cape Yasalya. The dashed line shows the boundaries of the ice body.

Massive ice beds mostly outcrop in steep erosional niches where the shoreline is relatively straight. However, they also exist under large thermocirques, spread along the shoreline, being the main reason of their formation. The ice is usually buried under the thawed grounds covering the flat bottom of the thermocirques. On the one hand, thawing of the ice, containing almost no sediments leads to considerable ground subsidence. On the other hand, the water saturates the sediments above the contact, while the remaining ice acts like an aquiclude, leading to massive landslides streaming into the Gulf. The size of the thermocirques varies from 10 to 20 to 100–150 m in diameter; fresh recently formed thermocirques occur along with older ones. In the north of the area, the plant cover of a relatively old landslide by mayweed (Matricaria inodora L.) suggested its formation about 2 years ago, after which the landslide remained inactive. The distribution of the thermocirques and outcrops in niches (Figure 4) suggests that the massive ice beds are present almost everywhere along high coasts near Cape Yasalya.

FIGURE 4

FIGURE 4. Thermocirque to the north of Cape Yasalya: (A) View from space in 2018 (Yandex Digital Maps, 2020 https://yandex.ru/maps/-/CCVtMVJd). (B) Field view in 2016; the cliff height is ∼15 m.

Rates of Coastal Erosion

During the 46-years time period between 1964 and 2010, the eastern erosional coasts of the Gulf of Kruzenstern retreated on average by 23.3 ± 9.2 m (Table 1), resulting in an average shoreline retreat rate of 0.5 ± 0.2 m yr⁻¹. For active thermoabrasional segments not covered by vegetation in 2016, mean shoreline retreat rates reached 0.9 ± 0.2 m yr⁻¹. The area north of Cape Yasalya had the highest shoreline retreat rates of up to 1.5–1.7 ± 0.2 m yr⁻¹ concentrated along 2 km of the shoreline facing to the northwest (Figure 5). Volumetric rates show more variability than the linear shoreline retreat rates; however, their values also increase toward the segments north of Cape Yasalya, reaching extremes up to 27.8 ± 3.1 m³ m⁻¹ yr⁻¹. The highest retreat of up to 80.0 ± 9.2 m was observed between the two largest outcrops of massive ice beds exposed in summer 2016 (Figure 6). Because of considerable planimetric retreat, the highest volumes of sediment were eroded from the terrace with heights of 10–15 m immediately north from Cape Yasalya, and not from the high coastal bluffs of 25–40 m immediately south of Yunetayakha River or south of Mordiyakha River (Figure 5).

TABLE 1

TABLE 1. Parameters of shoreline retreat for erosional coastal segments of the Gulf of Kruzenstern.

FIGURE 5

FIGURE 5. Rates of retreat of the eastern coast of Gulf of Kruzenstern (area 1, Figure 1C), for 1964–2010 and 2010–2019. The dots indicate rates of shoreline retreat measured at every 500 m. Gray dots (“no data”) correspond to areas covered by clouds where it was impossible to calculate retreat rates. “T10–15”—coastal plain (terrace) with heights of 10–15 m; “T25–40”—coastal plain (terrace) with heights of 25–40 m. The position of Figure 6 is shown as a white rectangle.

FIGURE 6

FIGURE 6. Shorelines of 1964 and 2010 within the area of the fastest erosion and massive ice beds’ exposure: (A) Corona image. (B) WorldView-1 image. 1—outcrops of massive ice beds in 2016, 2—outcrops of thermocirques in 2016, 3—shoreline of 2010. 4—Shoreline of 1964.

During the 9-years time period between 2010 and 2019, coastal erosion accelerated. The average retreat rate of the whole shoreline almost doubled: 0.5 ± 0.2 m yr⁻¹ in 1964–2010 compared to 0.9 ± 0.7 m yr⁻¹ in 2010–2019 (Table 1), although the increase in the uncertainty of the rates’ calculation was considerable because of a shorter time period between observations. Peak erosion rates also increased by 40% (2.4 ± 0.7 m yr⁻¹ compared to 1.7 ± 0.2 m yr⁻¹ in 1964–2010). Overall, in nine years, the coast retreated by 21.9 ± 6.3 m. Volumetric retreat rates almost doubled making up to 41.1 ± 11.6 m³ m⁻¹ yr⁻¹. At the same time, the spatial variability of shoreline retreat increased, evidenced by higher RMSE values in 2010–2019 compared to 1964–2010 (Table 1). Segments with shoreline retreat rates exceeding 2 m yr⁻¹ were documented not only north of Cape Yasalya, where most of the massive ice beds outcrop, but also to the south of it, and south of Mordiyakha River. The last two areas had relatively low erosion in 1964–2010, which became more active in the last nine years. The rates in 2010–2019 can be partly biased by their large uncertainty; however, a trend in the increase of both, the average and highest shoreline retreat rates, and considerable erosion at previously calm segments together suggest an intensification in erosion in the last decade.

During the entire 55-years observational period, from 1964 to 2019, mean erosion rates were 0.5 ± 0.1 m yr⁻¹ on average, with extreme values up to 1.6 ± 0.1 m yr⁻¹, reflecting the effect of smoothening the extreme shoreline retreat rates through time: the longer period that is analyzed, the lower the mean erosion rates generally are. The greatest distance of shoreline retreat of 86.9 ± 5.5 m was typical for segments north of Cape Yasalya, between the two largest exposures of massive ice beds noted in the field (Figure 6). Areas with outcrops of massive ice beds and thermocirques marked in 2016 are generally destroyed faster than the adjacent segments; however, the outcrops do not exactly coincide with peak rates (Figure 7). Higher cliffs generally retreat faster than the lower segments, which is seen both in linear and volumetric erosion rates.

FIGURE 7

FIGURE 7. Distribution of erosion rates in 1964–2019, cliff heights and position of outcrops of massive ice beds (summer 2016) along the shoreline. The y axes show (from left to right): bluff height, planimetric coastal erosion rates, volumetric coastal erosion rates; the x axis shows the number of transects from the south to the north. 1—linear erosion rates, 2—volumetric erosion rates, 3—cliff height, 4—outcrops of massive ice beds exposed in July 2016.

Evolution of Hydrometeorological Drivers of Coastal Dynamics

Along the coast of the Gulf of Kruzenstern, the air thawing index in 1979–2018 increased considerably (Figure 8A): its trend rose by more than 300 degree-days. Comparison of this parameter to Marre-Sale and Kharasavey shows that, according to the latitudes of these three points, temperatures at the Gulf of Kruzenstern are higher compared to Kharasavey and lower compared to Marre-Sale; therefore there are more positive degree days than at Kharasavey and less positive degree days than at Marre-Sale. From 1979 to 2000, the air thawing index showed no clear trend, fluctuating from 300 to slightly more than 800 degree-days. The coldest periods tend toward the end of the 1970s and the end of the 1990s; warm years were 1982–1985 and 1993–1995. From the beginning of the 2000s, a steady temperature rise was observed; the air thawing index reached its peak value of 1,358 degree-days in 2016.

FIGURE 8

FIGURE 8. (A) Air thawing index (I_t) according to observational data and ERA Interim reanalysis at the Gulf of Kruzenstern (1), Kharasavey (2), and Marre-Sale (3). (B) Sea ice-free period duration (IFD) at the Gulf of Kruzenstern eastern coast (1), Kharasavey (2), Marre-Sale (3) and outer side at Sharapovy Koshki Islands (4); (C) Wind-wave energy potential (WE) evolution at the Gulf of Kruzenstern eastern coast (1), Kharasavey (2), Marre-Sale (3) and outer side at Sharapovy Koshki Islands (4). The trends for the Gulf of Kruzenstern are shown as dotted lines together with their slope (K), determination coefficient (R²) and p-value; (D) map showing the four sites for which the hydrometeorological parameters were calculated.

The duration of the sea ice-free period also increased in the Gulf of Kruzenstern and adjacent areas (Figure 8B), as a result of both earlier breakup of the ice in spring and its later formation in autumn. This increase coincides with rising air temperatures in the last decades, although there could be a lag in ice conditions from temperature changes. Despite considerable fluctuations, a steady rise since the 1980s is seen at all sites. During this time, the sea ice-free period increased by more than 50 days. Contrarily to the air thawing index, the longest sea ice-free period was observed at Kharasavey and at the outer side of the Gulf of Kruzenstern (Sharapovy Koshki islands), situated in the northern and central part of the area (Figure 1). The longest sea ice-free period of 205 days was noted in 2012 at Kharasavey. At the outer side of the Gulf of Kruzenstern, near Sharapovy Koshki islands, the sea ice-free period on average lasted five to ten days longer compared to the coast inside the Gulf. In 1982, the difference exceeded 20 days. The enclosed lagoon becomes covered by ice faster than the open sea, and the ice holds for a longer time in the early summer. Many shallow areas freeze to the bottom, resulting in later ice breakup. In this way, the waves erode the coast of the Gulf of Kruzenstern for a shorter duration than they do at open-sea segments.

Fluctuations of the wind-wave energy potential (Figure 8C) generally agree with changes in the sea ice-free period duration and show a slight increase in storm intensity. Similarly to the ice-free period, the wind-wave energy potential at Kharasavey, the outer side of the Gulf of Kruzenstern (Sharapovy Koshki) and Marre-Sale shows an increase in 1995–1995 and between 2007 and 2013. The values in all three areas agree well; the wind-wave energy potential inside the Gulf of Kruzenstern is more than five times smaller than at the outer rim of the Sharapovy Koshki Islands. Such difference results from smaller wave fetch on the one hand, and shorter sea ice-free period on the other hand.

The total hydrometeorological forcing (TE) also shows a general increase since 1979 for all four sites (Figure 9). The graphs show similar distributions with the same highs and lows, although the values of these highs and lows differ from site to site. For all open-sea sites (Figures 9B–D), the increase of TE was higher than for the Gulf of Kruzenstern.

FIGURE 9

FIGURE 9. Evolution of the normalized air thawing index (I_t) and wind-wave energy (WE), and their total effect (TE) showing the hydrometeorological forcing of coastal erosion of the Gulf of Kruzenstern (A), Kharasavey (B), Marre-Sale (C), and Sharapovy Koshki Islands (D). The trend of the total effect is shown as a dotted line. See Figure 8D for the map of the four sites.

In 1992 and 1999, at all sites TE was lowest, resulting from relatively cool summers. The two positive peak years were 1994 and 2016. The reasons of increased TE were different: in 1994, the temperature was close to its average values, while the wind-wave energy showed a considerable increase, increasing mechanical erosion. On the contrary, in 2016, the wind-wave energy was slightly lower than usual, and TE mainly resulted from high air temperatures and thus intensified thawing of the frozen grounds.

Generally, the highest TE was recorded at Marre-Sale (Table 2) where its long-term mean was close to 1.3. The values for Kharasavey and the Sharapovy Koshki Islands outside the Gulf of Kruzenstern are slightly smaller (1.06 and 1.15 respectively), while inside the Gulf TE is 30% lower than at the Sharapovy Koshki Islands. The total forcing inside the Gulf is smaller than at open-sea locations because the Gulf of Kruzenstern receives 75% less wave energy as it is shielded by the Sharapovy Koshki islands. Therefore, variations in the air thawing index inside the Gulf have higher impact on coastal erosion relatively to other sites: the relation of I_t to the total effect is equal to 85%, while at other sites it makes 51–53% (Table 3). It should be taken into consideration, however, that the obtained ratios were based on the assumption that for open-sea segments, the thermal and wave factor have equal weights, which could in fact be more complex because of the ice content, cliff composition, morphology and other factors. Consequently, the percentage of coastal erosion attributed to the thermal denudation in the Gulf of Kruzenstern is a relative assessment in comparison to the other three sites.

TABLE 2

TABLE 2. Comparison of shoreline retreat at key sites of monitoring in the southwestern Kara Sea.

TABLE 3

TABLE 3. I_t and WE long-term values transformed to 0–1 scale by dividing them by the 0.975 maximum values for four sites, their sum (TE), and the relative role of thermal denudation in coastal destruction (I_t-t/TE).

Discussion

Drivers of Spatial Variability in Coastal Erosion

Results of shoreline retreat measurements using remotely sensed data show considerable spatial variability in erosion rates: from 0.1 to 1.6 ± 0.1 m yr⁻¹ or from 0.2 to 24.2 ± 2.3 m³ m⁻¹ yr⁻¹ (1964–2019) at different segments, despite similar climatic and wave conditions within the eastern coast of the Gulf. The patterns of erosion rate distribution, showing faster retreat of the segment north from Cape Yasalya compared to the adjacent segment south of Cape Yasalya, with similar cliff heights, sediment composition and permafrost properties, indicate that the orientation of the coastal segment toward the potentially longest fetch and the most common wind direction is important for erosion intensity. Despite the small wave fetch in the Gulf of Kruzenstern, resulting in relatively low wind-wave energy potential compared to Kharasavey, Marre-Sale or to the coast seawards from the Sharapovy Koshki Islands, the most actively retreating coastal segment is the one that faces westerly and northwesterly winds. This gives evidence that even small amounts of wave energy can considerably influence coastal dynamics, depending on the coastal segment exposure toward waves. Such effect on larger scales was noted for Baydaratskaya Bay: while the Yamal coast faces to the southwest, and its average erosion rates were 0.3 ± 0.16 m yr⁻¹ in 1968–2016, the Ural coast, exposed to northerly and northeasterly winds from the open sea was characterized by average shoreline retreat rates of 1.2 ± 0.15 m yr⁻¹ in 1964–2016 (Novikova et al., 2018). The orientation of the Ural coast toward the open sea leads to the reception of more wave energy than the Yamal coast (Shabanova et al., 2018). Our findings imply that even for small enclosed water areas, the orientation and exposure effect remains significant and can directly influence long-term coastal erosion patterns.

Another factor determining the spatial variability of erosion rates are permafrost properties (sediment composition and the presence of ground ice and massive ice beds). The massive ice beds exposed along the coast are similar to massive ice beds found near the Bovanenkovo settlement in terms of morphology and ice properties (Dubikov and Koreisha, 1964; Solomatin and Koniakhin, 1993; Vasil’chuk et al., 2009; Vasil’chuk et al., 2014). In Bovanenkovo, the massive ice beds are 1–2–28.5 m thick and are incorporated in silts of 2–12 m thickness or clays of 3–9 m thickness; the general composition of these silts and clays is similar to the ice containing dark gray silts of the Gulf of Kruzenstern. The ice is considered either as polygenetic: segregation, infiltration-segregation, injection-segregation, lacustrine, river and sea ice (Vasil’chuk, 2012) or buried glacial ice, relic from a land-based ice sheet (Solomatin and Koniakhin, 1993).

Despite varying interpretations of the origin of the massive ice beds, it is generally agreed that the massive ice beds are widely distributed under the lowlands in the central and western Yamal (see, e.g., Vasil’chuk, 2012). As the massive ice beds of the Gulf of Kruzenstern are similar to Bovanenkovo, they should also be largely present under the coastal plains of the Gulf, as confirmed by borehole data in the area (Baulin et al., 2003). The presence of numerous thermocirques confirms this assumption as their presence is indicative of ice-rich permafrost deposits. The distribution of erosion rates relative to the outcrops of the ground-ice (Figure 7) show that the peaks do not directly match with the exposures. However, according to field observations, segments of increased shoreline retreat rates generally correspond to areas with abundant ground-ice outcrops and thermocirques (Figure 7). As the thermoabrasional niches with the outcrops were actively collapsing in summer 2016, we suppose that the life span of such a niche is relatively short; they were carved and became buried many times in many locations along the coast between the acquisitions of the three satellite images of 1964, 2010, and 2019. Another important factor is the high ground-ice content of the silts and clays above the massive ice beds. A case study of the Kharasavey settlement gives evidence that segments where the coastal bluff is composed by ice-rich clays can retreat faster than segments with abundant massive ice beds (Belova et al., 2017). Not only the ice bodies themselves, but also silts and clays with ice-rich reticulate cryostructure can contribute to considerable rates of retreat of the segment north of Cape Yasalya (e.g., Baulin et al., 2003). Wide presence of massive ice beds generally enhances erosion in the region, although every single outcrop does not necessarily imply peak erosion rates in its every location.

In order to gain insight into regional coastal erosion patterns, we compared shoreline retreat rates from the Gulf of Kruzenstern to rates from coastal monitoring stations on the southwestern Kara coast (Table 2). The highest average shoreline retreat rates were recorded at the Marre-Sale station and on the Ural coast of the Baydaratskaya Bay, where also the largest outcrops of massive ice beds are present (Kritsuk et al., 2014; Forman et al., 2002; Slagoda et al., 2012; Vasil’chuk, 2012). The average rates of the Gulf of Kruzenstern shoreline retreat are smaller; however, they exceed average shoreline change rates for the Yamal coast of the Baydaratskaya Bay, which is facing the open sea but does not contain massive ice beds. Comparing the long-term maximal shoreline retreat rates, we can conclude that the peaks for the Gulf of Kruzenstern were close to the values for the Ural coast, Marre-Sale and Kharasavey, and exceeded values for the Yamal coast. The average shoreline retreat rates along the Gulf of Kruzenstern in 1964–2010 (0.5 ± 0.1 m yr⁻¹) were lower compared to the Ural coast (1.2 ± 0.15 m yr⁻¹ in 1964–2012, Novikova et al., 2018). However, in the last decade, these rates were comparable: 0.5 ± 0.3 m yr⁻¹ in 2012–2016 at the Ural coast (Novikova et al., 2018), against 0.9 ± 0.7 m yr⁻¹ in 2010–2019 at the Gulf of Kruzenstern. Higher rates at the Gulf of Kruzenstern in the most recent period partly result from the warm summer of 2016, the effect of which might not be included in the observations at the Ural coast. Nevertheless, due to the increases in the air thawing index, erosion of ice-rich coasts in the Gulf may occur at similar speed to the more exposed and wave-dependent open-sea sites in the region. At the same time, average rates over longer timescales exceed rates for the Yamal coast, but are lower in comparison to all the other open sea sites. This leads us to a conclusion that long-term dynamics are also influenced by the wind-wave energy and sea ice-free period duration.

For the Gulf of Kruzenstern, the presence of massive ice beds has a major impact on the erosional behavior of the coast, since rising air temperatures are leading to ground ice melt which in turn evokes that the coastal bluff faces degrade and slump. Even though wave development is hampered by the Sharapovy Koshki Islands which are sheltering the study site from the open sea, waves developing inside the Gulf are sufficient to remove the debris from the collapsed coastal bluffs. Similar patterns were noted by Vasiliev et al. (2006), suggesting that for gulfs and bays, sediment composition and permafrost properties are more important than the wave energy in these settings. Therefore, coastal erosion in the Gulf of Kruzenstern is mainly attributable to the high presence of ground ice, rather than to the exposure to high wave energy.

Another prominent feature of the Gulf of Kruzenstern is the increase in the spatial variability of erosion rates in 2010–2019 compared to 1964–2010. In other locations in the Arctic, such as Drew Point on the Beaufort coast of Arctic Alaska, the distribution of coastal erosion rates became more uniform along with their recent increase (Jones et al., 2009): different geomorphic levels composed by sediments of different ice content were eroding at nearly identical rates in 2002–2007, contrarily to the previous periods from 1955 to 2002. There seems to be a regional difference in these patterns, as other areas of long-term coastal erosion monitoring in the Kara Sea show more variable distribution of shoreline retreat rates in recent years, just like the Gulf of Kruzenstern. An example of this is the Ural coast of the Baydaratskaya Bay (Figure 1B; point 4), where shoreline retreat rates in 2005–2012 and 2012–2016 had higher spatial variability across different landforms compared to 1964–1988 (Novikova et al., 2018). Erosion became especially active on the 6–8 m terrace and got slower on the 12–16 m terrace in the recent years, while in 1964–1988, the 12–16 m terrace experienced higher retreat. This effect needs more investigations with higher temporal resolution than the present study, however, suggestions can be made that the increasing spatial variability of erosion at the Gulf of Kruzenstern is caused by the specific distribution of permafrost with large massive ice beds and relatively small and rare ice wedges. As the top of the massive ice beds typically has its highs and lows (Lim et al., 2020), their exposure usually differs from year to year, as the shoreline retreats. If the top of the massive ice is close to the surface, it crops out along the cliff, and erosion rates accelerate locally. If erosion rates are averaged over long-term periods (20 years or more), the effect of locally accelerated erosion due to ground ice presence is evened out, as massive ice beds are present in the sediments composing all topographic levels above 5–6 m. In case of shorter observational periods, spatial variations in massive ice presence along the coast can influence the speed of erosion and thus the distribution of erosion rates, leading to a high spatial variability.

Drivers of Temporal Variability in Coastal Erosion

The shoreline of the Gulf of Kruzenstern generally retreats faster than the shorelines of other gulfs and bays across the Kara Sea: 0.2–0.7 m yr⁻¹ in all bays and gulfs (Vasiliev et al., 2006) against 0.5 ± 0.1 to 0.9 ± 0.7 m yr⁻¹ in the Gulf of Kruzenstern. As previously shown, we suggest that this effect is caused by the presence of massive ice beds, which becomes especially important given that temperature variations provide 85% of the coastal erosion (Table 3) compared to 51–53% of it for open-sea sites. Because of this, temperature variations are especially important for the temporal evolution of erosion rates of the Gulf. The recent increase in TE in 2012–2017 (Figure 9) is a result of temperature rise, while the maximum of 2005–2009 seen at other sites (Figure 9 b, c, d) is mainly caused by enhanced wind-wave energy and is not seen in the curve of TE for the Gulf of Kruzenstern (Figure 9A). The high erosion rates between 2010 and 2019 can be mainly attributed to the recent increase in the air thawing index. Similar patterns were observed in the last decade at Kharasavey, where shoreline retreat accelerated between 2006 and 2016 (Belova et al., 2017), on the Yamal coast of the Baydaratskaya Bay, where the shoreline retreated three times faster than earlier in the last decade (Novikova et al., 2018), and at Marre-Sale station, where erosion rates increased by roughly three times in 2010–2013 (Kritsuk et al., 2014). On the contrary, the Ural coast of the Baydaratskaya Bay experienced a decrease in coastal erosion in 2012–2016 (Novikova et al., 2018).

The above mentioned increase in the spatial variability of erosion rates in 2010–2019 also partly results from the high role of thermal denudation driven by the air thawing index increase. Segments of the Gulf of Kruzenstern which contain massive ice beds but receive less wave energy started to retreat faster as a result of active thawing due to air temperature increases. Examples of these process are the segments south of Cape Yasalya and south of Mordiyakha River, which were relatively stable in 1964–2010 and retreated up to 21.9 ± 6.3 m in 2010–2019. With ongoing climate and permafrost warming (Biskaborn et al., 2019), it is reasonable to expect a further acceleration in coastal erosion along segments with high massive ice contents.

Conclusion

Multitemporal imagery analysis of a coastal segment of an enclosed bay (Gulf of Kruzenstern, Western Yamal, Kara Sea) showed erosion of coastal bluffs composed by ice-rich permafrost. The average rate of shoreline retreat from 1964 to 2010 was 0.5 ± 0.2 m yr⁻¹. The fastest erosion of up to 1.7 ± 0.2 m yr⁻¹ was noted for segments with abundant outcrops of massive ice beds and numerous thermocirques. In 2010–2019, coastal erosion accelerated, reaching 0.9 ± 0.7 m yr⁻¹ on average with peaks up to 2.4 ± 0.7 m yr⁻¹. This recent increase in erosion rates results from the increase in the air thawing index, especially noticeable in 2010–2019, along with increases in sea ice-free period duration and wind-wave energy, which together contribute to an intensification of hydrometeorological forcing in the recent decade. Compared to other sites of coastal dynamics monitoring at open-sea sites in the region (Marre-Sale, Kharasavey), the Gulf of Kruzenstern experiences 30% less hydrometeorological forcing. However, due to abundance of massive ground-ice beds exposed along the bluffs of the Gulf, its shoreline retreats at comparable rates, as the thermal factor is more important in driving erosion than the wind-wave energy factor, providing 85% of the variability in shoreline retreat. Based on these findings, projections of ongoing air temperature warming in the Arctic will likely result in a continuation of rapid coastal erosion in the Gulf of Kruzenstern, and generally along wave-sheltered shorelines with ice-rich permafrost. Such erosion is substantially less dependent on changes in wave dynamics: an increase in surface air temperature alone leads to higher coastal erosion rates in ice-rich permafrost bluffs in enclosed bays and lagoons in the Arctic.

Data Availability Statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

Author Contributions

Conceptualization, AB; methodology by AB, AN, and NS; validation by AB, AN, and NS; formal analysis by AN, BJ, and AB; investigation by AB, AN, and BJ; resources by AB; data curation by AN; writing—original draft preparation by AB; writing - review and editing by AB, NB, SM and BJ; visualization by AN, AB, and NS; supervision by SO; project administration by AB; funding acquisition by SO.

Funding

Studies were funded by the RSF Grant 16-17-00034. BJ was supported by the US National Science Foundation award OISE-1927553.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

We thank Digital Globe Foundation for the WorldView-1 satellite imagery. Additional geospatial support for this work provided by the Polar Geospatial Center under NSF-OPP awards 1043681 and 1559691. Investigations were made using the equipment obtained within the State Budget Theme АААА-А16-116032810055-0. We thank the anonymous reviewers and the editors for comments and corrections which lead to substantial changes in the manuscript.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/feart.2020.566227/full#supplementary-material.

References

Aleksandrov, Y., Bryazgin, N., Førland, E., Radionov, V., and Svyashchennikov, P. (2005). Seasonal, interannual and long-term variability of precipitation and snow depth in the region of the Barents and Kara seas. Polar Res. 24. doi:10.3402/polar.v24i1.6254

CrossRef Full Text | Google Scholar

Andersland, O. B., and Ladanyi, B. (2004). Frozen ground engineering. 2nd Edn. Hoboken, NJ: John Wiley & Sons.