Rapid and synchronous response of outlet glaciers to ocean warming on the Barents Sea coast, Novaya Zemlya

Abstract

The Arctic is a hotspot for climate warming, making it crucial to quantify the sea level rise contribution from its ice masses. Novaya Zemlya's ice caps are the largest glacier complex in Europe and are a major contributor to contemporary sea level rise. Here we show that Novaya Zemlya outlet glaciers on the Barents Sea coast respond rapidly and consistently to oceanic forcing at annual timescales, likely due to their exposure to Atlantic Water variability. Glaciers on the Kara Sea show more variable response, likely reflecting their reduced exposure to Atlantic Water. Data demonstrate that the pause in glacier retreat previously observed on Novaya Zemlya between 2013 and 2015 has not persisted and that these changes correspond to ocean temperature variability on the Barents Sea coast. We document a marked shift to warmer air and ocean temperatures, and reduced sea ice concentrations from 2005 onwards. Although we identify ocean warming as the primary trigger for glacier retreat, we suggest that multi-year thinning, driven by the shift towards warmer air temperatures since 2005, pre-conditioned Novaya Zemlya's glaciers to retreat. Despite commonality in the timing of outlet glacier retreat, the magnitude is highly variable during rapid retreat phases, which we attribute to glacier-specific factors.

1. Introduction

The Arctic has warmed at twice the global average over the past 50 years, meaning that its glaciers and ice caps are undergoing rapid change and they will contribute an estimated 19–25 cm to sea level rise by 2100 (SWIPA, 2017). The Russian High Arctic (RHA) has a glaciated area of ~52 000 km² (Millan and others, 2022) and lost over 1000 Gt of ice between 1961 and 2016, equating to ~3 mm of sea level rise (Zemp and others, 2019). The RHA accounts for approximately 11% of glacial ice globally, excluding the ice sheets (Dowdeswell and others, 1997; Radić and others, 2014; Millan and others, 2022), making it a major ice reservoir and key area for studying the impacts of climate change on the Arctic cryosphere. Globally, the RHA is a major source of glacial loss since, contributing 10.6 ± 1.7 Gt a⁻¹ to sea level between 2002 and 2016 (Wouters and others, 2019). During the 21st century, the RHA is forecast to lose between 21 ± 11% (RCP2.6) and 43 ± 11% of its mass (RCP8.5), relative to 2015 (Rounce and others, 2023), and is projected to contribute between 14 ± 10 (Rounce and others, 2023) and 26 ± 10 mm (Marzeion and others, 2020) to sea level rise by 2010, under the RCP8.5 scenario. Furthermore, the Barents Sea has been identified as a key area of oceanic warming (Lind and others, 2018), sea ice decline and air temperature increases (Tepes and others, 2021b), with January 2016 temperature anomalies reaching 7 °C above the 1981–2010 mean (SWIPA, 2017). Thus, it is vital to quantify the contribution of the RHA to contemporary and near-future global sea level rise and determine the processes driving the response of its ice masses to observed rapid warming.

The RHA is composed of three main ice masses: Novaya Zemlya, which is the largest in terms of glaciated area at 21 200 km² and is the dominant source of recent ice loss in the region (Moholdt and others, 2012; Melkonian and others, 2016; Tepes and others, 2021b; Jakob and Gourmelen, 2023); Severnaya Zemlya, which has an area of 16 700 km² and has recently exhibited major, rapid-onset changes in the dynamics of some of its outlet glaciers (Willis and others, 2015, 2018; Sánchez-Gámez and others, 2019); and Franz Josef Land with an area 12 700 km², where ice losses have recently accelerated (Zheng and others, 2018). Novaya Zemlya's mass balance was negative during the 2000s (Moholdt and others, 2012; Matsuo and others, 2013; Zhao and others, 2014), with losses increasing from 10 ± 5 Gt a⁻¹ for 2003–2009 to 14 ± 4 Gt a⁻¹ for 2010–2016, and mass balance being close to zero between 2009 and 2011 (Ciracì and others, 2018). Negative mass balance coincided with retreat of Novaya Zemlya's marine-terminating glaciers: retreat rates were significantly higher between 2000 and 2013 than any other time since the mid-1970s (Carr and others, 2014, 2017a). This period of accelerated retreat coincided with atmospheric warming and reduced sea ice concentrations and retreat rates were an order of magnitude higher on glaciers terminating in the ocean than those on land (Carr and others, 2014, 2017a). Retreat rates were also significantly higher on the Barents Sea coast compared to the Kara Sea, which likely reflects their differing atmospheric and oceanic conditions (Carr and others, 2014, 2017a): the Barents Sea coast is exposed to warm ocean currents and air masses from the north Atlantic, resulting in warmer ocean and air temperatures and higher precipitation rates (Zeeberg and others, 2001; Politova and others, 2012; Przybylak and others, 2016), whereas the Kara Sea coast is exposed to cool ocean waters and air masses emanating from the Arctic Ocean (Pavlov and others, 1995; Zeeberg and others, 2001). Subsequent to enhanced retreat between 2000 and 2013, many of Novaya Zemlya's outlet glaciers underwent a short period of reduced retreat between 2013 and 2015, with several glaciers on the Barents Sea coast advancing (Carr and others, 2017a). However, it was unclear whether this was a temporary hiatus or the start of a longer-term trend of reduced retreat (Carr and others, 2017a). This is particularly timely, with the ‘Atlantification’ of the Barents Sea (Tepes and others, 2021a, 2021b; Jakob and Gourmelen, 2023), whereby air and ocean temperatures are warming and sea ice is declining (Årthun and others, 2012; Polyakov and others, 2017; Barton and others, 2018; Stroeve and others, 2018).

