Exceptional retreat of Novaya Zemlya’s marine-terminating outlet glaciers between 2000 and 2013

. Novaya Zemlya (NVZ) has experienced rapid ice loss and accelerated marine-terminating glacier retreat during the past 2 decades. However, it is unknown whether this retreat is exceptional longer term and/or whether it has persisted since 2010. Investigating this is vital, as dynamic thinning may contribute substantially to ice loss from NVZ, but is not currently included in sea level rise predictions. Here, we use remotely sensed data to assess controls on NVZ glacier retreat between 1973/76 and 2015. Glaciers that terminate into lakes or the ocean receded 3.5 times faster than those that terminate on land. Between 2000 and 2013, retreat rates were signiﬁcantly higher on marine-terminating outlet glaciers than during the previous 27 years, and we observe widespread slowdown in retreat, and even advance, between 2013 and 2015. There were some common patterns in the timing of glacier retreat, but the magnitude varied between individual glaciers. Rapid retreat between 2000 and 2013 corresponds to a period of signiﬁcantly warmer air temperatures and reduced sea


Introduction
Glaciers and ice caps are the main cryospheric source of global sea level rise and contributed approximately −215 ± 26 Gt yr −1 between 2003 and 2009 . This ice loss is predicted to continue during the 21st century (Meier et al., 2007;Radić et al., 2014), and changes are expected to be particularly marked in the Arctic, where warming of up to 8 • C is forecast (IPCC, 2013). Outside of the Greenland Ice Sheet, the Russian high Arctic (RHA) accounts for approximately 20 % of Arctic glacier ice (Dowdeswell and Williams, 1997;Radić et al., 2014) and is, therefore, a major ice reservoir. It comprises three main archipelagos: Novaya Zemlya (NVZ; glacier area = 21 200 km 2 ), Severnaya Zemlya (16 700 km 2 ), and Franz Josef Land (12 700 km 2 ) (Moholdt et al., 2012). Between 2003 and, these glaciated regions lost ice at a rate of between 9.1 Gt a −1 (Moholdt et al., 2012) and 11 Gt a −1 , with over 80 % of mass loss coming from Novaya Zemlya (NVZ) (Moholdt et al., 2012). This much larger contribution from NVZ has been attributed to it experiencing longer melt seasons and high snowmelt variability between 1995 and 2011 (Zhao et al., 2014). More recent estimates suggest that the mass balance of the RHA was −6.9 ± 7.4 Gt between 2004(Matsuo and Heki, 2013 and that thinning rates increased to −0.40 ± 0.09 m a −1 between 2012/13 and 2014, compared to the long-term average of −0.23 ± 0.04 m a −1 (1952 and 2014) (Melkonian et al., 2016). The RHA is, therefore, following the Arctic-wide pattern of negative mass balance  and glacier retreat that has been observed in Greenland (Enderlin et al., 2014;McMillan et al., 2016), Svalbard (Moholdt et al., 2010a, b;Nuth et al., 2010), and the Canadian Arctic (Enderlin et al., 2014;McMillan et al., 2016). However, the RHA has been studied far less than Published by Copernicus Publications on behalf of the European Geosciences Union.
other Arctic regions, despite its large ice volumes. Furthermore, assessment of 21st-century glacier volume loss highlights the RHA as one of the largest sources of future ice loss and contribution to sea level rise, with an estimated loss of 20-28 mm of sea level rise equivalent by 2100 (Radić et al., 2014).
Arctic ice loss occurs via two main mechanisms: a net increase in surface melting, relative to surface accumulation, and accelerated discharge from marine-terminating outlet glaciers (e.g. Enderlin et al., 2014;. These marine-terminating outlets allow ice caps to respond rapidly to climatic change, both immediately through calving and frontal retreat (e.g. Blaszczyk et al., 2009;Carr et al., 2014;McNabb and Hock, 2014;Moon and Joughin, 2008) and also through long-term drawdown of inland ice, often referred to as "dynamic thinning" (e.g. Price et al., 2011;Pritchard et al., 2009). During the 2000s, widespread marine-terminating glacier retreat was observed across the Arctic (e.g. Blaszczyk et al., 2009;Howat et al., 2008;Mc-Nabb and Hock, 2014;Moon and Joughin, 2008;Nuth et al., 2007), and substantial retreat occurred on Novaya Zemlya between 2000 and 2010 (Carr et al., 2014): retreat rates increased markedly from around 2000 on the Barents Sea coast and from 2003 on the Kara Sea (Carr et al., 2014). Between 1992 and 2010, retreat rates on NVZ were an order of magnitude higher on marine-terminating glaciers (−52.1 m a −1 ) than on those terminating on land (−4.8 m a −1 ) (Carr et al., 2014), which mirrors patterns observed on other Arctic ice masses (e.g. Dowdeswell et al., 2008;Moon and Joughin, 2008;Pritchard et al., 2009;Sole et al., 2008) and was linked to changes in sea ice concentrations (Carr et al., 2014). However, the pattern of frontal-position changes on NVZ prior to 1992 is uncertain, and previous results indicate different trends, dependant on the study period: some studies suggest glaciers were comparatively stable or retreating slowly between 1964 and 1993 (Zeeberg and Forman, 2001), whilst others indicate large reductions in both the volume (Kotlyakov et al., 2010) and the length of the ice coast (Sharov, 2005) from ∼ 1950 to 2000, and longer-term retreat (Chizov et al., 1968;Koryakin, 2013;Shumsky, 1949). Consequently, it is difficult to contextualize the observed period of rapid retreat from ∼ 2000 until 2010 (Carr et al., 2014) and to determine if it was exceptional or part of an ongoing trend. Furthermore, it is unclear whether glacier retreat has continued to accelerate after 2010, and hence further increased its contribution to sea level rise, or whether it has persisted at a similar rate. This paper aims to address these limitations, by extending the time series of glacier frontal-position data on NVZ to include the period 1973/76 to 2015, which represents the limits of available satellite data.
Initially, surface elevation change data from NVZ suggested that there was no significant difference in thinning rates between marine-and land-terminating outlet glacier catchments between 2003 and 2009 (Moholdt et al., 2012). This contrasted markedly with results from Greenland (e.g. Price et al., 2011;Sole et al., 2008) but was similar to the Canadian Arctic, where the vast majority of recent ice loss occurred via increased surface melting (∼ 92 % of total ice loss), rather than accelerated glacier discharge (∼ 8 %) (Gardner et al., 2011). This implied that outlet glacier retreat was having a limited and/or delayed impact on inland ice or that available data were not adequately capturing surface elevation change in outlet glacier basins (Carr et al., 2014). More recent results demonstrate that thinning rates on marine-terminating glaciers on the Barents Sea coast are much higher than on their land-terminating neighbours, suggesting that glacier retreat and calving do promote inland, dynamic thinning (Melkonian et al., 2016). However, higher melt rates also contributed to surface lowering, evidenced by the concurrent increase in thinning observed on landterminating outlets (Melkonian et al., 2016). High rates of dynamic thinning have also been identified on Severnaya Zemlya, following the collapse of the Matusevich Ice Shelf in 2012 (Willis et al., 2015). Here, thinning rates increased to 3-4 times above the long-term average , following the ice-shelf collapse in summer 2012, and outlet glaciers feeding into the ice shelf accelerated by up to 200 % (Willis et al., 2015). The most recent evidence, therefore, suggests that NVZ and other Russian high Arctic ice masses are vulnerable to dynamic thinning, following glacier retreat and/or ice-shelf collapse. Consequently, it is important to understand the longer-term retreat history on NVZ in order to evaluate its impact on future dynamic thinning. Furthermore, we need to assess whether the high glacier retreat rates observed on NVZ during the 2000s have continued and/or increased, as this may lead to much larger losses in the future and may indicate that a step change in glacier behaviour occurred in ∼ 2000.
In this paper, we use remotely sensed data to assess glacier frontal-position change for all major (> 1 km wide) Novaya Zemlya outlet glaciers (Fig. 1). This includes all outlets from the ice cap of the northern island (hereafter referred to as the northern island ice cap for brevity) and its subsidiary ice fields (Fig. 1). We were unable to find the names of these subsidiary ice fields in the literature, so we name them Sub 1 and Sub 2 (Fig. 1). A total of 54 outlet glaciers were investigated, which allowed us to assess the impact of different glaciological, climatic and oceanic settings on retreat (Table S1 in the Supplement). Specifically, we assessed the impact of coast (Barents Sea versus Kara Sea on the northern ice mass), ice mass (northern island ice cap, Sub 1, or Sub 2), terminus type (marine-, lake-, and land-terminating), and latitude (Table 1). The two coasts of Novaya Zemlya are characterized by very different climatic and oceanic conditions: the Barents Sea coast is influenced by water from the North Atlantic (Loeng, 1991;Pfirman et al., 1994;Politova et al., 2012) and subject to Atlantic cyclonic systems (Zeeberg and Forman, 2001), which results in warmer air and ocean temperatures as well as higher precipitation (Przybylak and Wyszyński, 2016;Zeeberg and Forman, 2001 contrast, the Kara Sea coast is isolated from North Atlantic weather systems, by the topographic barrier of NVZ (Pavlov and Pfirman, 1995), and is subject to cold, Arctic-derived water, along with much higher sea ice concentrations (Zeeberg and Forman, 2001). We therefore aim to investigate whether these differing climatic and oceanic conditions lead to major differences in glacier retreat between the two coasts. Glaciers identified as surge type (Grant et al., 2009) were excluded from the retreat calculations and analysis. However, frontalposition data are presented separately for three glaciers that were actively surging during the study period. Glacier retreat was assessed from 1973/76 to 2015 in order to provide the greatest temporal coverage possible from satellite imagery. We use these data to address the following questions: 1. At multi-decadal timescales, is there a significant difference in glacier retreat rates according to (i) terminus type (land-, lake-or marine-terminating); (ii) coast  (Fig. 1). The northern island ice cap contains the vast majority of ice (19 841 km 2 ) and the majority of the main outlet glaciers (Fig. 1). The northern island also has two smaller ice fields, Sub 1 and Sub 2, which are much smaller in area (1010 and 705 km 2 respectively) and have far fewer, smaller outlet glaciers (Sub 1 = 4; Sub 2 = 5) ( Fig. 1). All glaciers that have been previously identified as surge type and those smaller than 1 km in width were excluded from our main analysis of glacier retreat rates and response to climate forcing. However, we also observed three glaciers surging during the study period: Anuchina (ANU), Mashigina (MAS), and Serp i Molot (SER) (Fig. 1). MAS and SER have been previously identified as surge type (Grant et al., 2009), but our data provide better constraints on the duration and timing of these surges. ANU was identified as potentially surge type, on the basis of looped moraines (Grant et al., 2009). Our study confirms it as surge type and provides information on the surge timing and duration. These three glaciers are not included in the assessment of NVZ glacier response to climate change, as surging can occur impudently of climate forcing (Meier and Post, 1969), but are discussed separately to improve our knowledge of NVZ surge characteristics. This resulted in a total of 54 outlet glaciers, which were located in a variety of settings and hence allowed us to assess spatial controls on glacier retreat (Table 1). Where available glacier names and World Glacier Inventory IDs are given in Table S1 in the Supplement, along with glacier acronyms used in this paper. The impact of coast could only be assessed for the main ice mass, as the glaciers on the smaller ice masses, Sub 1 and Sub 2, are located on the southern ice margin and so do not fall on either coast (Fig. 1).

Glacier frontal position
Outlet glacier frontal positions were acquired predominantly from Landsat imagery. These data have a spatial resolution of 30 m and were obtained freely via the United States Geological Survey (USGS) Global Visualization Viewer (Glo-Vis) (http://glovis.usgs.gov/). The frequency of available im-agery varied considerably during the study period. Data were available annually from 1999 to 2015 and between 1985 and 1989, although georeferencing issues during the latter time period meant that imagery needed to be re-coregistered manually using stable, off-ice locations as tie points. Prior to 1985, the only available Landsat scenes dated from 1973, and these also needed to be manually georeferenced. We verified all images that required georeferencing against Landsat 8 data, which should have the most accurate location data of the imagery time series. We did this by comparing the location of features that should be static between images (e.g. large rock fractures) and also checking for any unrealistic changes in the lateral glacier margins, over and above what could be expected by glacier melting. Any images where we saw changes in the location of static features above the image resolution were not used. As such, orthorectification was not required for these images, as the terrain is relatively gentle on NVZ, and our verification process showed that the images were co-located with the Landsat 8 imagery to within a pixel using just georeferencing. Hexagon KH-9 imagery was used to determine frontal positions in 1976 and 1977, but full coverage of the study area was not available for either year. The data resolution is 20 to 30 ft (∼ 6-9 m). The earliest common date for which we have frontal positions for all glaciers is 1986, and so we calculate total retreat rates for the period 1986-2015 and use these values to assess spatial variability in glacier recession across the study region. All glacier frontal positions are calculated relative to 1986 (i.e. the frontal position in 1986 = 0 m) to allow for direct comparison.
Due to the lack of Landsat imagery during the 1990s, we use synthetic aperture radar (SAR) image mode precision data during this period. The data were provided by the European Space Agency, and we use European Remote-sensing Satellite-1 (ERS-1) and ERS-2 products (https://earth.esa. int/web/guest/data-access/browse-data-products/-/asset_ publisher/y8Qb/content/sar-precision-image-product-1477). Following Carr et al. (2013b), the ERS scenes were first co-registered with Envisat imagery and then processed using the following steps: apply precise orbital state vectors; radiometric calibration; multi-look; and terrain correction. This gave an output resolution of 37.5 m, which is comparable to Landsat. For each year and data type, imagery was acquired as close as possible to 31 July to minimize the impact of seasonal variability. However, this is unlikely to substantially effect results, as previous studies suggest that seasonal variability in terminus position is very limited on NVZ (∼ 100 m a −1 ) (Carr et al., 2014) and is therefore much less than the interannual and inter-decadal variability we observe here. Glacier frontal-position change was calculated using the box method: the terminus was repeatedly digitized from successive images, within a fixed reference box, and the resultant change in area is divided by the reference box width to get frontal-position change (e.g. Moon and Joughin, 2008). Following previous studies (Carr et al., 2014), we The Cryosphere, 11, 2149-2174, 2017 www.the-cryosphere.net/11/2149/2017/ determined the frontal-position errors for marine-and laketerminating outlets glaciers by digitizing 10 sections of rock coastline from six images, evenly spread through the time series (1976, 1986, 2000, 2005, 2010, and 2015) and across NVZ. The resultant error was 17.5 m, which equates to a retreat rate error of 1.75 m a −1 at the decadal time intervals discussed here. The terminus is much harder to identify on land-terminating outlet glaciers due to the similarity between the debris-covered ice margins and the surrounding land, which adds an additional source of error. We quantified this by re-digitizing a sub-sample of six land-terminating glaciers in each of the six images, which were spread across NVZ. The additional error for land-terminating glaciers was 66.1 m, giving a total error of 68.4 m, which equates to a retreat rate error of 6.86 m a −1 for decadal intervals.

Climate and ocean data
Air temperature data were obtained from meteorological stations located on, and proximal to, Novaya Zemlya ( Fig. 1). Directly measured meteorological data are very sparse on NVZ, and there are large gaps in the time series for many stations. We use data from two stations, Malye Karmakuly (WMO ID: 20744) and E. K. Fedorova (WMO ID: 20946), as these are the closest stations to the study glaciers that have a comprehensive (although still not complete) record during the study period (Supplement Table S2). The data were obtained from the Hydrometeorological Information -World Data Centre Baseline Climatological Data Sets (http: //meteo.ru/english/climate/cl_data.php) and were provided at a monthly temporal resolution. For each station, we calculated meteorological seasonal means (December-February, March-May, June-August, September-November) in order to assess the timing of any changes in air temperature, as warming in certain seasons would have a different impact on glacier retreat rates. Seasonal and annual means were only calculated if values were available for all months. Due to data gaps, particularly from 2013 onwards (Supplement Table S2), we also assess changes in air temperature using ERA-Interim reanalysis data (http://www.ecmwf.int/ en/research/climate-reanalysis/era-interim). We use temperature data from the surface (2 m elevation) and 850 hPa pressure level, as these are likely to be a good proxy for meltwater availability (X. Fettweis, personal communication, 2017). We use the "monthly means of daily means" product for all months between 1979 and 2015. As with the meteorological stations, we calculate means for the meteorological seasons and annual means. Sea ice data were acquired from the Nimbus-7 SMMR and DMSP SSM/I-SSMIS Passive Microwave data set (https: //nsidc.org/data/docs/daac/nsidc0051_gsfc_seaice.gd.html). The data provide information on the percentage of the ocean covered by sea ice, and this is measured using brightness temperatures from microwave sensors. The data have a spatial resolution of 25 × 25 km, and we use the monthly-averaged product. This data set was selected due to its long temporal coverage, which extends from 26 October 1978 to 31 December 2015 and thus provides a consistent data set throughout our study period. NVZ glaciers are not located within long fjords and are relatively exposed to the open ocean ( Fig. 1). Consequently, sea ice conditions within 25 km of the glacier fronts (i.e. the data resolution) are likely to be reasonably representative of the overall sea ice trends experienced by the glaciers, particularly at the decadal timescales assessed here. However, it should be noted that the data cannot provide detailed information on sea ice conditions specific to each glacier front, but they are used here as they comprise the only data set available for the entire study period. Monthly sea ice concentrations were sampled from the grid squares closest to the study glaciers and were split according to coast (i.e. Barents Sea and Kara Sea). From the monthly data, we calculated seasonal means and the number of ice-free months, which we define as the number of months where the mean monthly sea ice cover is less than 10 %.
Data on the North Atlantic Oscillation (NAO) were obtained from the Climatic Research Unit (https://crudata. uea.ac.uk/cru/data/nao/), and the monthly product was used. This records the normalized pressure difference between Iceland and the Azores (Hurrell, 1995). Arctic Oscillation (AO) data were acquired from the Climate Prediction Center (http://www.cpc.noaa.gov/products/ precip/CWlink/daily_ao_index/teleconnections.shtml). The AO is characterized by winds at 55 • N, which circulate anticlockwise around the Arctic (e.g. Higgins et al., 2000;Zhou et al., 2001). The AO index is calculated by projecting the AO loading pattern onto the daily anomaly 1000 mbar height field, at 20-90 • N latitude (Zhou et al., 2001). The Atlantic Multidecadal Oscillation (AMO) is a mode of variability associated with averaged, de-trended sea surface temperatures (SSTs) in the North Atlantic and varies over timescales of 60 to 80 years (Drinkwater et al., 2013;Sutton and Hodson, 2005). Monthly data were downloaded from the National Oceanic and Atmospheric Administration (https://www.esrl. noaa.gov/psd/data/timeseries/AMO/).
We use ocean temperature data from the "Climatological Atlas of the Nordic Seas and Northern North Atlantic" (Hurrell, 1995;Korablev et al., 2014) (https://www.nodc.noaa. gov/OC5/nordic-seas/). The atlas compiles data from over 500 000 oceanographic stations, located across the Nordic Seas, between 1900 and 2012. It provides gridded climatologies of water temperature, salinity, and density, at a range of depths (surface to 3500 m), for the region bounded by 83.875 to 71.875 • N and 47.125 • W to 57.875 • E. Here, we use data from the surface and 100 m depth to capture changes in ocean temperatures at different depths: surface warming may influence glacier behaviour through changes in sea ice and/or undercutting at the waterline (Benn et al., 2007), whereas warming in the deeper layers can enhance sub-aqueous melting (Sutherland et al., 2013) 100 m was chosen, as it is the deepest level that includes the majority of the continental shelf immediately offshore of Novaya Zemlya. Further details of the data set production and error values are given in Korablev et al. (2014). We use the decadal ocean temperature product to identify broadscale changes, which is provided at the following time intervals: 1971-1980, 1981-1990, 1991-2000, and 2001-2012. We use the decadal product as there are few observations offshore of Novaya Zemlya during the 2000s, whereas the data coverage is much denser in the 1980s and 1990s (a full inventory of the number and location of observations for each month and year is provided here: https://www.nodc. noaa.gov/OC5/nordic-seas/atlas/inventory.html). As a result, maps of temperature changes in the 2000s are produced using comparatively data few points, meaning that they may not be representative of conditions in the region and that directly comparing data at a shorter temporal resolution (e.g. annual data) may be inaccurate. Furthermore, the input data were measured offshore of Novaya Zemlya and not within the glacier fjords. Consequently, there is uncertainty over the extent to which offshore warming is transmitted to the glacier front and/or the degree of modification due to complexities in the circulation and water properties within glacial fjords.
We therefore use decadal-scale data to gain an overview of oceanic changes in the region, but we do not attempt to use them for detailed analysis of the impact of ocean warming at the glacier front, nor for statistical testing.

Statistical analysis
We used a Kruskal-Wallis test to investigate statistical differences in total retreat rate  for the different categories of outlet glacier within our study population, i.e. terminus type (marine-, land-, and lake-terminating), coast (Barents Sea and Kara Sea), and ice mass (northern island ice cap, Sub 1, and Sub 2). The Kruskal-Wallis test is a nonparametric version of the one-way ANOVA (analysis of variance) test and analyses the variance using the ranks of the data values, as opposed to the actual data. Consequently, it does not assume normality in the data, which is required here, as Kolmogorov-Smirnov tests indicate that total retreat rate  is not normally distributed for any of the glacier categories (e.g. terminus type). This is also the case when we test for normality at each of the four time intervals discussed below (1973/76-1986, 1986-2000, 2000-2013, and 2013-2015). The Kruskal-Wallis test gives a p value for the null hypothesis that two or more data samples come from the same population. As such, a large p value suggests it is likely that the samples come from the same population, whereas a small value indicates that this is unlikely. We follow convention and use a significance value of 0.05, meaning that a p value of less than or equal to 0.05 indicates that the data samples are significantly different.
We assessed the influence of glacier latitude on total retreat rate , using simple linear regression. This fits a line to the data points and gives an R 2 value and a p value for this relationship. The R 2 value indicates how well the line describes the data: if all points fell exactly on the line, the R 2 would equal 1, whereas if the points were randomly distributed about the line, the R 2 would equal 0. The pvalue tests the null hypothesis that the regression coefficient is equal to zero, i.e. that the predictor variable (e.g. glacier catchment size) has no relationship to the response variable (e.g. total glacier retreat rate). A p value of 0.05 or less therefore indicates that the null hypothesis can be rejected and that the predictor variable is related to the response variable (e.g. glacier latitude is related to glacier retreat rate). The residuals for these regressions were normally distributed. However, we also regressed catchment area against total retreat rate, and the regression residuals were not normally distributed, indicating that it is not appropriate to use regression in this case. Consequently, we used Spearman's rank correlation coefficient, which is non-parametric and therefore does not require the data to be normally distributed. Catchments were obtained from Moholdt et al. (2012).
Wilcoxon tests were used to assess significant differences in mean glacier retreat rates between four time intervals : 1973/76-1986, 1986-2000, 2000-2013, and 2013-2015. These intervals were chosen through manual assessment of apparent breaks in the data. For each interval, data were split according to terminus type (marine, land, and lake), and marine-terminating glaciers were further subdivided by coast (Barents Sea and Kara Sea). For each category, we then used the Wilcoxon test to determine whether mean retreat rates for all of the glaciers during one time period (e.g. 1986-2000) were significantly different from those for another time period (e.g. 2000-2013). The Wilcoxon test was selected as it is non-parametric and our retreat data are not normally distributed, and it is suitable for testing statistical difference between data from two time periods (Miles et al., 2013). As with the Kruskal-Wallis test, a p value of less than or equal to 0.05 is taken as significant and indicates that the two time periods are significantly different. We also used the Wilcoxon test to identify any significant differences in mean air temperatures and sea ice conditions for the same time intervals as glacier retreat to allow for direct comparison. For the first time interval (1973/76-1986), we use air temperature data from 1976 to 1986 from the meteorological stations, but the sea ice and ERA-Interim data are only available from 1979. The statistical analysis was done separately for sea ice on the Barents Sea and Kara Sea coast and using meteorological data from Malye Karmakuly and E. K. Fedorova (Fig. 1). ERA-Interim data were analysed as a whole, as the spatial resolution of the data does not allow us to distinguish between the two coasts. In each case, we compared seasonal means for each year of a certain time period with the seasonal means for the other time period (e.g. 1976-1985 versus 2000-2012). For the sea ice data, we used calendar seasons (January-March, April-June, July-September, October-December), which fits with the Arctic sea ice min- The Cryosphere, 11, 2149-2174, 2017 www.the-cryosphere.net/11/2149/2017/ ima in September and maxima in March. For the air temperature data, meteorological seasons (December-February, March-May, June-August, September-November) are more appropriate. We also tested mean annual air temperatures and the number of sea-ice-free months.
In order to further investigate the temporal pattern of retreat on Novaya Zemlya, we use statistical change-point analysis (Eckley et al., 2011). We applied this to our frontalposition data for marine-and lake-terminating glaciers, and to the sea ice and air temperature data. Land-terminating glaciers are not included, due to the much higher error margins compared to any trends, which could lead to erroneous change-points being identified. Change-point analysis allows us to automatically identify significant changes in the time series data and whether there has been a shift from one mode of behaviour to another (e.g. from slower to more rapid retreat) (Eckley et al., 2011). Formally, a change-point is a point in time where the statistical properties of prior data are different from the statistical properties of subsequent data; the data between two change-points are a segment. There are various ways that one can determine when a change-point should occur, but the most appropriate approach for our data is to consider changes in regression.
In order to automate the process, we use the cpt.reg function in the R EnvCpt package (Killick et al., 2016) with a minimum number of four data points between changes. This function uses the pruned exact linear time (PELT) algorithm (Killick et al., 2012) from the change-point package (Killick and Eckley, 2015) for fast and exact detection of multiple changes. The function returns change-point locations and estimates of the intercept and slope of the regression lines between changes. We give the algorithm no information on when we might be expecting a change or how large it may be, allowing it to automatically determine statistically different parts of the data. In this way, we use the analysis to determine whether, and when, retreat rates change significantly on each of the marine-and lake-terminating glaciers on NVZ, and whether there are any significant breaks in our sea ice and air temperature data. We also apply the change-point analysis to the number of sea-ice-free months, but as the data do not contain a trend, we identify breaks using significant changes in the mean, rather than a change in regression. Thus, we can identify any common behaviour between glaciers, determine the timing of any common changes, and compare this to any significant changes in atmospheric temperatures and sea ice concentrations.

Spatial controls on glacier retreat
The Kruskal-Wallis test was used to identify significant differences in total retreat rate  for glaciers located in different settings. First, terminus type was in-vestigated. Results demonstrated that total retreat rates  were significantly higher on lake-and marineterminating glaciers than those terminating on land, at a very high confidence interval (< 0.001) (Fig. 2). Retreat rates were 3.5 times higher on glaciers terminating in water (lake = −49.1 m a −1 ; marine = −46.9 m a −1 ) than those ending on land (−13.8 m a −1 ) (Fig. 2). In contrast, there was no significant difference between lake-and marineterminating glaciers (Fig. 2). Next, we assessed the role of coastal setting (i.e. Barents Sea versus Kara Sea) as climatic and oceanic conditions differ markedly between the two coasts. When comparing glaciers with the same terminus type, there was no significant difference in retreat rates between the two coasts ( Fig. 2: p value = 0.178 for marineterminating glaciers, and p value = 1 for land-terminating glaciers). Retreat rates on land-terminating glaciers were very similar on both coasts: Barents Sea = −6.5 m a −1 , and Kara Sea = −9.0 m a −1 (Fig. 2). For marine-terminating outlets, retreat rates being higher on the Barents Sea confirmed that the significant difference in total retreat rates between land-and marine-terminating glaciers persists when individual coasts are considered (Fig. 2). Finally, we tested for differences in retreat rate between the ice masses of Novaya Zemlya, specifically the northern island ice cap, which is by far the largest, and the two smaller subsidiary ice fields, Sub 1 and Sub 2. Here, we found no significant difference in retreat rates between the ice masses (Fig. 2). Retreat rates were highest on Sub 2, followed by the northern island ice cap, and lowest on Sub 1 (Fig. 2). Our results therefore demonstrate that the only significant difference in total retreat rates  relates to glacier terminus type, with land-terminating outlets retreating 3.5 times slower than those ending in lakes or the ocean (Fig. 2).
We used simple linear regression to assess the relationship between total retreat rate  and latitude, as there is a strong north-south gradient in climatic conditions on NVZ, but no significant linear relationship was apparent (R 2 = 0.001, p = 0.819) (Fig. 3). However, if we divide the glaciers according to terminus type, total retreat rate shows a significant positive relationship for land-terminating glaciers (R 2 = 0.363, p = 0.023), although the R 2 value is comparatively small (Fig. 3). This indicates that more southerly landterminating outlets are retreating more rapidly than those in the north. Conversely, total retreat rate for lake-terminating glaciers has a significant inverse relationship with total retreat rate (R 2 = 0.811, p = 0.014), suggesting that glaciers at high latitudes retreat more rapidly (Fig. 3). No linear relationship is apparent between latitude and total retreat rate for marine-terminating glaciers, and the data show considerable scatter, particularly in the north (Fig. 3). We find no significant relationship between catchment area and total retreat rate (ρ = −0.149, p = 0.339), which demonstrates that observed retreat patterns are not simply a function of glacier size (i.e. that larger glacier retreat more simply because they are bigger).
www.the-cryosphere.net/11/2149/2017/ The Cryosphere, 11, 2149-2174, 2017 Figure 2. Box plots and Kruskal-Wallis test results for different glacier terminus settings for (a) terminus type; (b) coast and terminus: L stands for land-terminating, and m for marine-terminating; and (c) ice mass, specifically the northern island ice cap and subsidiary ice fields 1 and 2. See Fig. 1 for ice mass locations. In all cases, total retreat rate  is used to test for significant differences between the classes. Mean total retreat rates for each class are given on each plot, below the associated box plot. For each box plot, the red central line represents the median; the blue lines represent the upper and lower quartile; red crosses are outliers (a value more than 1.5 times the interquartile range above/below the interquartile values); and the black lines are the whiskers, which extend from the interquartile ranges to the maximum values that are not classed as outliers. P values for each Kruskal-Wallis test are given on the right of the plot.

Temporal change
Based on an initial assessment of the temporal pattern of retreat for individual glaciers, we manually identified major break points in the data and divided glacier retreat rates into four time intervals : 1973/76 to 1986, 1986 to 2000, 2000 to 2013, and 2013 to 2015 (Fig. 4). Data were separated according to terminus type and, in the case of marineterminating glaciers, according to coast. We then used the Wilcoxon test to evaluate the statistical difference between these time periods for each category (Table 2). For land-and lake-terminating glaciers, there were no significant differences in retreat rates between any of the time periods ( Fig. 4; Table 2). Indeed, retreat rates on lake-terminating glaciers were remarkably consistent between 1986 and 2015, both  Table 2). Between 2000 and 2013, retreat rates were much higher than at any other time (−85.4 m a −1 ). Conversely, the average frontal-position change between 2013 and 2015 was positive, giving a mean advance of +11.6 m a −1 (Fig. 4). On the Kara Sea coast, marine-terminating outlet glacier retreat rates were significantly higher between 2000 and 2013 than any other time period (−64.8 m a −1 ) ( Fig. 4; Table 2). Retreat rates reduced substantially during the period 2013-2015 The Cryosphere, 11, 2149-2174, 2017 www.the-cryosphere.net/11/2149/2017/  Table 2. Wilcoxon test results, used to assess significant differences in retreat rates between each manually identified time interval (1976-1986, 1986-2000, 2000-2013, 2013, 2015). Retreat rate data were tested separately for each terminus type, and marine-terminating glaciers were further sub-divided by coast. Following convention, p values of < 0.05 are considered significant and are highlighted in bold.
On both the Barents Sea and Kara Sea coasts, the temporal pattern of marine-terminating outlet glacier retreat showed large variability, both between individual glaciers and over time (Fig. 5). Following our initial analysis, we used change-point analysis to further assess the temporal patterns of glacier retreat, by identifying the timing of significant breaks in the data. On the Barents Sea coast, five glaciers underwent a significant change in retreat rate from the early 1990s onwards (Fig. 6). Of these, retreat rates on four glaciers (MAK, TAI2, VEL, and VIZ; see Fig. 1 for glacier locations and names) subsequently increased, whereas retreat was slower on INO between 1989 and 2006. The most widespread step change on the Barents Sea coast occurred in the early 2000s, after which nine glaciers retreated more rapidly (Fig. 6). A second widespread change in glacier retreat rates occurred in the mid-2000s, which was also the second change-point for four glaciers (Fig. 6). Of these eight glaciers, only VOE retreated more slowly after the mid-2000s change-point. On the Kara Sea coast, we see a broadly similar temporal pattern, with two glaciers showing a significant change in retreat rate from the early 1990s and again in 2005 and 2007 (Fig. 6). In the case of MG, retreat rates were higher after each breakpoint, whereas for SHU1 retreat rates were lower between the 1990s and mid-2000s. Four glaciers began to retreat more rapidly from 2000 onwards, and five other glaciers showed a www.the-cryosphere.net/11/2149/2017/ The Cryosphere, 11, 2149-2174, 2017  1973/76-1986, 1986-2000, 2000-2013, and 2013-2015. (a) Retreat rates were calculated separately for different terminus types, and marine-terminating glaciers were further sub-divided into those terminating into the Barents Sea versus the Kara Sea. Wide bars represent mean values, and thin bars represent the total range (i.e. minimum and maximum values) within each category. significant change in retreat rates beginning between 2005 and 2010 (Fig. 6), with VER being the only glacier to show a reduction in retreat rates after this change (Fig. 6). Focusing on lake-terminating glaciers, a significant change in retreat rates began between 2006 and 2008 on all but one glacier, which began to retreat more rapidly from 2004 onwards (Fig. 6).
In the ERA-Interim reanalysis data, mean annual air temperatures increased significantly between 1986-1999 and 2000-2012 at both the surface and 850 hPa pressure level (Table 3). Winter (surface) and autumn (850 hPa) temperatures also warmed significantly between these time intervals (Table 3). Surface air temperatures were significantly warmer in 2013-2015 than in 1986-1999, in winter and annually (Table 3). No significant differences in air temper-atures were observed at either height between 2000-2012 and 2013-2015 for any season (Table 3). Surface air temperatures were comparable between 2000-2012 and 2013-2015 in winter and autumn, and somewhat warmer in spring (+ 2.6 • C) and summer (+0.7 • C) in 2013-2015 (Fig. 4). At 850 m height, winter (−0.7 • C) and autumn temperatures were slightly cooler (−0.7 • C), and summer temperatures were warmer (+0.8 • C) in 2013-2015 than in 2000-2012 ( Fig. 4). At the regional scale, warmer surface air temperatures penetrate further into the Barents Sea and the southern Kara Sea with each time step (Supplement Fig. S1). We observed a similar, although less marked, northward progression of the isotherms at 850 hPa level (Supplement Fig. S1).
On the Barents Sea coast, sea ice concentrations during all seasons were significantly lower in 2000-2012 than in 1976-1985 or 1986-1999, as was the number of ice-free months ( Fig. 7; Table 4). Between 1976Between -1985Between and 2000Between -2012 winter sea ice concentrations reduced from 68 to 35 %, mean spring values declined from 59 to 28 %, and mean autumn averages fell from 27 to 7 % (Fig. 7). Mean summer sea ice concentrations reduced slightly, from 12 to 5 % (Fig. 7). Over the same time interval, the number of ice-free months increased from 3.0 to 6.9 (Fig. 7). Summer sea ice concentrations on the Barents Sea coast reduced significantly between 2000-2012 and 2013-2015, but no significant change was observed in any other month, nor in the number of ice-free www.the-cryosphere.net/11/2149/2017/ The Cryosphere, 11, 2149-2174, 2017 Table 3. P values for Wilcoxon tests for significant differences in mean seasonal and mean annual air temperatures, for the periods 1976-1985, 1986-1999, 2000-2013, and 2013-2015 Table 4. P values for Wilcoxon tests for significant differences in mean seasonal sea ice concentrations and the number of ice-free months, for the periods 1976-1985, 1986-1999, and 2000-2013 Figure 6. Results of the change-point analysis for glacier retreat rates and climatic controls. Red dots indicate the start of a significantly different period in the time series data, and grey dots represent the end of the previous period, with grey dashed lines connecting the two. This is done to account for missing data: we know that the change-point occurred between the grey and the red dot, and that the new phase of behaviour occurred from the red dot onwards, but not the exact timing of the change. Blue dots show the start of a second significant change in the time series. Frontal-position data were analysed separately for marine-terminating outlets on the Barents Sea coast (a), Kara Sea coast (b), and lake-terminating glaciers (c). (d) Change-point results for seasonal means in air temperatures and sea ice, and the number of ice-free months. Only climatic variables that demonstrated change-points are shown. months ( Fig. 7; Table 4). With exception of winter, sea ice concentrations were significantly lower in 2013-2015 than in 1976-1985or 1986-1999; Table 4). As on the Barents Sea coast, sea ice concentrations on the Kara Sea were significantly lower in all seasons in 2000than in 1976-1985or 1986-1999; Table 4). Summer mean sea ice concentrations declined from 25 % in 1976-1985 to 13 % in 2000-2012 (Fig. 7). Over the same time interval, autumn mean concentrations reduced from 56 to 33 %, spring values declined from 87 to 73 %, and winter values decreased from 87 to 79 % (Fig. 7). The number of ice-free months also reduced from 1. 6 (1976-1985) to 3.0 (2000-2012) (Fig. 7). No significant differences were apparent between seasonal sea ice concentrations and the number of ice-free months in 2013-2015 and any other time period, with the exception of summer sea ice concentrations between 1976-1985 and 2013-2015 (Table 4).
Focusing on the change-point analysis, we see a significant change in air temperatures at E. K. Fedorova from 2008 onwards, after which air temperatures increased markedly (Fig. 6). On the Barents Sea coast, we observe significant breaks in summer sea ice concentrations at 2000 and 2008: www.the-cryosphere.net/11/2149/2017/ The Cryosphere, 11, 2149-2174, 2017  1973/76-1986, 1986-2000, 2000-2013, and 2013-2015.  Fig. S3). The number of ice-free months increased significantly on both the Kara Sea (from 2003) and Barents Sea (from 2005) (Fig. 6). Between 1970 and1989, the summer and annual NAO index were largely positive, with a few years of negative values (Fig. 8a). From 1989 to 1994, values were all positive, followed by strongly negative values in 1995 (Fig. 8a). Subsequently, the summer and annual NAO index remained weakly negative between 1999 and 2012, with values becoming increasingly negative in the final 5 years of this period (Fig. 8a). In 2013, the NAO index became strongly positive, particularly during summer, and values were also positive in 2015 and 2016 (Fig. 8a). The AO index follows an overall similar pattern to the NAO until ∼ 2000, although shifts are less distinct: the index is generally negative until 1988, followed by 5 years of more positive values. In the 2000s, the AO index fluctuates between positive and negative, and more negative summer values are observed in 2009, 2011, 2014, and 2015 (Fig. 8b). The AMO was generally negative from 1970 to 2000, although values fluctuated and were positive around 1990 (Fig. 8c). Subsequently, the AMO entered a positive phase from 2000 onwards (Fig. 8c).
At the broad spatial scale, data indicate that surface ocean temperatures have warmed in the Barents Sea over time (Fig. 9). Warming was particularly marked in the area extending approximately 100 km offshore of the Barents Sea  (Fig. 9), although it should be noted that data are much sparser for the latter period. The Kara Sea also warmed over the study period, with temperatures increasing from 0-2 • C in 1971-1980to 4-5 • C in 2001). Although input data are comparatively sparse for 2001-2012, it appears that ocean temperatures have warmed in both the Barents and Kara seas at each time step, suggesting there may be a broadscale warming trend in the region. At 100 m depth, the data suggest that warmer ocean water extends substantially during the study period, on both the Barents Sea and Kara Sea coasts (Fig. 9).

Glacier surging
During the study period, we observed three glaciers surging: ANU, MAS, and SER (Fig. 1). These were excluded from the analysis of glacier retreat rates and are discussed separately here. ANU has previously been identified as possibly surge type, based on the presence of looped moraine (Grant et al., 2009). Here, we identify an active surge phase, on the basis of a number of characteristics identified from satellite imagery and following the classification of Grant et al. (2009): rapid frontal advance, heavy crevassing, and dig-itate terminus. High flow speeds are also evident close to the terminus (Melkonian et al., 2016), which is consistent with the active phase of surging. Our results show that advance began in 2008 and was ongoing in 2015, with the glacier advancing 683 m during this period (Fig. 10a). MAS was previously confirmed as surge type (Grant et al., 2009), and our data suggest that its active phase persisted between 1989 and 2007 (Fig. 10a). The imagery indicates that surging on MAS originated from the eastern limb of the glacier, which may be partially fed by the neighbouring glacier ( Fig. 10b-f). The exact timing of this tributary surge is uncertain, but imagery from 1985 (Fig. 10c) shows limited evidence of surging, whereas a number of surge indicators are clearly visible by 1988, including looped moraines and rapid advance (Fig. 10d), suggesting it began in the late 1980s. The tributary glacier then advanced into the eastern margin of the main outlet of MAS, causing it to advance, and produced heavy crevassing on the eastern portion of its terminus (Fig. 10d and e). The main terminus of MAS reached its maximum extent for the study period in 2007, and the tributary continued advancing from the 1980s until 2007 (Fig. 10f).  1970-1981, 1981-1990, 1991-2000, and 2001-2012. These intervals were chosen to match as closely as possible with the glacier frontal-position data and other data sets. Note that data coverage was substantially lower for 2001-2012 than for other time periods. Further details on data coverage are available here: https://www.nodc.noaa.gov/OC5/ nordic-seas/.

Spatial patterns of glacier retreat
Our results demonstrate that retreat rates on marineterminating outlet glaciers (−46.9 m a −1 ) were more than 3 times higher than those on land (−13.8 m a −1 ) between 1986 and 2015 (Fig. 2). This is consistent with previous shorter-term studies from Greenland (Moon and Joughin, 2008;Sole et al., 2008) and Svalbard (Dowdeswell et al., 2008), which demonstrated an order-of-magnitude difference between marine-and land-terminating glaciers. It also confirms that the differences in retreat rates, relating to terminus type, observed between 1992 and 2010 on NVZ (Carr et al., 2014) persist at multi-decadal timescales. Recent results suggest that marine-terminating glacier retreat and/or ice tongue collapse can cause dynamic thinning in the RHA (Melkonian et al., 2016;Willis et al., 2015), meaning that these long-term differences in retreat rates may lead to substantially higher thinning rates in marine-terminating basins at multi-decadal timescales. The Russian high Arctic is forecast to be the third-largest source of ice volume loss by 2100 outside of the ice sheets (Radić and Hock, 2011). However, these estimates only account for surface mass balance, not ice dynamics, meaning that they may underestimate 21st-century ice loss for the RHA. Consequently, dynamic changes associated with marine-terminating outlet glacier retreat on NVZ need to be taken into account in order to accurately forecast its near-future ice loss and sea level rise contribution. Our data showed no significant difference in total retreat rates for marine-terminating (−46.9 m a −1 ) and laketerminating glaciers (−49.1 m a −1 ). This contrasts with results from Patagonia, which were obtained during a similar time period (mid-1980s to 2010/11) and showed that lake-terminating outlet glaciers retreated significantly more rapidly than those ending in the ocean (Sakakibara and Sugiyama, 2014). For example, marine-terminating outlets retreat at an average rate of −37.8 m a −1 between 2000 and 2010/11, whereas lake-terminating glaciers receded at −80.8 m a −1 (Sakakibara and Sugiyama, 2014). Laketerminating glacier retreat on NVZ also differs from Patagonia in that retreat rates are remarkably consistent between individual glaciers and remained similar over time (Figs. 4 and 5). Conversely, frontal-position changes in Patagonia showed major spatial variations, and retreat rates on several lake-terminating glaciers changed substantially between the two halves of the study period (mid-1980s-2000 and 2000-2010/11) (Sakakibara and Sugiyama, 2014).
One potential explanation for the common behaviour of the lake-terminating outlet glaciers on NVZ is that retreat may be dynamically controlled and sustained by a series of feedbacks once it has begun. As observed on large Greenlandic tidewater glaciers, initial retreat may bring the terminus close to floatation, leading to faster flow and thinning, which promote further increases in calving and retreat (e.g. Howat et al., 2007;Hughes, 1986;Joughin et al., 2004;Meier and Post, 1987;Nick et al., 2009). This has been suggested as a potential mechanism for the rapid recession for Upsala Glacier in Patagonia (Sakakibara and Sugiyama, 2014) and Yakutat Glacier, Alaska (Trüssel et al., 2013). However, rapid retreat was not observed on all lake-terminating glaciers in Patagonia (Sakakibara and Sugiyama, 2014), and the potential for these feedbacks to develop depends on basal topography (e.g. Carr et al., 2015;Porter et al., 2014;Rignot et al., 2016). Consequently, the basal topography would need to be similar for each of the NVZ glaciers to explain the very similar retreat patterns, which is not implausible but perhaps unlikely. Alternatively, it may be that the proglacial lakes act as a buffer for atmospheric warming, due the greater thermal conductivity of water relative to air, and so reduce variability in retreat rates. Furthermore, lake-terminating glaciers are not subject to variations in sea ice and ocean temperatures, which may account for their more consistent retreat rates, compared to marine-terminating glaciers (Figs. 4 and 5). In order to differentiate between these two explanations, data on lake temperature changes during the study period and lake bathymetry would be required. However, neither are currently available, and we highlight this as an important area for further research, given the rapid recession observed on these lake-terminating glaciers.
For the period between 1986 and 2015, we find no significant difference in retreat rates between the Barents Sea and Kara Sea coasts (Fig. 2). This is contrary to the results of a previous shorter-term study, which showed that retreat rates on the Barents Sea coast were significantly higher than on the Kara Sea between 1992 and 2010 (Carr et al., 2014) and the higher thinning rates observed on marine outlets on the Barents Sea coast (Melkonian et al., 2016). Furthermore, there are substantial differences in climatic and oceanic conditions on the two coasts (Figs. 4 and 7) (Pfirman et al., 1994;Politova et al., 2012;Przybylak and Wyszyński, 2016;Zeeberg and Forman, 2001), so we would expect to see significant differences in outlet glacier retreat rates. This indicates that longer-term glacier retreat rates on NVZ may relate to much broader, regional-scale climatic change, which is supported by the widespread recession of glaciers across the Arctic during the past 2 decades (e.g. Blaszczyk et al., 2009;Carr et al., 2014;Howat and Eddy, 2011;Jensen et al., 2016;Moon and Joughin, 2008). One potential overarching control on NVZ frontal positions is fluctuations in the NAO, which covaries with Northern Hemisphere air temperatures, Arctic sea ice, The Cryosphere, 11, 2149-2174, 2017 www.the-cryosphere.net/11/2149/2017/ and North Atlantic Ocean temperatures (Hurrell, 1995;Hurrell et al., 2003;IPCC, 2013). More recent work has also recognized the influence of the AMO on oceanic and atmospheric conditions in the Barents Sea and broader North Atlantic (Drinkwater et al., 2013;Oziel et al., 2016). Our data suggest that the major phases of frontal-position change on NVZ correspond to changes in the NAO and AMO (Fig. 8; Sect. 4.2): rapid retreat between 2000 and 2013 coincides with a weakly negative NAO and positive AMO, following almost 3 decades characterized by a generally positive NAO and negative AMO (Fig. 8). As such, these large-scale changes may overwhelm smaller-scale spatial variations between the two coasts of NVZ when retreat is considered on multi-decadal timescales. Marine-terminating outlet glacier retreat rates do not show a linear relationship with latitude, and there is considerable scatter when the two variables are regressed (Fig. 3). This may be due to the influence of fjord geometry on glacier response to climatic forcing (Carr et al., 2014) and the capacity for warmer ocean waters to access the calving fronts. In contrast, southerly land-terminating outlets retreat more rapidly than those in the north, which we attribute to the substantial latitudinal air temperature gradient on NVZ (Zeeberg and Forman, 2001). Conversely, lake-terminating glaciers retreat more rapidly at more northerly latitudes (Fig. 3), which we speculate may relate to the bathymetry and basal topography of individual glaciers, but data are not currently available to confirm this.

Temporal patterns
Our results show that retreat rates on marine-terminating outlet glaciers on NVZ were significantly higher between 2000 and 2013 than during the preceding 27 years (Fig. 4). At the same time, land-terminating outlets experienced much lower retreat rates and did not change significantly during the study period (Figs. 4 and 5). This is consistent with studies from elsewhere in the Arctic, which identified the 2000s as a period of elevated retreat on marine-terminating glaciers (e.g. Blaszczyk et al., 2009;Howat and Eddy, 2011;Jensen et al., 2016;Moon and Joughin, 2008) and increasing ice loss (e.g. Gardner et al., 2013;Lenaerts et al., 2013;Moholdt et al., 2012;Nuth et al., 2010;Shepherd et al., 2012). As discussed above, recent evidence suggests that glacier retreat in the Russian high Arctic can trigger substantial dynamic thinning and ice acceleration (Melkonian et al., 2016;Willis et al., 2015), but it not currently incorporated into predictions of 21st-century ice loss from the region (Radić and Hock, 2011). Consequently, the period of higher retreat rates during the 2000s may have a much longer-term impact on ice losses from NVZ, and this needs to be quantified and incorporated into forecasts of ice loss and sea level rise prediction.
Within the decadal patterns of glacier retreat, we observe clusters in the timing of significant changes in marineterminating glacier retreat rates (Fig. 6). Specifically, we see breaks in the frontal-position time series on both the Barents Sea and Kara Sea coasts in the early 1990s, ∼ 2000, and the mid-2000s (Fig. 6). This demonstrates some synchronicity in changes in glacier behaviour around NVZ, although it is not ubiquitous (Fig. 6). The timing of these changes coincides with those observed in Greenland, where the onset of widespread retreat and acceleration in south-east Greenland began in ∼ 2000 (e.g. Howat et al., 2008;Moon and Joughin, 2008;Seale et al., 2011) and occurred from the mid-2000s onwards in the north-west (e.g. Carr et al., 2013b;Howat and Eddy, 2011;Jensen et al., 2016;McFadden et al., 2011;Moon et al., 2012). Whilst these changes could be coincidental, they may also relate to broad, regional-scale changes observed in the North Atlantic region during the 2000s (Beszczynska-Möller et al., 2012;Hanna et al., 2013Hanna et al., , 2012Holliday et al., 2008;Sutherland et al., 2013). Data demonstrate that the NAO was weakly negative from the mid-1990s until 2012, in contrast to strongly positive conditions in the late 1980s and early 1990s, and the AMO was persistently positive from 2000 onwards, following 3 decades of overall positive conditions (Fig. 8). These changes coincide with increases in glacier retreat rates, sea ice decline, and atmospheric warming in NVZ between 2000 and 2013 ( Figs. 4 and 7).
Between the 1950s and mid-1990s, positive phases of the NAO were associated with the influx of warm Atlantic water into the Barents Sea (Hurrell, 1995;Loeng, 1991) and increased penetration of Atlantic cyclones and air masses into the region, which lead to elevated air temperatures and precipitation (Zeeberg and Forman, 2001). Conversely, negative NAO phases were associated with cooler oceanic and atmospheric conditions in the Barents Sea (Zeeberg and Forman, 2001). During this period, therefore, the impact of the NAO was opposite in the Barents Sea and in western portions of the Atlantic-influenced Arctic (e.g. the Labrador Sea) (Drinkwater et al., 2013;Oziel et al., 2016). However, since the mid-1990s, changes in the Barents Sea and the western Atlantic Arctic have been in phase, and warming and sea ice reductions have been widespread across both regions (Drinkwater et al., 2013;Oziel et al., 2016). As such, increased glacier retreat rates on NVZ during the 2000s (Figs. 4 and 5) may have resulted from the switch to a weaker, and predominantly negative, NAO phase from the mid-1990s (Fig. 8), which would promote warmer air and ocean temperatures, and reduced sea ice, as we observe in our data (Figs. 4 and 7). Previous studies have suggested a 3-5-year lag between NAO shifts and changes in conditions on NVZ, due to the time required for Atlantic water to transit into the Barents Sea (Belkin et al., 1998;Zeeberg and Forman, 2001), which is consistent with the onset of retreat in ∼ 2000 (Figs. 4 and 8). However, it has recently been suggested that the NAO's role may have reduced since the mid-1990s and that the AMO may be the dominant influence on warming in the North Atlantic (Drinkwater et al., 2013;Oziel et al., 2016). The AMO is thought to promote blocking of highwww.the-cryosphere.net/11/2149/2017/ The Cryosphere, 11, 2149-2174, 2017 pressure systems by westerly winds, which changes the wind field (Häkkinen et al., 2011). This allows warm water to penetrate further into the Barents Sea and other Nordic Seas, leading to atmospheric and oceanic warming during periods with a weakly negative NAO (Häkkinen et al., 2011). As such, rapid retreat on NVZ between 2000 and 2013 may have resulted from the combined effects of a weaker, more negative NAO from the mid-1990s and a more positive AMO from 2000 onwards (Fig. 8). This suggests that synoptic climatic patterns may be an important control on glacier retreat rates on NVZ and that the recent relationship between the NAO and glacier change on NVZ contrasts with that observed during the 20th century (Zeeberg and Forman, 2001). Following higher retreat rates in the 2000s, our data indicate that marine-terminating glacier retreat slowed from 2013 onwards on both the Barents Sea and Kara Sea coasts, with several glaciers beginning to re-advance (Figs. 4 and 5). Our data demonstrate that marine-terminating glaciers on NVZ have previously undergone a step-like pattern of retreat, with short (1-2 year) pauses in retreat (Fig. 5). Thus, it is unclear whether this reduction in retreat rates is another temporary pause, before continued retreat, or the beginning of a new phase of reduced retreat rates. One possible explanation for reduced retreat rates on both coasts of NVZ is the stronger NAO values observed from the late 2000s onwards: winter 2009/10 had the most negative NAO for 200 years (Delworth et al., 2016;Osborn, 2011), and values were strongly positive in 2013 (Fig. 8a). This is consistent with the 3-5-year lag required for NAO-related changes in Atlantic water inflow to reach NVZ (Zeeberg and Forman, 2001), and so we speculate that reduced glacier retreat rates from 2013 onwards (Figs. 4 and 5) may relate to an increase in the influence of the NAO, relative to the AMO, from the late 2000s (Fig. 8). Evidence indicates that the impact of the NAO in the Barents Sea is now in phase with the western North Atlantic (Drinkwater et al., 2013;Oziel et al., 2016), and so a more positive NAO could lead to cooler conditions on NVZ and hence glacier advance. However, the relationship between large-scale features, such as the NAO and AMO; ocean conditions; and glacier behaviour is complex (Drinkwater et al., 2013;Oziel et al., 2016), and the period of glacier advance/reduced retreat on NVZ has lasted only 2 years. Consequently, further monitoring is required to determine whether this represents a longer-term trend or a short-term change and to confirm its relationship to synoptic climatic patterns.
Despite the changes in the NAO and AMO, our data show no significant change in sea ice concentrations, nor the length of the ice-free season, between 2000-2012 and 2013-2015 on either the Barents Sea or Kara Sea coast (Table 4; Fig. 7). Likewise, we see no significant change in winter (January-March) air temperatures at E. K. Fedorova (Table 3; Fig. 4) nor in the ERA-Interim data during any season (Table 3; Fig. 4). Although not significant, we see summer warming of 0.7 • C (surface) and 0.8 • C (850 hPa pressure level) in the ERA-Interim data (Fig. 4), which is the opposite of what we would expect if reductions in air temperatures and surface melt were driving the slowdown in retreat rates. As such, reduced retreat rates do not seem to be directly linked to short-term changes in sea ice or air temperatures. They are also unlikely to result from changes in surface mass balance, as the response time of NVZ glaciers is likely to be slow: they have long catchments (∼ 40 km) and slow flow speeds (predominantly < 200 m a −1 ; Melkonian et al., 2016) and are likely to be polythermal. Furthermore, thinning rates between 2012 and 2013/14 averaged 0.4 m a −1 across the ice cap and reached up to 5 m a −1 close to the glacier termini (Melkonian et al., 2016), meaning that even a positive surface mass balance is very unlikely to deliver sufficient ice quickly enough to promote advance and/or substantially lower retreat rates. Instead, this may be a response to oceanic changes, which we cannot detect from available data; it may reflect a lagged response to forcing; and/or it may relate to more localized, glacier-specific factors. We suggest that the latter is unlikely, given the widespread and synchronous nature of the observed reduction in retreat rates (Figs. 4 and 5). Future work should monitor retreat rates to determine whether reduced retreat is persistent or is a short-term interruption to overall glacier retreat and collect more extensive oceanic data to assess its impact on this change. Furthermore, detailed data are also required to determine how short-term frontalposition fluctuations relate to changes in ice velocities and/or surface elevation.
Although we observe some common behaviour, in terms of the approximate timing and general trend in retreat, there is still substantial variability in the magnitude of retreat between individual marine-terminating glaciers (Figs. 4 and 5). Furthermore, not all glaciers shared common change-points, and certain outlets showed a different temporal pattern of retreat to the majority of the study population . For example, INO retreated more slowly between 1989 and 2006 than during the 1970s and 1980s. We attribute these differences to glacier-specific factors and, in particular, the fjord bathymetry and basal topography of individual glaciers. Previous studies have highlighted the impact of fjord width on retreat rates on NVZ (Carr et al., 2014) and basal topography on marine-terminating glacier behaviour elsewhere (e.g. Carr et al., 2015;Porter et al., 2014;Rignot et al., 2016). This may result from the influence of fjord geometry on the stresses acting on the glacier once it begins to retreat: as a fjord widens, lateral resistive stresses will reduce, and the ice must thin to conserve mass, making it more vulnerable to calving (Echelmeyer et al., 1994;Raymond, 1996;van der Veen, 1998a, b), whilst retreat into progressively deeper water can cause feedbacks to develop between thinning, floatation, and retreat (e.g. Joughin and Alley, 2011;Joughin et al., 2008;Schoof, 2007). Thus, retreat into a deeper and/or wider fjord may promote higher retreat rates on a given glacier, even with common climatic forcing. In addition, differences in fjord bathymetry may determine whether warmer Atlantic water can access the glacier front (Porter et al., 2014;Rignot The Cryosphere, 11, 2149-2174 www.the-cryosphere.net/11/2149/2017/ et al., 2016), which could promote further variations between glaciers. This highlights the need to collect basal topographic data for NVZ outlet glaciers, which is currently very limited but a potentially key control on ice loss rates.

Climatic and oceanic controls
Our data demonstrate that air temperatures were very substantially warmer between 2000 and 2012 than during the preceding decades and that sea ice concentrations were also much lower on both the Barents Sea and Kara Sea coasts during this period (Figs. 4 and 7). This is consistent with the atmospheric warming reported across the Arctic during the 2000s (e.g. Carr et al., 2013a;Hanna et al., 2013Hanna et al., , 2012Mernild et al., 2013) and the well-documented decline in Arctic sea ice (Comiso et al., 2008;Kwok and Rothrock, 2009;Park et al., 2015). As such, the decadal patterns of marine-terminating outlet glacier retreat correspond to decadal-scale climatic change on NVZ (Figs. 4 and 7), and exceptional retreat during the 2000s coincided with significantly warmer air temperatures and lower sea ice concentrations (Tables 2 and 3). Interestingly, step changes in the air temperature and sea ice data identified by the change-point analysis did not correspond to significant changes in outlet glacier retreat rates (Fig. 6), suggesting that such changes may not substantially influence retreat rates or that the relationship may be more complex, e.g. due to lags in glacier response.
The much lower retreat rates on land-terminating outlets (Fig. 4) may indicate an oceanic driver for retreat rates on marine-terminating glaciers. Previous studies have identified sea ice loss as a potentially important control on NVZ retreat rates (Carr et al., 2014), which fits with observed correspondence between sea ice loss and retreat, but it is unclear whether the two variables simply co-vary or whether sea ice can drive ice loss, by extending the duration of seasonally high calving rates (e.g. Amundson et al., 2010;Miles et al., 2013;Moon et al., 2015). The available ocean data indicate that temperatures were substantially warmer during the 2000s (Fig. 9), which would provide a plausible mechanism for widespread retreat on both coasts of NVZ (Fig. 4). However, oceanic data for the 2000s is sparse in the Barents and Kara seas, compared to previous decades, so it is difficult to ascertain the magnitude and spatial distribution of warming and to link it directly with glacier retreat patterns. Laketerminating glaciers are not affected by changes in sea ice or ocean temperatures but could be influenced by air temperatures. However, despite much higher air temperatures in the 2000s, mean retreat rates on lake-terminating outlet glaciers were similar for each decade of the study (Fig. 4), suggesting that the relationship is not straightforward. Instead, the presence of lakes may at least partly disconnect these glaciers from climatic forcing, by buffering the effects of air temperatures changes and/or by sustaining dynamic changes, follow-ing initial retreat (Sakakibara and Sugiyama, 2014;Trüssel et al., 2013).

Glacier surging
During the study period, we identify three actively surging glaciers, based on various lines of glaciological and geomorphological evidence (Copland et al., 2003;Grant et al., 2009), including terminus advance (Fig. 10). Frontal advance persisted for 18 years on MAS and 15 years on SER, whilst ANU began to advance in 2008, and this continued until the end of the study period (Fig. 10a). This is comparatively long for surge-type glaciers, which usually undergo short active phases over time frames of months to years (Dowdeswell et al., 1991;Raymond, 1987). For comparison, surges on Tunabreen, Spitzbergen, last only ∼ 2 years , and Basin 3 on Austfonna underwent major changes in its dynamic behaviour in just a few years (Dunse et al., 2015). Surges elsewhere can occur even more rapidly: the entire surge cycle of Variegated Glacier in Alaska takes approximately 1-2 decades, and the active phase persists for only a few months (e.g. Bindschadler et al., 1977;Eisen et al., 2005;Kamb, 1987;Kamb et al., 1985;Raymond, 1987). Furthermore, the magnitude of advance on these three glaciers is in the order of a few hundred metres, which is smaller than advances associated with surges on Tunabreen (1.4 km) and Kongsvegen (2 km)  and much less than the many kilometres of advance observed on Alaskan surgetype glaciers, such as Variegated Glacier (Bindschadler et al., 1977;Eisen et al., 2005). Consequently, the active phase on NVZ appears to be long, in comparison to other regions, and terminus advance is more limited, which may provide insight into the mechanism(s) driving surging here and may indicate that these glaciers are located towards one end of the climatic envelope required for surging in the Arctic .
During the active phase of the NVZ surge glaciers, we observe large sediment plumes emanating from the glacier terminus (Fig. 10g), which indicates that at least part of the glacier bed is warm-based during the surge. Together with the comparatively long surge interval, this supports the idea that changes in thermal regime may drive glacier surging on NVZ, as hypothesized for certain Svalbard glaciers (Dunse et al., 2015;Murray et al., 2003;. In addition, the surge of MAS appears to have been triggered by a tributary glacier surging into its lateral margin (Fig. 10bf). This demonstrates an alternative mechanism for surging, aside from changes in the thermal regime and/or hydrology conditions of the glacier, which has not been widely observed but will depend strongly on the local glaciological and topographical setting of the glacier. The data presented here focus only on frontal advance and glaciological/geomorphological evidence, whereas information on ice velocities is also an important indicator of surging . Consequently, information on velocity and surface elevation www.the-cryosphere.net/11/2149/2017/ The Cryosphere, 11, 2149-2174, 2017 changes are needed to further investigate the surge cycle and its possible controls on NVZ. This is important, as NVZ is thought to have conditions that are highly conducive to glacier surging  but has a long surge interval. We therefore want to ensure that we can disentangle surge behaviour and the impacts of climate change on NVZ.

Conclusions
At multi-decadal timescales, terminus type remains a major overarching determinant of outlet glacier retreat rates on NVZ. As observed elsewhere in the Arctic, land-terminating outlets retreated far more slowly than those ending in the ocean. However, we see no significant difference in retreat rates between ocean-and lake-terminating glaciers, which contrasts with findings in Patagonia. Retreat rates on laketerminating glaciers were remarkably consistent between glaciers and over time, which may result from the buffering effect of lake temperature and/or the impact of lake bathymetry, which could facilitate rapid retreat that is largely independent of climate forcing, after an initial trigger. We cannot differentiate between these two scenarios with currently available data. Retreat rates on marine-terminating glaciers were exceptional between 2000 and 2013, compared to previous decades. However, retreat slowed on the vast majority of ocean-terminating glaciers from 2013 onwards, and several glaciers advanced, particularly on the Barents Sea coast. It is unclear whether this represents a temporary pause or a longer-term change, but it should be monitored in the future, given the potential for outlet glaciers to drive dynamic ice loss from NVZ. The onset of higher retreat rates coincides with a more negative, weaker phase of the NAO and a more positive AMO, whilst reduced retreat rates follow stronger NAO years. This suggests that synoptic atmospheric and oceanic patterns may influence NVZ glacier behaviour at decadal timescales. Marine-terminating glaciers showed some common patterns in terms of the onset of rapid retreat (1990s, ∼ 2000 and mid 2000s) but showed substantial variation in the magnitude of retreat, which we attribute to glacierspecific factors. Glacier retreat corresponded with decadalscale climate patterns: between 2000 and 2013, air temperatures were significantly warmer than the previous decades and sea ice concentrations were significantly lower. Available data indicate oceanic warming, which could potentially explain why retreat rates on marine-terminating glaciers far exceed those ending on land, but data are comparatively sparse from 2000 onwards, making their relationship to glacier retreat rate difficult to evaluate. The surge phase on NVZ glaciers appears to be comparatively long and warrants further investigation to separate its impact on ice dynamics from that of climate-induced change and to determine the potential mechanism(s) driving these long surges. Recent results suggest that outlet glaciers can trigger dynamic losses on NVZ, but these processes are not yet included in estimates of the region's contribution to sea level rise. As such, it is vital to determine the longer-term impacts of exceptional glacier retreat during the 2000s and to monitor the near-future behaviour of these outlets.
Data availability. The primary dataset created by the paper is the glacier frontal position data, which are provided in the Supplement. Shapefiles of the data can be provided on request to the lead author. Other datasets (e.g. climate data) are available online, and the sources are given in the paper.