Seasonal-scale variations in terminus position, and the factors controlling them, have not previously been analysed on Novaya Zemlya. However, work from the Greenland Ice Sheet has shown that analysing seasonal controls can provide important insights into the factors driving interannual glacier behaviour (e.g. Carr and others, 2013; Moon and others, 2014; Fried and others, 2018; Catania and others, 2020). On the Greenland Ice Sheet, previous work has demonstrated that the onset of glacier retreat coincides with the beginning of the melt season in late spring / early summer (e.g. Sohn and others, 1998; Moon and others, 2014; Fried and others, 2018; Vijay and others, 2019; Davison and others, 2020), which may drive retreat through a combination of hydrofracturing and/or increased basal melting, due to enhanced plume flow (e.g. Vieli and others, 2011; Carr and others, 2013; Straneo and others, 2013). Similarly, ocean temperatures display seasonal cycles, and can therefore seasonally-enhance basal melting and undercutting of the terminus in the warmer months (e.g. Benn and others, 2007; Straneo and others, 2012; O'Leary and others, 2013; Straneo and others, 2013; Todd and others, 2019). Finally, an ice mélange (a semi-rigid mix of icebergs and sea ice) forms seasonally at ice fronts in Novaya Zemlya and the Greenland Ice Sheet and can strongly suppress the movement of icebergs away from the terminus (e.g. Sohn and others, 1998; Amundson and others, 2010; Carr and others, 2013; Moon and others, 2015). Thus, an earlier onset and/or late formation of the ice mélange can result in higher calving rates and/or a longer calving season and hence net recession of the terminus (Robel, 2017; Amundson and others, 2018; Todd and others, 2019). Previous work on Novaya Zemlya has suggested that glacier retreat corresponds with warmer air temperatures and reduced sea ice concentrations at decadal timescales, and that ocean warming may play a role (Carr and others, 2017a), but the role of each forcing factor remains uncertain, particularly at seasonal timescales. To address this knowledge gap, we construct the first time series of seasonal frontal position variation on Novaya Zemlya for a subset of our study glaciers and analyse these variations in relation to seasonal-scale climatic forcing factors.

Changes in the dynamics of marine-terminating outlet glaciers represent an important ice loss mechanism from Arctic ice caps and ice sheets (e.g. Mouginot and others, 2019; IMBIE, 2020; Kochtitzky and others, 2022), both through immediate losses from iceberg calving and longer-term draw-down of inland ice, or ‘dynamic thinning’ (e.g. Pritchard and others, 2009; Price and others, 2011). Across the Barents and Kara Sea region, changes in ice discharge accounted for 43% of ice mass loss between 2010 and 2020, with the remaining 57% resulting from negative surface mass balance (Jakob and Gourmelen, 2023). On Novaya Zemlya, SMB accounts for the majority of mass loss (89%), while increased discharge account for 76 and 100% of mass loss for Franz Josef Land and Severnaya Zemlya respectively (Jakob and Gourmelen, 2023). Rapid thinning rates have been observed on Novaya Zemlya, particularly at low elevations (Tepes and others, 2021a, 2021b), but it is unclear whether observed thinning on Novaya Zemlya is a response to enhanced surface melt, which has then promoted glacier retreat and terminus ungrounding, or whether outlet glacier retreat has caused thinning, due to draw down of inland ice. Thus, this paper aims to systematically investigate the relationship between terminus retreat and surface elevation change on Novaya Zemlya.

Here, we use remotely sensed data to: (i) determine whether reduced retreat rates observed on Novaya Zemlya between 2013 and 2015 have persisted or whether the enhanced retreat observed between 2000 and 2012 has resumed; (ii) quantify seasonal-scale variations in frontal positions for selected glaciers; (iii) assess the relationship between frontal position change and external forcing at seasonal to interannual timescales and; (iv) determine whether glacier retreat corresponds to thinning and the direction of this relationship, i.e. does retreat precede or follow thinning. We determined frontal position changes between 2015 and 2020 for all major glaciers (i.e. >1 km wide) identified by Carr and others, 2017a. To investigate seasonal-scale variations in frontal position and their relationship with forcing, we also mapped terminus positions for six glaciers on each coast, at weekly to monthly temporal resolution between 2010 and 2020. To do this, we selected glaciers that had the biggest change in retreat rates between 2013–2015 and 2015–2020, thus enabling us to assess the switch in retreat behaviour in more detail. Surge-type glaciers were excluded (Grant and others, 2009), leaving a total of 53 outlets (Fig. 1), which were distributed across the three main ice caps of Novaya Zemlya (northern island ice cap, subsidiary ice mass 1 and subsidiary ice mass 2), and between the Barents and the Kara Sea coasts. Retreat rates were analysed in relation to: (1) air temperature, determined from ERA-5 reanalysis data and meteorological station data; (2) sea ice and (3) ocean potential temperature data from the Global Ocean Physics Reanalysis dataset GLORYS12V1; and (4) surface elevation change from CryoSat-2 interferometric radar altimetry.

Figure 1.

Location map showing the location of (a) Novaya Zemlya within the Russian High Arctic and (b) to location of the main ice masses and study glaciers on Novaya Zemlya. Glaciers are symbolised according to terminus type and surge-type glaciers are identified with a star. Three letter IDs are given for each marine-terminating glacier and a full list of glacier IDs and names are given in Supplementary Table 1.

[Figure omitted. See PDF]

2. Methods

2.1. Glacier frontal positions

Outlet glacier frontal positions were determined for the 53 study glaciers between 2015 and 2020, which were >1 km in width (Fig. 1). The outlets are located on Novaya Zemlya's main ice mass (Fig. 1), the northern ice cap, also known as Severny Island ice cap (n = 45), and its two smaller ice masses, Sub 1 (n = 4) and Sub 2 (n = 4). We include glaciers of different terminus types, specifically, ocean-terminating (n = 32), land-terminating (n = 15) and lake terminating (n = 6). On the northern ice cap, outlet glaciers flow from the main divide, either west towards the Barents Sea coast (n = 27) or east towards the Kara Sea coast (n = 18; Fig. 1). Glaciers on Sub 1 and Sub 2 were not categorised as Barents or Kara Sea, as they predominantly flow north/south from their ice divides (Fig. 1).

We mapped terminus position for all 53 study glaciers between 2015 and 2020 from Landsat 8 imagery, which was acquired from the USGS Earth Explorer and has a spatial resolution of 30 m. To avoid seasonal variability and ensure comparability with data from Carr and others (2017a), imagery was obtained from as close to 31 July as possible for each year. Following previous studies, glacier termini were digitised for each time interval using the box method, whereby a reference box was created for each glacier, which falls between its bounding fjord walls. The terminus was then digitised from each image, so that any changes in area are due to changes in frontal position, which was calculated by dividing area change by the reference box width. This enables uneven retreat across the terminus to be captured. To investigate seasonal scale terminus variations and the period of reduced retreat rates between 2013 and 2015, we produced a high temporal resolution time series of frontal positions for six glaciers on both the Barents Sea (CHA, INO, KRI, NOR, SH, VIJ, VIS) and the Kara Sea coasts (KRO, OGA, POL, SHU1, VER, VYL2). To construct the most comprehensive time series available, we used Imagery from Landsat 7 and 8 and Sentinel 1 and 2. Data availability ranged from a few images per year between 2010 and 2015 to approximately weekly resolution from 2016 onwards, when Landsat 8 and Sentinel data were available. Average retreat rates for the 12 glaciers were 90.2 m a⁻¹ between 2015 and 2020, which equates to 1.7 m per week, meaning that weekly to monthly resolution data are likely sufficient to capture the majority of short-term terminus change.

All imagery used to map terminus positions was assessed for potential geolocation errors, by comparing features that should not move between images (e.g. bedrock ridges), and no discernible differences in geolocation were observed. Thus, the primary source of error in our terminus positions comes from digitising errors associated with correct identification of the glacier margin. To assess these errors, we re-digitised the terminus positions for three marine-terminating, three lake terminating and three land-terminating glaciers for each time interval between 2015 and 2020, without reference to the original mapping. We then calculated the average divergence between the frontal positions for each image date, for each glacier. The resultant errors for marine-terminating glaciers was 26.8, 35.4 m for lake-terminating and 77.6 m for land-terminating glaciers, equating to errors in retreat rate of 5.4, 7.1 and 15.5 m a⁻¹ for the period 2015 to 2020 respectively.

Frontal positions from 2015 to 2020 extend the previous dataset created by Carr and others (2017a), which documented frontal position changes on Novaya Zemlya between 1973/6 and 2015. Here we use the same study glaciers and reference boxes, to allow for direct comparison and to construct a coherent, consistent record through to 2020. The Carr and others, 2017a dataset was created from a range of satellite imagery, including: Landsat 1–5, Landsat 7 and Landsat 8 imagery from the 1970s through to the 2010s; declassified Hexagon HK-9 imagery from 1976 and 1977; and Synthetic Aperture Radar (SAR) Image Mode Precision data from the European Remote-sensing Satellite-1 (ERS-1) and ERS-2 for the 1990s. Full details of the specific imagery used are available in Carr and others, 2017a. Frontal position errors for the 1973/6–2015 datasets were calculated at 17.5 m for marine- and lake-terminating glaciers, equating to an error in retreat rates of 3.5 m a⁻¹ per five-year period, which is comparable to our errors of 5.4 and 7.1 m a⁻¹ for 2015 to 2020. For land-terminating glaciers, the error in the Carr and others (2017a) dataset was 68.4 m, equating to a five-year error in retreat rates of 13.7 m a⁻¹, which is also comparable to the error in retreat rate of 15.5 m a⁻¹ for the 2015–2020 data.

2.2. Surface elevation change data

Surface elevation data were sampled from CryoSat-2 interferometric radar altimetry data for the period August 2010 to August 2020. Data processing follows the methods outlined in Jakob and Gourmelen (2023). Swath processing was applied to the CryoSat-2 data, which increases the quantity of data retrieved, compared to conventional processing approaches, and increases the temporal and spatial resolution across complex topography (Foresta and others, 2016; Gourmelen and others, 2018; Jakob and others, 2021; Jakob and Gourmelen, 2023). In this study we use two surface elevation change datasets; (1) a gridded product of 10-year elevation change rates and (2) a glacier-wide monthly time series of elevation changes. For the 10-year surface elevation change dataset a plane-fit algorithm was applied to the swath points to derive a gridded elevation change raster at a spatial resolution of 500 m (Jakob and Gourmelen, 2023). To obtain the monthly elevation change time series the swath data are gridded into 2 km pixels using a 90-day moving time window (Jakob and Gourmelen, 2023). The RGI v6.0 outlines are used to generate integrated elevation change time series over the study glaciers. We sampled elevation change from within a 5 km buffer of the July 2010 terminus positions for each study glacier. Buffers were then clipped to the RGI v6.0 outlines, so that the elevation change was sampled for the area within each outlet glacier basin with 5 km of the July 2010 terminus and within the glacier catchment, as defined in the RGI v6.0.

2.3. Climatic and oceanic data

Data on air temperatures were obtained from the ERA5 reanalysis dataset ‘ERA5 monthly averaged data on single levels from 1979 to present’ (Hersbach and others, 2019), which have a spatial resolution of 0.25° and cover the period 1979 to present. The 2 m temperature product was used, and data were sampled from the gird cell closest to Novaya Zemlya's coast, using the product's ocean vs land mask. Average values were calculated for the Barents and Kara Sea coasts. For both the meteorological station data and the reanalysis data, we calculated meteorological seasonal means for Winter (December-February), Spring (March-May); Summer (June-August); and Autumn (September-November). Seasonal and annual means were only calculated if data were available for all months.

Sub-surface ocean temperatures and sea ice data were obtained from the Global Ocean Physics Reanalysis dataset GLORYS12V1 (Fernandez and others, 2018), provided by the Copernicus Marine Environment Monitoring Service (CMEMS). The dataset has a horizontal resolution of 0.083° and 50 vertical levels in the ocean and extends from 1 January 1993 to 31 December 2019. We determined sea ice concentration for the Barents and Kara seas, by sampling from the grid squares immediately offshore of each coast, using the land/ocean mask. The number of ice-free months (i.e. months where the concentration was zero) were also calculated. We also sampled potential ocean temperature data from the GLORYS12V1 dataset, at depths of 5 m (i.e. near-surface water) and 92 m, with 92 m being the greatest depth attainable immediately offshore of Novaya Zemlya and thus likely to be most representative of the ocean conditions influencing glacier behaviour (Carr and others, 2017a).

2.4. Statistical analysis

We used a Wilcoxon test to identify significant differences between our study time periods, for glacier retreat rates, air temperatures, sea ice concentrations, number of sea ice free months and ocean temperatures at depths of 5 and 92 m. Wilcoxon tests were used, as the test is nonparametric and is suitable for assessing statistical differences between two time periods of data (Miles and others, 2013; Carr and others, 2017b). Following statistical convention, a p-value of less than or equal to 0.05 is taken as significant. The Matlab function xcorr was used to calculate the cross-correlation between the time series of surface elevation change and relative frontal position change for each glacier between 2011 and 2020. Cross-correlation determines the similarity between one discrete time series (x) and another discrete time series (y), at a series of lagged intervals, i.e. it calculates the correlation between the two time series for x _{time =1} and y _{time =1}, then x _{time =1} and y _{time =2}, then x _{time =1} and y _{time =3}, until the final value of y is reached in the time series. In our analysis, x was the surface elevation time series and y was the relative frontal position time series, so that:

• (i) A lag of −1 years evaluates the correlation between glacier frontal position (y) and surface elevation change in the previous year (x), i.e. glacier frontal position (y) at time 0, is correlated with surface elevation change (x) at time −1.

• (ii) A lag of + 1 years evaluates the correlation between surface elevation change (x) and glacier frontal position (y) in the previous year, i.e. surface elevation change (x) at time 0, is correlated with glacier frontal position (y) at time −1.

• (iii) A lag of 0 evaluates the correlation between surface elevation change and relative frontal position for the same year, i.e. surface elevation change (x) at time 0, is correlated with glacier frontal position (y) at time 0.

The correlation coefficient for each cross-correlation was normalised to give a value between 0 and 1, with 1 indicating the highest correlation. This is to ensure comparability between different glaciers. As such, we use the cross correlation to determine the timing of the maximum correlation and to thus identify when the maximum correlation occurs between surface elevation change and glacier frontal position.

3. Results

3.1. Interannual changes in glacier frontal position

Results show that 92% (49 out of 53 glaciers) of our study glaciers retreated between 2015 and 2020, across all terminus types (Figs 1, 2). The highest retreat rates occurred on marine-terminating outlet glaciers located on the Barents Sea coast, where rates averaged −116.9 m a⁻¹ and reached a maximum of −297.0 m a⁻¹ at KRI (Fig. 2). Retreat rates on these glaciers were significantly higher than during the 2013–2015 pause, when termini underwent a slight net advance on average (+11.4 m a⁻¹) and were comparable to those observed during the previous phase of marine-terminating outlet glacier retreat on Novaya Zemlya between 2000 and 2013 (−107.6 m a⁻¹; Table 1; Fig. 2). Mean retreat rates from 2015–2020 on the Barents Sea coast were approximately five times greater than those observed prior to the 2000s (Table 1; Fig. 2). The range in rates of frontal position change between individual marine-terminating outlet glaciers on the Barents Sea Coast was much higher during periods of rapid retreat, reaching 299.4 m a¹ (2000–2013) and 294.2 m a¹ (2015–2020), compared to 160.5 m a¹ (1976–1986), 73.5 m a¹ (1986–2000) and 108.5 m a¹ (2013–2015).

Figure 2.

Box chart showing median retreat rates of Novaya Zemlya outlet glaciers, divided by study time period and terminus type. Retreat rates are plotted for the five study time periods: 1973/76–1986, 1986–2000, 2000–2013 and 2013–2015 and 2015–2020. Data are split according to terminus type (marine, land and lake) and marine-terminating glaciers are further sub-divided by coast. The median retreat rate is shown as the line inside the box and the upper (lower) quartile is shown as the top (bottom) edges of the box. The lines extending beyond the boxes show the full range in values and any outliers are indicated by circles.

[Figure omitted. See PDF]

Table 1. P-values from Wilcoxon test results, to test for statistical difference between the study time periods (1973/76–1986, 1986–2000, 2000–2013 and 2013–2015 and 2015–2020) for glacier retreat

Marine – Kara	76–86	86–00	00–13	13–15	15–20
76–86		0.3408	0.0404	0.2366	0.0194
86–00			0.0120	0.9770	0.0102
00–13				0.0061	0.8399
13–15					0.0056
15–20
Marine – Barents	76–86	86–00	00–13	13–15	15–20
76–86		0.4233	0.0000	0.0056	0.0000
86–00			0.0007	0.0008	0.0000
00–13				0.0000	0.7360
13–15					0.0000
15–20
Land	76–86	86–00	00–13	13–15	15–20
76–86		0.4763	0.1029	0.6961	0.0063
86–00			0.6295	0.2603	0.0769
00–13				0.1826	0.1982
13–15					0.0660
15–20
Lake	76–86	86–00	00–13	13–15	15–20
76–86		0.2403	0.3939	0.4848	0.1017
86–00			0.7316	0.5887	0.5887
00–13				0.5887	0.2403
13–15					0.3095
15–20

Data were divided by terminus type (marine-, land -and lake-terminating glaciers) and the marine-terminating glacier category was further sub divided by coast (Barents and Kara). Significant p-values <0.05 are highlighted in bold and p-values are given to four decimal places.

Marine-terminating outlets on the Kara Sea coast followed a similar, although more muted, pattern to those on the Barents Sea: retreat rates were highest in 2015–2020, at −64.6 m a⁻¹, and this was a significant increase from 2013–2015, when retreat rates averaged −18.2 m a⁻¹ (Table 1; Fig. 2). As on the Barents Sea coast, retreat rates between 2015 and 2020 were not significantly different from those between 2000 and 2013 (−58.4 m a⁻¹) but were three times higher than between 1973/76 and 2000 (Table 1; Fig. 2). Variability in retreat rates between glaciers was also much higher during periods of rapid retreat, with a range of 172.7 m a¹ for 2000–2013 and 208.0 m a¹ for 2015–2020, compared to 32.6 m a¹ (1976–1986), 64.7 m a¹ (1986–2000) and 69.3 m a¹ (2013–2015).

In contrast to their marine-terminating neighbours, the lake-terminating outlets showed no significant differences in their retreat rates between any of the study time periods: we did not observe significantly higher retreat rates in either 2000–2013 or 2015–2020 nor significantly reduced retreat in 2013–2015 (Fig. 2, Fig. S1; Table 1). Indeed, mean retreat rates were notably consistent between time periods, ranging from −30.6 m a¹ in 1976–1986 to −61.7 m a¹ in 2015–2020 (Fig. 2), and all lake-terminating outlets followed a linear pattern of terminus retreat between 1976 and 2020 (Fig. S1). The range in retreat rates was highest in 2013–2015 (90.8 m a⁻¹) and lowest in 2015–2020 (66.9 m a⁻¹), meaning that the greatest range in retreat rates did not correspond to the periods of highest retreat rates. Land-terminating outlet glaciers experienced the lowest retreat rates of any terminus type, ranging between −7.1 m a⁻¹ in 1976–1986 and −32.1 m a⁻¹ between 2015 and 2020 (Fig. 2, Fig. S1). Retreat rates were higher in 2015–2020 (−32.1 m a⁻¹) than in any other time period, but this difference was not statistically significant, and no significant differences were found between retreat rates for any of the time periods (Table 1). As for lake-terminating glaciers, variability in retreat rates for land-terminating outlet glaciers was highest for 2013–2015 (−137.9 m a⁻¹) and showed no clear correspondence to mean retreat rates (Fig. 2).

3.2. Seasonal cycles of frontal position change

Across the six study glaciers on the Barents Sea coast, seasonal retreat began between mid-June and mid-July each year and most often started in July, for the period 2016 to 2020 (Figs 3a, 4a; Table S2). The most seasonally retreated position showed more temporal and spatial variability and generally occurred between November and December (Figs 3b, 4a; Table S2). On all six of the glaciers, larger than average seasonal retreats occurred in both 2016 and 2020 (Fig. 4a). The duration of seasonal retreat was not greater than usual during these events, meaning that the additional retreat occurred through higher rates of retreat (Fig. 4a). For example, KRI retreated 832 m between 6 July and 27 December 2016, compared to just 182 m between 7 July and 22 December 2017 (Fig. 4a). Retreat rates in 2016 were between 5 (VIJ) and 10 (VIS) times higher than the average for each glacier between 2015 and 2020.

Figure 3.

Histogram showing the month of the year for seasonal changes in frontal position for selected glaciers on the Barents Sea coast (a & b) and Kara Sea coast (c & d). Frequency is shown for the date of seasonal retreat onset (a & c) and for the date of the end of seasonal retreat (b & d).

[Figure omitted. See PDF]

Figure 4.

High temporal resolution frontal position of selected study glaciers on the Barents Sea coast, shown with anomalies in forcing factors. (a) Frontal positions are relative to 7 July 2010 and are colour-coded by glacier name. Glaciers are in order of latitude, with the first glacier (CHA) being furthest north and the last glacier (NOR) being furthest south. Glacier locations are shown in Figure 1. (b) Air temperature anomalies (°C) from ERA5 data, for the Barents Sea coast. Anomalies were calculated relative to the 1979–2020 mean (i.e. the duration of the ERA5 data), for all seasons and annual values. (c) Sea ice concentration anomalies (%) for the Barents Sea coast from the Global Ocean Physics Reanalysis dataset GLORYS12V1. Anomalies were calculated relative to the 1993–2020 mean (i.e. the duration of the GLORYS12V1 data), for all seasons and annual values. (d & e) Ocean temperature anomalies (°C) for the Barents Sea coast GLORYS12V1, at depths of 5 (d) and 92 m (e). Anomalies were calculated relative to the 1993–2020 mean and are shown for all seasons and annual values.

[Figure omitted. See PDF]

On the Kara Sea coast, seasonal retreat usually began in June or July (Figs 3c, 5a; Table S3). The glaciers usually reached their seasonally retreated position in late October, with mid-October and early November also being common (Figs 3d, 5a; Table S3). Changes in the magnitude of seasonal retreat varied between the six study glaciers. In 2018, a larger seasonal retreat was observed on KRO (823 m) and on POL (552 m), which, in both cases, was driven by greater retreat rates, rather than a longer season (Figs 5a, b). POL also retreated substantially in 2020 (465 m). The other four study glaciers on the Kara Sea coast underwent larger retreats in 2019 and 2020, and ROZE only established a clear seasonal cycle in these years (Figs a, b).

Figure 5.

High temporal resolution frontal position of selected study glacier on the Kara Sea coast, shown with anomalies in forcing factors. (a) Frontal positions are relative to 7th July 2010 and are colour-coded by glacier name. Glaciers are in order of latitude, with the first glacier (ROZE) being furthest north and the last glacier (POL) being furthest south. Glacier locations are shown in Figure 1. (b) Air temperature anomalies (°C) from ERA5 data, for the Kara Sea coast. Anomalies were calculated relative to the 1979–2020 mean (i.e. the duration of the ERA5 data), for all seasons and annual values. (c) Sea ice concentration anomalies (%) for the Kara Sea coast from the Global Ocean Physics Reanalysis dataset GLORYS12V1. Anomalies were calculated relative to the 1993–2020 mean (i.e. the duration of the GLORYS12V1 data), for all seasons and annual values. (d & e) Ocean temperature anomalies (°C) for the Kara Sea coast GLORYS12V1, at depths of 5 (d) and 92 m (e). Anomalies were calculated relative to the 1993–2020 mean and are shown for all seasons and annual values.

[Figure omitted. See PDF]

3.3. Climate and ocean forcing

Between 1990 and 2020, air temperatures from ERA5 reanalysis data followed a similar pattern on the Barents and Kara Sea coasts and showed large interannual variability (Figs 6b, e). Overall, the record was dominated by positive air temperature anomalies from 2005 onwards with negative anomalies beforehand (Figs 6c, f). This is supported by Wilcoxon test results, which showed that both the Barents and Kara Sea coasts were significantly warmer in 2015–2020 than in 1986–1999, for all seasons and annually (Table 2). The air temperature record on both coasts was characterised by 1 to 3 year long periods of warmer or cooler temperatures: warmer air temperatures were observed in 2011, 2012, 2016 and 2020, while cooler temperatures occurred in 1998, 2004, 2009 and 2013 to 2015 (Figs 4b, 5b, 6b, e). However, there was no significant difference between air temperatures in 2015–2020 and 2013–2014 during any season, while summer and annual mean temperatures in 2015–2020 were significantly warmer than those in 2000–2012 (Table 2).

Figure 6.

Relative frontal position change of marine-terminating outlet glaciers and air temperature (T) between 1990 and 2020, for the Barents Sea coast (left-hand panels) and the Kara Sea Coast (right-hand panels). Frontal positions (a and d) are relative to 1986, which is the earliest date common to all glaciers, and are colour-coded according to glacier name, in alphabetical order. Glacier locations are shown in Figure 1. Air temperatures are from ERA5 data and are for the Barents Sea coast (b) and Kara Sea coast (e), for the period 1990 to 2020. Data are coloured coded by annual and meteorological season averages. Air temperature anomalies were calculated for the Barents Sea coast (c) and Kara Sea coast (f) for the annual and seasonal means and are relative to the 1979–2020 mean (i.e. the duration of the ERA5 data).

[Figure omitted. See PDF]

Table 2. P-values from Wilcoxon test results, to test for statistical difference between 2015–2020 and the other study time periods (1979–1985, 1986–1999, 2000–2012 and 2013–2014) for air temperatures from ERA5 reanalysis data

ERA5 – Barents	79–85	86–99	00–12	13–14
Spring	0.0111	0.0335	0.0813	0.6667
Summer	0.5350	0.0335	0.0265	0.6667
Autumn	0.1282	0.0041	0.0684	0.2222
Winter	0.1807	0.0009	0.2346	0.8889
Annual	0.0175	0.0004	0.0394	0.2222
ERA5 – Kara	79–85	86–99	00–12	13–14
Spring	0.0111	0.0277	0.1130	0.6667
Summer	0.3176	0.0277	0.0394	0.6667
Autumn	0.0973	0.0012	0.0684	0.1111
Winter	0.1375	0.0015	0.2048	0.6667
Annual	0.0111	0.0003	0.0394	0.2222

Significant p-values <0.05 are highlighted in bold and p-values are given to four decimal places.

On both the Barents and Kara Sea coasts, sea ice concentrations were lower than the 1993–2019 mean across all seasons and annually from 2005 onwards (Figs 7c, g). There was also large interannual variability in sea ice concentrations along both coasts, particularly during winter and spring (Figs 7b, f). The number of ice-free months on both coasts increased over time: on the Barents Sea (Kara Sea) coast, there were an average of 5 (2.5) ice free months between 1993 and 2005, compared to 6.4 (4.4) months between 2006 and 2019 (Figs 7d, h). The Barents Sea coast was sea ice free during the summer for the entire study period (Fig. 7d) and the Kara Sea became ice free in summer from 2005 onwards (Fig. 7h). Wilcoxon tests showed no significant difference in sea ice concentrations or the number of ice-free months in 2015–2019, compared to either 2000–2012 or 2013–2014, with the exception of autumn on the Barents Sea coast (Table 3). As with air temperatures, sea ice concentrations were characterised by short-term (1–3 years) periods of lower or higher concentrations. On both coasts, sea ice concentrations were comparatively high in 1998 and low in 2000–2001 (Figs 7b, f). Sea ice concentrations then generally declined over the 2000s, but with large interannual variability, before reaching the lowest concentrations observed during the study period in 2012 (Figs 4b, 5b, 7b, f). Subsequently, sea ice concentrations on both coasts increased between 2013 and 2015 (Figs 7b, f) and reached another low point in 2016 (Figs 4b, 5b, 7b, f).

Figure 7.

Relative frontal position change of marine-terminating outlet glaciers and sea ice conditions between 1990 and 2020, for the Barents Sea coast (left-hand panels) and the Kara Sea Coast (right-hand panels). Frontal positions (a and d) are relative to 1986, which is the earliest date common to all glaciers, and are colour-coded according to glacier name, in alphabetical order. Glacier locations are shown in Figure 1. Sea ice concentrations (%) are from the Global Ocean Physics Reanalysis dataset GLORYS12V1 and are shown for the Barents Sea coast (b) and Kara Sea coast (e), for the period 1993 to 2019. Data are coloured coded by annual and meteorological season averages. Sea ice concentration anomalies (%) were calculated for the Barents Sea coast (c) and Kara Sea coast (f) for the annual and seasonal means and are relative to the 1993–2019 mean (i.e. the duration of the GLORYS12V1 data). The average number of sea ice free months for each year is displayed for the Barents Sea (d) and Kara Sea (h).

[Figure omitted. See PDF]

Table 3. P-values from Wilcoxon test results, to test for statistical difference between the study time periods (1993–1999, 2000–2012 and 2013–2014 and 2015–2020) for sea ice concentrations and the number of ice-free months

Sea ice conc. – Barents	93–99	00–12	13–14
Spring	0.1490	0.1433	0.5714
Summer	1	1	1
Autumn	0.0025	0.0376	0.3810
Winter	0.0101	0.1433	0.5714
Annual	0.0025	0.0754	0.5714
Sea ice conc. – Kara
Spring	0.3434	1	0.5714
Summer	0.1414	0.7010	1
Autumn	0.0480	0.2460	0.3810
Winter	0.1061	0.3873	0.5714
Annual	0.0177	0.1734	1
Ice free months
Barents	0.1010	0.2745	0.8571
Kara	0.0101	0.1902	0.8571

Data were divided by coast. Significant p-values <0.05 are highlighted in bold and p-values are given to four decimal places.

Potential ocean temperatures on both coasts and at both depths (5 m and 92 m) showed mostly positive anomalies after 2010 and negative anomalies beforehand (Figs 8c, e, h, j). There was large interannual variability, which was particularly marked on the Barents Sea coast, and this variability increased substantially on both coasts after 2010 (Figs 8b, d, g, i). For example, on the Barents Sea coast at 92 m depth, potential temperature warmed by 0.40 °C between 1993 and 1995, then cooled by 0.61 °C by 1997 (Fig. 8d). In contrast, temperatures cooled by 1.21 °C between 2012 and 2014 and then warmed by 1.35 °C between 2014 and 2016 (Fig. 8d). On the Barents Sea, ocean temperatures at 5 at 92 m depth were significantly warmer in 2015–2019 for all seasons (except winter 2000–2012) than for 1993–1999 and 2000–2012, but were not significantly warmer than ocean temperatures in 2013–2014 (Table 4). For the Kara Sea, ocean temperatures at 5 m depth were significantly warmer in 2015–2019 than for 1993–1999 for summer, spring and annual values, while at 92 m depth, 2015–2019 temperatures were warmer than those in 1993–1999 (all seasons except spring) and 2000–2012 (in summer autumn and annual data). Ocean temperatures were characterised by 1 to 4 year periods of warming and cooling, which were particularly marked on the Barents Sea and after 2010 (Figs 8b, d). On the Barents Sea, ocean temperatures at both depths were particularly warm in 2012 and 2016 in all seasons (Figs 8b, d) and the warmest anomalies were observed in 2016, particularly in summer and autumn (Figs 4d, e). At both depths, ocean temperatures were comparatively cool in 2010, 2014 and 2019 (Figs 8b, d), with 2014 showing slightly cooler than average conditions (Figs 4d, e). On the Kara Sea coast, ocean temperatures were also warmest in 2012 and 2016 (Figs 8g, i), particularly during the summer (Figs 4d, e). At both depths, ocean temperatures were cool in 214, 2017 and 2018 (Figs 4d, e, 8g, i).

Figure 8.

Relative frontal position change of marine-terminating outlet glaciers and ocean temperature (T) between 1990 and 2020, for the Barents Sea coast (left-hand panels) and the Kara Sea Coast (right-hand panels). Frontal positions (a and d) are relative to 1986, which is the earliest date common to all glaciers, and are colour-coded according to glacier name, in alphabetical order. Glacier locations are shown in Figure 1. Ocean temperatures are from the Global Ocean Physics Reanalysis dataset GLORYS12V1, for the period 1993–2019. Both absolute ocean temperatures and anomalies relative to the 1993–2019 mean (i.e. the duration of the GLORYS12V1 data) are shown for depths of 5 m and 92 m depth. Data are coloured coded by annual and meteorological season averages.

[Figure omitted. See PDF]

Table 4. P-values from Wilcoxon test results, to test for statistical difference between the study time periods (1993–1999, 2000–2012 and 2013–2014 and 2015–2020) for potential ocean temperature

Ocean temperature – Barents: 5 m	93–99	00–12	13–14
Spring	0.0101	0.0350	0.1905
Summer	0.0177	0.0264	0.8571
Autumn	0.0051	0.0350	0.3810
Winter	0.0025	0.0593	0.5714
Annual	0.0025	0.0068	0.8571
Ocean temperature – Kara: 5 m	93–99	00–12	13–14
Spring	1	0.3873	0.5714
Summer	0.0303	0.1172	0.5714
Autumn	0.0101	0.1433	0.5714
Winter	0.5303	0.9241	1
Annual	0.0480	0.1734	0.5714
Ocean temperature – Barents: 92 m	93–99	00–12	13–14
Spring	0.0303	0.0098	0.3810
Summer	0.0051	0.0068	0.1905
Autumn	0.0051	0.0350	0.1905
Winter	0.0480	0.0593	0.5714
Annual	0.0051	0.0068	0.1905
Ocean temperature – Kara: 92 m	93–99	00–12	13–14
Spring	0.2677	0.6331	0.3810
Summer	0.0101	0.0350	0.8571
Autumn	0.0025	0.0264	0.5714
Winter	0.0101	0.5028	0.8571
Annual	0.0177	0.0194	0.8571

Data were divided by coast (Barents and Kara coasts) and test statistics were calculated for depths of 5 and 92 m. Significant p-values <0.05 are highlighted in bold and p-values are given to four decimal places.

3.4. Glacier surface elevation change

Results show a significant correlation between outlet glacier retreat rates for the period 2015–2020 and surface thinning rates between 2010 and 2020 (Fig. 9). Specifically, regressing 2015–2020 retreat against surface elevation change within 5 km of the terminus gave an R² value of 0.46 for all glaciers and 0.36 for marine terminating outlets only, both with a p-value of 0.00, to two decimal places. Despite an overall significant correlation, there was a spread in the relationship between thinning and retreat rates. For example, mean retreat rates on KRA1, KRI, KRO and VIS were high compared to thinning rates, i.e. the four glaciers were substantially below the linear regression line in Figure 9. The four glaciers were all marine terminating, and three glaciers were located close together on central portion of Novaya Zemlya's Barents Sea coast (KRA1, KRI and VIS) and one was located midway along Novaya Zemlya's Kara Sea coast (KRO; Fig. 1). Conversely, four glaciers exhibited retreat rates that were low compared to their thinning rates, i.e. they sat above the regression line in Figure 9, specifically: RYK and KRA2, located on the central Barents Sea coast, INO, located on the northern portion of the Barents Sea coast, and VYL1, located on the central Kara Sea coast (Fig. 1).

Figure 9.

Mean retreat rate for 2015–2020 regressed against mean thinning rate for Aug 2010 to Aug 2020 for: (a) All glacier terminus types (i.e. marine-, land- and lake-terminating glaciers); and (b) marine terminating glaciers only. Thinning rates were sampled from within 5 km from each glacier terminus, as of July 2010, and within each glacier catchment, as defined in the Randolph Glacier Inventory (RGI) v6.0. Thinning rates were calculated across Novaya Zemlya for the period August 2010 to August 2020 and the average thinning rate for each glacier catchment was calculated, using the RGI v6.0 outlines.

[Figure omitted. See PDF]

We used cross-correlation to evaluate the relationship between time series of elevation change and frontal position change and to identify any lags and/or leads within this relationship (Fig. 10). For most marine-terminating outlet glaciers, relative frontal position showed strong correlation with the previous ~5 years of surface elevation change and the maximum correlation occurred with zero lag: I.e. the strength of the correlation (shown by the length of the stick on the y-axis in Fig. 10) was greatest when the lag value on the x axis was between 0 and −5 years (i.e. the correlation between frontal position and thinning in the same year, through to thinning five years previously). This suggests that thinning over the previous ~5 years may drive frontal position change and that thinning in the same year has the most significant influence. At some glaciers (e.g. KRA1 and VEL), the correlation between surface elevation change and frontal position change decreases rapidly as the lag increases i.e. the influence of past surface elevation change on glacier retreat reduces rapidly with time. At other outlets (e.g. GLA, INO and VYL1), the decline in correlation with increasing lag is less marked, suggesting that previous thinning impacts glacier frontal position for longer. Focusing on positive lags, i.e. the relationship between current thinning rates and past frontal position change, the correlation reduces rapidly with time, suggesting that the impacts of changes in glacier frontal position on surface elevation change are shorter-lived than the impacts of thinning on retreat. Finally, a small number of glaciers (e.g. BRO and VIJ) show a different relationship between changes in frontal position and surface elevation, whereby there is no peak in correlation at the zero-lag time and surface elevation change shows a consistent correlation with frontal position for all negative lags. Conversely, there is an inverse correlation between surface elevation change and relative frontal position after the zero lag. This suggests that frontal position change at these glaciers is approximately equally influenced by surface elevation change in the previous years, and that retreat (advance) may drive thickening (thinning), which is the inverse of the relationship observed at the other glaciers.

Figure 10.

Cross correlation results for surface elevation change between 2011 and 2020 vs relative frontal position change between 2011 and 2020. Stick height (y-axis) indicates the strength of the correlation at each lag time in years (x axis) between surface elevation change and relative frontal position change. Cross correlation evaluates correlation at all possible lags (indicated by a positive value on the x axis) and leads (indicated by a negative value on the x axis) between the two datasets. For example, the value at −1 on the x axis gives the correlation between frontal position change and thinning one year previous. Conversely, the value at 1 on the x axis gives the correlation between thinning and frontal position change one year previous. Finally, the value at x = 0 gives the correlation between frontal position change and thinning during the same year, i.e. no lag in either direction.

[Figure omitted. See PDF]

4. Discussion

4.1. Have reduced retreat rates observed on Novaya Zemlya between 2013 and 2015 persisted?

The pause in marine-terminating outlet glacier retreat observed on Novaya Zemlya between 2013 and 2015 (Carr and others, 2017a) has not persisted and, instead, rapid retreat has resumed (Fig. 2): retreat rates between 2015 and 2020 were significantly higher than those observed pre-2000 (Table 1) and were comparable to Novaya Zemlya's previous phase of rapid retreat between in 2000–2012 (Table 1; Fig. 2) (Carr and others, 2017a). Retreat rates on land- and lake-terminating glaciers on Novaya Zemlya showed no significant difference between any of the time periods (Table 1; Fig. 2; Fig. S1), suggesting that periods of enhanced glacier retreat on Novaya Zemlya are primarily due to changes on marine-terminating outlet glaciers, as has been widely documented on the Greenland Ice Sheet (e.g. Howat and others, 2008; Moon and others, 2008; Joughin and others, 2010; Murray and others, 2015).

On the Barents Sea Coast, we see remarkable commonality in the timing of changes in marine-terminating outlet glacier frontal positions. Our data show a dramatic reduction in retreat rates for all marine-terminating glaciers between 2013 and 2015 (Figs 2, 6a) and large retreats in 2016: retreat rates on our six focus glaciers in 2016 were between 5 and 10 times greater than the 2015–2020 average (Fig. 6a). Furthermore, detailed seasonal data highlights the common characteristics of the 2016 retreat event between the six glaciers: retreat began and ended on very similar dates and occurred through higher magnitude seasonal ice loss, as opposed to a longer seasonal retreat duration (Fig. 4a). Thus, marine-terminating outlets on the Barents Sea coast show far more commonality in retreat patterns than observed elsewhere across the Arctic (e.g. Murray and others, 2015; Carr and others, 2017b; Bunce and others, 2018; Catania and others, 2018; Fried and others, 2018). Furthermore, they appear to respond rapidly to external forcing, as exemplified by the 2016 retreat event (Figs 4a, d, e). We attribute this to the very short fjords in Novaya Zemlya (Fig. 1), which directly expose the glaciers to variability in adjacent water mass characteristics in the Barents Sea and thus to Atlantic Water (Zeeberg and others, 2001; Polyakov and others, 2017). Finally, our data suggest that periods of frontal position change on the Barents Sea coast are dominated by short-lived (1–2 years) ‘steps’, as opposed to longer-term trends (Fig. 4a). These events are only distinguishable from the high temporal resolution data available from ~2013 onwards but may also have occurred during the period of enhanced retreat between 2000 and 2013.

On the Kara Sea coast, the response of marine-terminating glaciers was less consistent than on the Barents Sea, in both our interannual (Figs 6a, b) and high temporal resolution (Figs 4a, 5a) data. For example, larger seasonal retreats (Figs 5a, b) were observed on certain glaciers in 2018 (KRO and POL) and in 2019 and 2020 on others (ROZE, SHU1, VER, VYL2). Glacier location did not appear to impact the timing of this retreat event: SHU1 and VYL2 retreated in 2019 and 2020 and are located between KRO and POL, which retreated in 2018 (Fig. 1). We suggest that the more limited commonality in patterns of frontal position is due to the glaciers' location on the Kara Sea coast, which is more isolated from Atlantic Water influences and subject to colder water masses from the Arctic (Pavlov and others, 1995; Zeeberg and others, 2001). Thus, our data indicate that glacier-specific factors, such as fjord length and depth, play a greater role on the Kara Sea coast and that ocean warming is not such a dominant driver of retreat as on the Barents Sea coast. Nevertheless, retreat on the Kara Sea coast was still dominated by larger seasonal retreats in specific years (Fig. 4a), indicating some commonality in response to external forcing.

Focusing on seasonal timescales, the timing of the onset of seasonal retreat on Novaya Zemlya is comparable between glaciers and coasts, with retreat most often starting in July, followed by June (Fig. 3, Tables S2, S3). This suggests Novaya Zemlya glaciers respond consistently to forcing at seasonal, as well as annual, timescales. The date on which seasonal retreat ended was more variable on both coasts than the date of onset, but generally fell within a two-month window (Fig. 3, Tables S2, S3), and seasonal retreat usually ended earlier on the Kara Sea coast (most often in October) than on the Barents Sea (most often in November or December). The longer retreat season on the Barents Sea coast most likely reflects the warmer ocean and air temperatures experienced on the Barents Sea coast, due to its exposure to Atlantic water and air masses (e.g. Zeeberg and others, 2001; Politova and others, 2012; Przybylak and others, 2016; Polyakov and others, 2017). This is supported by patterns observed in Greenland, whereby outlet glaciers located in the central- and north-west (i.e. further from the north-Atlantic and higher latitude) tend to retreat from June or July (e.g. Moon and others, 2015; Fried and others, 2018; Carr and others, 2023), compared to May onwards for those located further south, in south-east Greenland (Howat and others, 2010; Schild and others, 2013). Similarly, the variability in the timing of the end of retreat was greater than that for the start of retreat for two north-west Greenland glaciers (Carr and others, 2023).

4.2. What is the relationship between frontal position change and external forcing at seasonal to interannual timescales?

Overall, our data demonstrate a marked shift to warmer air (Figs 6c, f) and ocean (Figs 8c, e, h, j) temperatures and reduced sea ice concentrations (Figs 7c, g) from 2005 onwards. Following this shift to warmer conditions, these climatic factors were characterised by short-term (1–3 year) variations (Figs 6–8) and the magnitude of ocean variability increased markedly after 2010 (Figs 8c, h). These results are consistent with other reports of marked warming in the Barents and Kara Sea region since 2010, which has been termed ‘Atlantification’ (Ma and others, 2017; Tepes and others, 2021a, 2021b; Jakob and Gourmelen, 2023), and has been linked to increased storminess and enhanced Atlantic Water inflow and upwelling (Polyakov and others, 2017; Barton and others, 2018; Lind and others, 2018). We suggest that this shift to warmer climatic and oceanic conditions was responsible for the much higher retreat rates observed on both coasts of Novaya Zemlya from the 2000s onwards (Figs 6–8), which also coincided with elevated glacier velocities on Barents Sea marine-terminating outlets (Melkonian and others, 2016).

On the Barents Sea coast, our data demonstrate a rapid and consistent glacier response to ocean temperatures. Specifically, in 2016, enhanced seasonal retreat (Figs 4a, 8a) coincided with the warmest ocean temperatures in our record, at depths of both 5 and 92 m, and this anomaly was the largest we observed across all climatic/oceanic datasets (Figs 4, 8). The role of ocean temperatures in driving the 2016 retreat is further supported by the warmest anomalies occurring in the summer and autumn (Figs 4d, e), which is consistent with seasonal retreat rates being higher, rather than seasonal retreat persisting until later in the year (Fig. 4a). In addition, the 2016 event is not observed in the land- or lake-terminating glaciers on the Barents Sea coast (Fig. S1), further supporting an oceanic trigger. Comparatively low sea ice concentrations were also observed in 2016, but we suggest they did not play a primary role in driving glacier retreat, since: (i) sea ice concentration anomalies were not as exceptional as those for ocean temperatures (Figs 4c–e) and (ii) studies from Greenland suggest that reduced sea ice concentration may elongate the calving season, by delaying the formation of the seasonal ice mélange and/or causing earlier break up (e.g. Amundson and others, 2010; Robel, 2017; Amundson and others, 2018; Todd and others, 2019). However, this does not correspond with the observation that retreat in 2016 was high in magnitude, but not in duration (Fig. 4a). Although air temperatures were generally warmer from the 2000s onwards, 2016 was not anomalously warm on the Barents Sea coast (Fig. 4b), and high retreat rates in 2016 persisted until November, long after the end of seasonal melt on Novaya Zemlya in mid-August (Zhao and others, 2014). This suggests that air temperatures are not directly driving marine-terminating outlet glacier retreat. However, they may play a role via the enhancement of ocean melting due to plume activity (e.g. Vieli and others, 2011; Straneo and others, 2012; Straneo and others, 2013; Catania and others, 2020) and/or via their longer-term impacts on ice thinning rates and surface mass balance (discussed in Section 4.3). Enhanced glacier retreat also coincided with high rates of mass loss from the RHA as a whole (Jakob and Gourmelen, 2023), suggesting that the observed changes in oceanic forcing had widespread impacts across the region and were the primary trigger of observed retreat.

High resolution frontal position data were not available for the pause in retreat between 2013 and 2015. However, ocean temperatures were anomalously cool at both 5 and 92 m depth in 2014 (Figs 4d, e) and this was the only ‘cool’ anomaly from 2007 onwards (Figs 8c, e). Thus, we suggest that ocean temperatures were likely responsible for the 2013–2015 pause, as well as the high retreat rates in 2016. This is likely due to a combination of enhanced basal melting across the glacier fronts, due to warmer ocean temperatures, and undercutting of the glacier termini leading to enhanced calving, (e.g. Benn and others, 2007; Straneo and others, 2012; O'Leary and others, 2013; Straneo and others, 2013; Todd and others, 2019) given that the glaciers do not have substantial floating sections (Carr and others, 2014). These cooler ocean temperatures and reduced retreat rates also coincided with low rates of ice loss across the Russian High Arctic as a whole in 2014 and 2015, (Wouters and others, 2019; Hugonnet and others, 2021; Jakob and Gourmelen, 2023). Overall, our data suggest that ocean temperatures are the primary control on marine-terminating outlet glacier retreat in the Barents Sea and that the glacier termini respond very rapidly (with ~1 year) to changes in forcing. Furthermore, previous work has documented a 1.5-year response time between changes in ocean and/or atmospheric forcing and thinning rates on Novaya Zemlya (Tepes and others, 2021b). This suggests that the Novaya Zemlya is highly responsive to climate forcing, via both terminus retreat and glacier thinning, particularly on the Barents Sea coast.

On the Kara Sea coast, the relationship between forcing factors and glacier retreat is less clear than for the Barents Sea. As on the Barents Sea, our data show a shift to warmer air and ocean temperatures and reduced sea ice concentrations from the mid-2000s onwards (Figs 6–8). However, the main phases of retreat in the high-resolution terminus data occur in 2018, 2019 and 2020, when ocean and air temperatures on the Kara Sea coast were comparable to those in 2013 and 2014, when glacier retreat was substantially reduced (Fig. 5). Furthermore, warm ocean temperatures in 2016 do not correspond to enhanced glacier retreat on the Kara Sea coast (Fig. 5). We suggest that this difference in glacier responsiveness is due to the Kara Sea glaciers being more remote from warm Atlantic Water than those on the Barents Sea coast. This is consistent with previous work which has suggested that ocean forcing explains 59% of the difference in mass loss between the Barents and Kara Sea coasts, between 2010 and 2018 (Tepes and others, 2021b). Thus, we suggest that ocean warming is likely to substantially influence retreat on the Kara Sea coast, but its role is less clear and less immediate than on the Barents Sea.

The importance of ocean forcing in driving Novaya Zemlya glacier retreat agrees with previous findings from three major Svalbard glaciers, which suggested that ocean temperatures were the dominant control on frontal ablation and found limited correlation between air temperature and frontal ablation (Luckman and others, 2015; Holmes and others, 2019). Furthermore, there was spatial variability in the relationship between ocean temperature and frontal ablation, which was attributed to differences in Atlantic Water access to the glacier termini (Luckman and others, 2015), and we suggest a similar explanation for the observed differences in responsiveness to ocean forcing between the Barents and the Kara Sea coasts observed in our study (Figs 4, 5, 8). Results from LeConte Glacier, Alaska, have also highlighted the substantial contribution of submarine melt rates to total ice loss and demonstrated that seasonal changes in terminus position tracked changes in water temperatures (Motyka and others, 2003), as we observed on the Barents Sea coast glaciers (Fig. 4). In Greenland, previous work has underscored the role of warm Atlantic Water in driving glacier retreat and contributing to spatial variability in glacier behaviour (e.g. Rignot and others, 2010; Morlighem and others, 2017; Catania and others, 2018; Catania and others, 2020). However, the response to ocean forcing we observed on the Barents Sea coast is both more consistent between individual glaciers and faster than observed in Greenland (e.g. Murray and others, 2015; Carr and others, 2017b; Bunce and others, 2018; Catania and others, 2018; Fried and others, 2018). Furthermore, the apparent sensitivity of Novaya Zemlya outlet glaciers to oceanic forcing contrasts with findings in the Canadian Arctic, where regional retreat showed no correlation with ocean forcing and, instead, correlated with melt rates (Cook and others, 2019). Thus, we suggest that Novaya Zemlya glaciers, particularly on the Barents Sea coast, respond to ocean forcing, as observed in Svalbard, Alaska and Greenland, but that this response is quicker and more consistent than observed in Greenland. This has important implications of Novaya Zemlya's near-future contribution to sea level rise, as its glaciers may rapidly retreat if ocean temperatures continue to warm, whereas these glaciers could, in principle, reduce their contribution to sea level rise in the unlikely event that ocean temperatures cool.

4.3. Does glacier retreat correspond to thinning and what is the direction of this relationship?

Overall, our results suggest that retreat was primarily triggered by warmer ocean temperatures and that episodes of ocean warming determined the timing of major retreat phases, particularly on the Barents Sea coast (Figs 4, 5). However, analysis of thinning data suggests that multi-year thinning may have preconditioned glaciers to retreat in response to this ocean forcing and may also have contributed to the variation in response between individual glaciers. Specifically, our results show a strong relationship between glacier retreat and thinning (Fig. 9). Furthermore, the data suggest that thinning during the preceding ~5 years promotes retreat, whereas retreat only impacts thinning within ~1 year (Fig. 10). We therefore suggest that multi-year thinning, potentially in response to warmer air temperatures since 2005 (Fig. 6) (Tepes and others, 2021b) has preconditioned glaciers on Novaya Zemlya to retreat in response to oceanic warming, by bringing the glaciers closer to flotation and making them more vulnerable to full-thickness fracture and calving (e.g. van der Veen, 1996, 1998, 2002, 2007; Pfeffer, 2007; Joughin and others, 2011; Vieli and others, 2011). Our data suggest that glacier retreat correlates with thinning within ~1 year (Fig. 10), most likely as a result of a reduction in buttressing due to loss of a section of the grounded terminus, which in turn leads to increased ice flow and dynamic thinning (e.g. Joughin and others, 2004). However, thinning data are only available between 2011 and 2020, meaning that we may be missing the longer-term impacts of glacier retreat on dynamic thinning, which can operate over multidecadal timescales on the Greenland Ice Sheet (e.g. Pritchard and others, 2009; Price and others, 2011). Furthermore, our data document the relationship between terminus position, surface elevation change and climatic and oceanic forcing. However, they do not provide insight into the mechanisms by which this ice loss occurred, i.e. the extent to which changes can be attributed to surface mass balance vs calving. Unfortunately, accurate data on basal topography and offshore bathymetry for Novaya Zemlya are very scarce, which limits our ability to investigate calving rates from Novaya Zemlya's outlet glaciers. A potential avenue for future research is to generate basal topography by inverting from high resolution ice velocity data, although the lack of directly measured data to verify these ice thicknesses may limit their reliability and hence our abilities to accurately determine calving rates.

Although we observe commonality in the timing of marine-terminating outlet glacier retreat, particularly on the Barents Sea coast, the magnitude of this retreat varies greatly and this variability is highest during periods of rapid retreat, i.e. 2000–2012 and 2015–2020 (Figs 2, 6a, b). This indicates that glacier specific factors are influencing the magnitude of glacier retreat, after retreat is initially triggered by ocean warming. Previous work on Novaya Zemlya has identified the role of fjord width in modulating retreat on Novaya Zemlya (Carr and others, 2014, 2017a), and the role of bed topography has been highlighted in numerous studies on the Greenland Ice Sheet (e.g. Porter and others, 2014; Carr and others, 2015; Catania and others, 2018; Hill and others, 2018), although Novaya Zemlya's basal topography is largely unquantified, making its role difficult to quantify. However, a further potential contributor to variability in retreat rates between individual glaciers is variability in thinning rates: despite an overall significant correlation between glacier thinning and retreat, we observed spread in this relationship (Fig. 9). We suggest that this spread is due to the specific fjord bathymetry and/or subglacial topography of each glacier. For example, glaciers that have much higher rates of retreat compared to thinning (E.g. KRA, KRI, KRO and VIS; Fig. 9) may have already thinned to floatation and/or off pinning points, meaning that further thinning has limited impact on retreat. Conversely, the impact of thinning on glacier retreat may be lower on specific glaciers (e.g. INO, KRA2, RYK, VYL1) if thinning is partially offset by e.g. glacier retreat into narrower and/or shallower fjords. As such, while Novaya Zemlya glaciers showed commonality in the timing of phases of frontal position change, which we suggest was triggered by ocean forcing, the magnitude of their response is likely modulated by glacier specific factors and multi-year thinning patterns.

5. Conclusions

Our data demonstrate that the pause in retreat on Novaya Zemlya's marine-terminating outlet glaciers between 2013 and 2015 has not persisted and that rapid retreat rates observed in the 2000s have resumed. We see remarkable consistency in the timing of frontal positional change on the Barents Sea coast, both at seasonal and interannual timescales, which we attribute to the glaciers' short fjords and exposure to variability in Atlantic Water conditions immediately offshore. In 2016, glaciers on the Barents Sea coast experienced a large and synchronous retreat event, which occurred via enhanced rates of ice loss, rather than an extended period of seasonal retreat. We attribute this event to anonymously warm ocean temperatures during summer and autumn 2016. We suggest that air temperatures did not directly drive this retreat event, as they were not anomalously warm in 2016 and high retreat rates persisted long past the end of the melt season on Novaya Zemlya. We also tentatively attribute the 2013–2015 pause in retreat to cool ocean temperatures in 2014. On the Kara Sea, interannual changes in glacier frontal position were less synchronous than on the Barents Sea and the response to forcing was less clear, most likely due to their comparative isolation from warm waters in the North Atlantic. Our data suggests that multi-year thinning may condition Novaya Zemlya's outlet glaciers to retreat, in response to periods of warmer ocean temperatures, meaning that warmer air temperatures since the 2000s may be indirectly driving glacier retreat via thinning. Glacier frontal position change predominantly impacts thinning within 1 year, most likely due to the immediate impact of loss of ice buttressing. However, thinning data were only available for 2011–2020, and so longer-term thinning data are needed to investigate whether retreat can cause multi-decadal thinning on Novaya Zemlya, as observed in Greenland. Although glacier frontal position changes showed commonality in timing, the magnitude of the subsequent retreats varied greatly, likely due to glacier specific factors, such as basal topography, fjord bathymetry and thinning rates. Overall, our results show that marine-terminating outlet glaciers on the Barents Sea coast of Novaya Zemlya respond rapidly and synchronously to ocean forcing and may be pre-conditioned for retreat by multi-annual surface thinning, likely due to multi-annual air temperature increases observed in the region since 2005.