Mass Balance of Novaya Zemlya Archipelago, Russian High Arctic, Using Time-Variable Gravity from GRACE and Altimetry Data from ICESat and CryoSat-2

We examine the mass balance of the glaciers in the Novaya Zemlya Archipelago, located in the Russian High Arctic using time series of time-variable gravity from the NASA/DLR Gravity Recovery and Climate Experiment (GRACE) mission, laser altimetry data from the NASA Ice Cloud and land Elevation Satellite (ICESat) mission, and radar altimetry data from the European Space Agency (ESA) CryoSat-2 mission. We present a new algorithm for detecting changes in glacier elevation from these satellite altimetry data and evaluate its performance in the case of Novaya Zemlya by comparing the results with GRACE. We find that the mass loss of Novaya Zemlya glaciers increased from 10 ± 5 Gt/year over 2003–2009 to 14 ± 4 Gt/year over 2010–2016, with a brief period of near-zero mass balance between 2009 and 2011. The results are consistent across the gravimetric and altimetric methods. Furthermore, the analysis of elevation change from CryoSat-2 indicates that the mass loss occurs at elevation below 700 m, where the highest thinning rates are found. We also find that marine-terminating glaciers in Novaya Zemlya are thinning significantly faster than land-terminating glaciers, which indicates an important role of ice dynamics of marine-terminating glaciers. We posit that the glacier changes have been caused by changes in atmospheric and ocean temperatures. We find that the increase in mass loss after 2010 is associated with a warming in air temperatures, which increased the surface melt rates. There is not enough information on the ocean temperature at the front of the glaciers to conclude on the role of the ocean, but we posit that the temperature of subsurface ocean waters must have increased during the observation period.


Introduction
In the twentieth century, with a marked increase in global atmospheric temperatures at the end of the Little Ice Age, the mass balance of mountain glaciers and ice caps (GIC) has been predominantly negative [1][2][3][4]. Available observations suggest that the shrinkage process continued during the first decades of the 21st century [4,5]. Recently, Reager et al. [6] estimated that glaciers outside Greenland and Antarctica have been losing mass at a rate of 253 ± 35 Gt/year between 2002 and 2014. This result places GIC as a major contributor to present and near-term eustatic sea-level rise (SLR). More than 70% of the total glacier mass loss takes place in the Arctic [7], where albedo feedback caused by decreasing sea ice and snow cover yields warming rates 2-3 times larger than the global average [8][9][10]. This study

Study Region
The RHA includes 51,592 km 2 of GIC, which is 13% of the total glaciated area outside of Greenland and Antarctica [23]. The glaciated areas of the RHA are distributed across four regions: Novaya Zemlya (22,379 km 2 ), Severnaya Zemlya (16,382 km 2 ), Franz Josef Land (12,756 km 2 ), and Ushakon Island (359 km 2 ). NZEM, lies between the Barents and Kara Seas and consists of two major islands, Severny Island and Yuzhny Island, along with several smaller islands. Severny and Yuzhny Islands are separated by the Matochkin Strait, a 600 m wide channel that stretches 100 km and connects the Barents and Kara Seas. An axial mountain range extends over the entire length of the NZEM with a maximum elevation of 1340 m on Yuzhny Island and 1596 m on Severny Island [22]. 92% of the total glacier area (20,784.4 km 2 ) is located on Severny Island and concentrated in a single large ice cap (Northern Icefield) that covers 45% of the island and constitutes the largest body of ice in the Eurasian Continent [23,24] (Figure 1).
The Northern Icefield is more than 400 km long, with a mean east-west width of 80 km [24]. The ice cap has an average elevation of 800 m and an average thickness of 400-450 m [25]. 71% of the Northern Icefield area (16,064 km 2 ) is drained by marine-terminating glaciers while most of the glaciers separated from the main ice cap terminate on land. Glaciers were classified as cold-based [26]; but recent studies suggest a wider range of thermal regime and a transition to warm-based/polythermal glaciers [22,27]. A total of 32 glaciers, predominantly located along the western margin of the Northern Icefield, has been identified as potential surge-type [27].
The climate in the region is determined by a complex interaction between multiple oceanic and atmospheric forcings. On the North-West, mild temperatures are favored by the advection of warm and salty water transported by the North Cape current from the Atlantic Ocean [28][29][30]. Annual precipitation (400 mm/year) is favored by the influx of air masses rich in moisture by the Atlantic Cyclone [28,30]. Atmospheric conditions become gradually drier and colder toward the Southeast. The central mountain chain provides an orographic barrier for the eastward penetration of the Atlantic cyclonic system [28] and precipitation gradually decreases toward the South. On the Kara Sea coast, annual precipitation amounts to 250 mm/year. Here, lower atmospheric temperatures, especially in Winter and Spring, are controlled by semi-permanent sea ice and cold Arctic water [28]. January and February are the coldest months of the year with minimum air temperatures decreasing from the North-West (−15 • C) to the South-East (−21 • C). Maximum temperature (+2 • C) is recorded in August, with no significant difference between the two coasts. Maximum precipitation is between September and October. April and May are the driest months of the year [28,31].

Prior Studies
NZEM's present-day average temperatures are 3 • C higher than in the 19th century [31]. Atmospheric warming throughout the first decades of the 20th century contributed to glacier shrinkage [32], during which glaciers in the region lost mass with an average rate between 0.16 and 0.30 m of water equivalent per year [3,13,22,33,34]. The process was not uniform in time, but consisted of a warm first phase between the 1920s and the 1960s with extremely fast terminus retreat rates (greater than 300 m/year), followed by a second phase between the 1960s and the 1990s characterized by less negative mass balance and lower retreat rates (50-150 m/year) [22,32]. Zeeberg and Forman [32] observed that the decadal variability in glacier mass balance on NZEM is linked to long-term shifting of atmospheric and oceanic circulation patterns caused by the North Atlantic Oscillation (NAO). During a positive phase of the NAO, sea-level pressure at high latitudes is below normal, while above normal conditions are registered at mid-latitudes. An increased pressure gradient over the North Atlantic contributes to an intensified and more zonally oriented jet stream. When the NAO is in a negative phase, a decreased pressure gradient between high and mid-latitudes weakens the jet stream which becomes more meridionally oriented [35]. When the NAO is in a prolonged positive phase, the straightening of the polar vortex contributes to an increased flux of air masses from the Atlantic Ocean, resulting in elevated winter precipitation over NZEM. A stronger North Atlantic pressure gradient intensifies the advection of warm water from the Atlantic into the Barents Sea [32]. The Atlantic Water (AW) releases heat into the atmosphere and yields warmer than normal winter and summer temperatures. During extended negative NAO phases, the jet stream weakens and atmospheric conditions in the region become drier and colder. Between the 1920's and 1950's, the NAO was predominantly positive and higher than normal atmospheric and ocean temperatures favored fast glacial recession. After 1960, the NAO turned to negative, i.e., colder atmospheric conditions and lower ocean temperatures, which slowed down the glacial retreat [32]. Several studies have shown that the glaciers of NZEM have been experiencing a state of prolonged negative mass balance [6,11,22,36]. Moholdt et al. [11] used satellite laser altimetry from ICESat and gravity data from GRACE between October 2003 and October 2009 to estimate a glacier loss of 7.6 ± 1.2 to 5.8 ± 3.0 Gt/year, with rapid glacier thinning at elevations below 500 m (0.9 m/year). The analysis of elevation change, however, revealed no significant difference between frontal thinning of marine-and land-terminating glaciers (0.94 m vs. 0.89 m), but significant differences between glaciers on the North-West (Barents Sea) (0.45 m/year) vs. South-East (Kara Sea) (0.25 m/year), consistent with SMB dominance observed in other Arctic regions [16][17][18][19]. Zhao et al. [37] used multiple active and passive microwave data between 1995 and 2011 to report a progressive increase in melt season length, consistent with increased atmospheric temperatures. Increased snowmelt and summer runoff contributed to the high ice thinning rates observed by Moholdt et al. [11]. Yet, new studies suggest that increased ice discharge from marine-terminating glaciers may be significant [13,14,36,38]. Carr et al. [14] documented NZEM glacier retreat between 1973 and 2015 for both marine-and land-terminating glaciers. Marine-terminating glaciers retreated 3.5 faster than land-terminating glaciers (46.9 vs. 13.8 m/year) and faster in the Barents Sea vs. Kara Sea (55.9 vs. 37.2 m/year). Glacier recession accelerated after year 2000 on the Barents Sea and after 2003 on the Kara Sea [36,38]. Marine glaciers retreat rates peaked at 85.4 m/year for the Barents Sea and 64.8 m/year for the Kara Sea. These retreats coincide with warmer atmospheric and ocean temperatures and lower sea ice concentrations compared to the previous decades. These results were evaluated by Melkonian et al. [13] who mapped ice velocity between 1999 and 2014 and elevation change between 1952 and 2013/2014 to find the highest ice speeds and thinning rates at the front of marine-terminating glaciers facing the Barents Sea. Melkonian et al. [13] documented an acceleration in ice velocity and thinning rates after year 2000. Steady increases in marine-terminating glacier velocity were observed by Sun et al. [39] between 2002 and 2015 and by Strozzi et al. [15] between 2008 and 2016. All these findings indicate an increasingly negative regional glacier mass balance after the year 2000.

Time-Variable Gravity
The Gravity Recovery and Climate Experiment (GRACE) is a joint operation by NASA and the German Aerospace Center (DLR) to measure changes in the Earth gravity field at a monthly time scale [40]. Here, we use 156 monthly GRACE Release-5 (RL05) gravity field solutions provided by the Center for Space Research at the University of Texas (CSR) for the period April 2002-August 2016 [40,41]. Monthly gravity field anomalies from CSR are distributed in the form of fully normalized spherical harmonics up to degree and order 60. We replace the C 20 coefficients with more accurate estimates from the satellite laser ranging missions [42]. We account for the variation of degree-1 using coefficients calculated from a combination of GRACE coefficients and ocean model outputs following Swenson et al. [43]. We include an additional pole tide correction to remove the long-period pole tide signals not included in the standard corrections [44]. We remove the effects from Glacial Isostatic Adjustment (GIA), i.e., the viscoelastic response by the solid Earth from changes in ice loading since the last glacial maximum (LGM), using coefficients from A et al. [45]. To isolate gravity changes related to ice mass change, we subtract the signal contribution generated by changes in terrestrial hydrology from the monthly terrestrial water storage estimates from GRACE (i.e., groundwater, soil moisture and snow cover). We estimate these components using the average outputs from two land surface models: (1) the Community Land Surface Model (CLM) version 4.5 [46] and (2) the Global Land Data Assimilation System 2 (GLDAS-2) model, version NOAH-3.3 [47].
We calculate time series of glacier mass change using the least squares fit mascon approach presented in Jacob et al. [12]. We cover all regions with a glacier area >100 km 2 with one or more mascons. For each mascon we calculate a set of Stokes coefficients associated with a unitary mass equivalent to 1 cm of water uniformly distributed over the mascon surface [48]. We calculate the monthly mass anomaly associated with each mascon by simultaneously fitting the mascons Stokes coefficients to the GRACE monthly coefficients corrected for GIA and the Terrestrial Hydrology [12].
We fit the regional time series to an annual, semiannual, linear trend, and constant components in order to investigate the presence of long-term trends in the glacier mass balance. We add a quadratic term to assess potential acceleration in mass loss. The total uncertainty in mass loss rate is the sum in quadrature the GRACE of measurement errors, error in GIA correction, error in hydrology correction, leakage error caused by the assumption of uniform mass change within each mascon, ocean/atmospheric mass correction errors, and the statistical uncertainty of the model fit. To compare the quadratic and linear models, we use a variant of the Akaike Information Criterion, AIC c , with small sample size datasets [49] that accounts for the goodness of the fit and number of parameters, and identifies which model best fits the data. We calculate GRACE measurement errors following Wahr et al. [50]. GIA and hydrology uncertainties are calculated as in Jacob et al. [12] and Gardner et al. [51]. To evaluate the leakage error on the trend introduced by assuming a uniform mass distribution within each mascon, we use a Monte Carlo approach. We generate 10,000 synthetic configurations by distributing the monthly mass changes, (from the mascon inversion described above), non-uniformly within each mascon. The 10,000 random configurations are generated using pseudo-random weights extracted from a Gaussian distribution with standard deviation equal to 1 and zero mean. At each iteration, we apply the mascon fit to the synthetic data and calculate the regional leakage error as the difference between the trend obtained applying the inversion to GRACE and the trend obtained using the randomly generated synthetic data. Our final estimate of the uncertainty associated with the leakage error is set equal to 2 times the standard deviation of all errors measured in the 10,000 iterations. The leakage error associated with our mascon configuration is taken at the 2σ error. Trends uncertainties due to the ocean mass correction are calculated following Velicogna and Wahr [52].
We investigate the interannual variability of the glacier mass anomaly by smoothing the regional time series with a 13-month moving average filter. We define a glaciological year as the time period separating two consecutive ablation seasons, i.e., the 1st September in year n and 31st August in year n + 1. We calculate the regional glacier annual mass balance as the difference between values from the smoothed time series in two consecutive September months.

ICESat Altimetry
The NASA Ice, Cloud, and Land Elevation Satellite (ICESat) was an Earth orbiting laser altimeter designed to measure elevation changes of ice sheets and arctic ice caps mission operative between February 2003 and October 2009. The main instrument carried by the satellite was the Geoscience Laser Altimeter System (GLAS), a single-beam laser altimeter operating at a wavelength of 1064 nm [53]. The satellite flew at an elevation of 600 km with an orbit inclination of 94 degrees. GLAS took measurements at 40-Hz, therefore sampling the Earth's surface elevation every 172 m along the satellite flight direction with a footprint size of almost 70 m [54]. Here, we use level-2 Global Land Surface Altimetry HDF5 (GLAH14) release 34, which include quality attributes and elevation corrections for each footprint. The data quality attributes include a wave saturation flag to indicate saturation of the sensor when recording the returned pulse and a correction of a potential bias in the extracted elevation [55,56]. We apply the saturation and the intercampaign bias corrections. The Gaussian-Centroid correction has been already applied to this data release [57,58]. We use glacier outlines from the 6th release of the Randolph Glacier Inventory [23] to isolate elevation measurements over the NZEM (For more details, see Appendix A). Following Moholdt et al. [11] and Gardner et al. [16], we do not apply the data parameter-based cloud filtering [59], as it tends to be too strict and to discard a high number of elevation measurements that could be otherwise used to measure elevation changes. Over Novaya Zemlya, the application of the data-culling procedure reduces the number of available measurements by 65%.
We evaluate glacier elevation changes employing an along-track plane-fit technique following Moholdt et al. [60]. This approach evaluates ice elevation change rates and terrain slope simultaneously by least squares fitting a time-variable plane function to a set of measurements within a spatial range of 700 m. The plane function is presented in Equation (1) where x, y, t represent easting, northing and time, respectively.
The plane fit is applied at multiple locations considering planes with centroids equally spaced along each satellite track. The spacing is chosen equal to one half the size of the plane radius (350 m), consecutive planes are therefore partially overlapping [11]. To reduce the effect of erroneous elevation measurements, the least squares inversion is applied iteratively at each centroid location discarding all measurements with a residual value larger than 3 times the standard deviation of all residuals (3-σ Clip). This procedure is repeated iteratively until all observations have a residual contained within the 3-σ interval [61,62]. In addition, we consider only measurements containing at least 15 observations, from more than 4 different sub-tracks, and distributed on a temporal interval of at least 2 years. We reject elevation changes derived from planes with an estimated slope greater than 10 • .

CryoSat-2
The ESA CryoSat-2 satellite system was launched in April 2010 to study changes affecting the arctic regions [63]. CryoSat-2 payload consists of a highly sophisticated radar altimeter suitable to study land and sea ice. The use of radar technology allows the altimeter to work in all weather conditions with no limitations related to cloud coverage, which is persistent in the Arctic and highly affects laser altimeters. CryoSat-2 flies at an elevation of 717 km with an orbit inclination equal to 92 • that allows an almost unprecedented coverage of polar regions. The payload of the CryoSat-2 mission consists in the SAR/Interferometric Radar Altimeter (SIRAL). SIRAL is classified as a pulse limited delay-doppler altimeter (working in the Ku band-13.573 GHz, wavelength ∼2.2 cm) with a footprint size of 1.65 × 1.65 km. Instead of transmitting a single pulse per time, SIRAL transmits a burst of 64 pulses every 50 microsecond. The returning echoes are affected by small shift in phase due to the movement of the satellite during the interval between the burst transmission and acquisition (delay-doppler effect) [64]. The phase shift between the different echoes is used to increase the altimeter resolution along the direction of flight of the satellite (azimuth direction) up to a maximum of 400 m on a flat terrain [65]. The small size of the radar footprint makes possible for the first time the application of satellite radar altimetry in regions such as the ice sheet margins and small arctic ice caps where the presence of steep slopes and discontinuous surfaces limited the application of previous microwave altimeters [63].
Over NZEM, SIRAL operates in SAR Interferometry mode (SARIn). In this configuration, a secondary antenna that receives each backscattered signal from a slightly different position with respect to the primary one is activated. The signals received by the two antennas are characterized by a shift in phase due to the different paths linking them to the ground scatterer. Analyzing the difference in phase delay between the two signals, SIRAL detects not only the range distance between the satellite and the ground scatterer but also its location in the across-track plane [63,66,67]. Major limitations to the application of CryoSat-2 to estimate elevation changes over GIC are associated with a radar footprint one order of magnitude larger than for a LiDAR system. Radar performance significantly deteriorates over terrains with slopes higher than 1 • . However, Wang et al. [68] demonstrated that even in these operative conditions, the SIRAL accuracy is still one order of magnitude higher than the accuracy of prior radar altimeters e.g., ERS-1/2 and Envisat. A second limitation of CryoSat-2 is related to the temporal variability in scattering properties of ice and snow. As reported in [69][70][71][72] changes in snowpack reflection introduce episodic artifacts that influence the estimation of elevation change especially over short time periods. Elevation measurements are not aligned along the flight direction but scattered around it. After multiple passages over the same region, the elevation measurements are distributed over a large area, complicating the application of the classic cross-over or co-linear analysis used with other altimeters. The non-uniform distribution of the returning echoes introduces biases in the derivation of ice elevation changes. Elevation measurements are mainly at higher elevations, while depressions and low elevation regions are undersampled [69].
In this study, we employ Level-2 Baseline-C SARin mode elevation data provided by the European Space Agency [73] and freely available at (ftp://science-pds.cryosat.esa.int). We locate elevation measurements over ice using the recently released glacier inventory by Rastner et al. [24] (For more details, see Appendix A). We map elevation changes over the glaciers of NZEM using the plane fit employing all available elevation measurement locations as plane centroids as in Wouters et al. [74]. Given a single centroid, we locate all elevation data within 1000 m from the plane centroid. We least-squares fit the time-variable plane function of the elevation observations. We repeat the fit iteratively until all elevation measurements are within the 3-σ interval. Similar to our ICESat analysis, we least-squares fit the time-variable plane function to the elevation observations from CryoSat-2. Given the larger glacier area considered in the CryoSat-2 case, we repeat the inversion using two more model functions: We obtain our final estimate selecting the elevation trend value associated with the model providing the best fit to the considered elevation measurements. We use the R 2 adj as a model selection criterion (i.e., we chose the dh/dt value associated with the model function characterized by the largest R 2 adj ). We use the plane-fit residual to calculate regional monthly time series of ice elevation change following Wouters et al. [74] and Noël et al. [75] (For more details, see Appendix B).

Averaging Procedure
We average the elevation change measurements obtained with ICESat and in CryoSat-2 on a 1 km grid. The grid is defined on the standard NSIDC/North Polar Stereographic Projection (EPSG: 3413) (https://nsidc.org/data/polar-stereo/ps_grids.html). The averaging procedure is applied as in Wouters et al. [74] following the steps presented below: • If the considered grid cell contains a single measurement, we assign this value to the entire grid cell; • If the considered grid cell contains multiple measurements, we set its elevation change rate equal to the measurements mean. We use a 3-σ iterative filter in order to discard erroneous observations that could influence our final elevation change estimate.

Spatial Extrapolation
The evaluation of the total regional ice volume change requires the extrapolation of elevation change rates at the locations where measurements are not available. In this study, we update the approach presented in Moholdt et al. [11] by designing an extrapolation scheme that considers not only the relation between elevation change and absolute elevation, but also its spatial variability. We extrapolate elevation change rates for grid cells with no assigned elevation change value by applying the following steps. Given a grid cell with no value: • we find all the measurements available within a distance smaller than or equal to a selected search radius and we evaluate the quadratic polynomial function that provides the best parametrization of the relation between the considered elevation change estimates and their mean elevation.

•
We assign the considered grid cell an elevation change value calculated using the parameters from the quadratic polynomial fit and the average elevation of the grid cell provided by ArcticDEM.
The ArcticDEM is a high-quality DEM obtained using optical stereo imagery from the DigitalGlobe constellation and high-performance software from the National Geospatial-Intelligence Agency (NGA) and the National Science Foundation (NSF) [76], similar to what was produced for Greenland [77] but extended to the entire Arctic.
To define a realistic distance to be used in the extrapolation procedure, we evaluate the semivariogram of the elevation changes measured using the plane fit. We find an autocorrelation distance (lag-correlation) equal to 50 km. We therefore assume that ice elevation changes remain correlated for distances smaller than 50 km and we set the search radius equal to one half of the estimated autocorrelation distance. Further details regarding the spatial extrapolation algorithm are provided in Appendix C.

From Volume to Mass Change
The ice volume change for each grid point is calculated by multiplying the relative mean elevation change value with the grid cell area. The total ice volume change is calculated by summing the contributions from all grid cells. The lack of data regarding possible changes in ice elevation due to snow/firn compaction complicates the conversion of ice volume change to mass volume change [78]. Following Moholdt et al. [11], we calculate the total ice mass change by multiplying the total volume change with the density of ice (0.917 g/cm 3 ) and consider the effect of this assumption in our error budget. We provide our final estimates considering uncertainty terms related to: elevation change measurement and extrapolation error; error in the considered glacier area; error associated with the non-uniform distribution of the elevation change measurements on the glacier surface; error associated with the volume to mass conversion. A detailed description of the calculation of the different error components is provided in Appendix D.

Atmospheric Temperatures and Total Precipitation
We analyze surface temperatures using ERA-Interim reanalysis data [79] available at (http:// www.ecmwf.int/en/research/climate-reanalysis/era-interim). We employ "monthly means of daily means" of 2-m temperature. ERA-Interim has a resolution of 0.75 • that is coarse versus the size of the region of interest. We therefore compare the ERA-Interim linearly interpolated data to monthly temperature measurements obtained from the Hydrometeorological Information-World Data Centre Baseline Climatological Data Sets (http://meteo.ru/english/climate/temp.php) at two meteorological stations available in the region: Malye Karamakuly (59 E, 70.4 N-WMO ID: 20744) and E.K. Fedorova (52.7 E, 72.3 N-WMO Id: 20946) (See Figure 1b). The comparison shows an excellent agreement between the two datasets. In both cases, we find a root mean squared difference below 0.7 • C and a correlation >0.98. We therefore pursue the regional analysis using only the monthly temperatures from ERA-Interim. We evaluate the seasonal mean temperatures and temporal anomalies in reference the time period 1979-2002 in order to investigate climate variability in the region. We also use Synoptic Monthly Means of Total Precipitation to evaluate accumulation variability during the period under analysis. We do not employ data from meteorological station in this case since these data are considered biased and underestimate solid precipitation [11,80].  Table 1  . In red, regional glacier mass change using ICESat between October 2003 and October 2009 (range of variability in mass change is red dash). In green, regional mass change time series calculated using CryoSat-2 between July 2010 and July 2016 with vertical error bars. Table 2 shows the ice mass change trend (and its components) measured by GRACE for the entire mission This comparison shows that after 2010, the NZEM glaciers lost mass at a rate of 14.3 ± 4 Gt/year, hence significantly higher than for the entire period. However, when analyzing the mass change time series over the entire period under observation, the AIC c indicates that a quadratic fit is not a significantly better fit to the regional time series than a linear model. This means that despite the increased glacier mass loss after 2010, the entire time series does not exhibit the characteristics of a statistically significant accelerated mass loss. These results are confirmed by the R 2 adj criterion calculated for the two models. We find the same value 0.8 for both the linear and quadratic models. Between 2003 and 2009, we calculate a mass loss of 9.9 ± 5 Gt/year with a significant positive acceleration equal to +3 ± 1 Gt/year 2 , indicating a slowdown in mass loss after 2008. The regional glacier mass balance becomes positive for two consecutive years. For the years 2009 and 2010, we calculate a total annual mass balance of +4.5 and +9.8 Gt (See Table 1). After 2010, NZEM glaciers lose mass at an increasing rate. Between 2002 and 2017, glacier mass balance was negative for 11 of the 13 considered glaciological years, contributing to a cumulative ice loss of 124 Gt equivalent to 0.32 mm of SLR. Elevation change estimates obtained using ICESat (Figure 3) are consistent with Moholdt et al. [11], with a mass loss of 7.8 ± 5 Gt/year (versus 7.6 ± 1.2 Gt/year). ICESat measurements show that the mass loss is associated with pronounced glacier thinning at elevations below 500 m. During the same time period, GRACE measures a glacier mass loss equal of 9.9 ± 5 Gt/year, which agrees with error bounds.

Results
Elevation change estimates from CryoSat-2 in Figure 4 yield a mass loss of 13.3 ± 5 Gt/year, in agreement with mass change estimates obtained using GRACE, hence increasing confidence in the CryoSat-2 analysis. CryoSat-2 returns valid observations for 67% of the total glacier area. These values are one order of magnitude higher than the 9% obtained with ICESat. As for ICESat, the mass loss is associated with marked ice thinning at low elevation, with maximum thinning at glacier fronts. More stable conditions are observed above 700 m on both the Barents and Kara Sea coasts.  To assess the impact of coast and terminus type in glacier thinning rates, we compare thinning rates measured at low elevations (below 500 m) for the following subregions: (1) Marine-Terminating Glaciers-Barents Sea (8,703 km 2 ); (2) Marine-Terminating Glaciers-Kara Sea (7,016.4 km 2 ); (3) Land-Terminating Glaciers-Barents Sea (2,260.9 km 2 ); (4) Land-Terminating Glaciers-Kara Sea (2,277.9 km 2 ); and (5) Minor Glaciers (1,612.3 km 2 ). Results are presented in Table 3. We employ the Kruskal-Wallis test [81] to evaluate if the mean elevation rate of the five subregions is significantly different. We find that thinning rates on the Barents Sea are significantly higher than thinning rates than on the Kara Sea. In both cases, tidewater glaciers exhibit thinning rates at low elevation significantly higher than thinning rates of land-terminating glaciers on the same coast. This result is consistent with Melkonian et al. [13] and Carr et al. [14]. The high thinning rates at low elevation for tidewater glaciers facing the Barents and the Kara Sea suggest that ice dynamics plays a significant role in the glacier mass balance in the region.  CryoSat-2 provide a more uniform coverage of the drainage basins. Scatter plots for marine-and land-terminating glaciers facing the Kara Sea (Figures 5 and 6) show a pattern of increasing thinning rates with decreasing elevation. Figures 7 and 8 for the Barents Sea are discontinuous and more affected by noise. In almost all glaciers, measurements at elevations between 400 and 600 m are missing. On the Barents Sea coast, glaciers are characterized by a more complex morphology with more rugged surfaces and steeper slopes, together with a smaller average glacier size. Marine-terminating glaciers facing the Kara and Barents Seas have average areas of 559 km 2 and 306 km 2 respectively. Land-terminating glaciers facing the Kara and Barents Seas have average areas of 14.8 km 2 and 6.3 km 2 respectively. However, there are more glaciers facing the Barents Sea, and the total glaciated area on the Barents Sea side is greater than that on the Kara Sea side. These differences limit the ability of CryoSat-2 to measure elevation changes in these regions. However, in all considered cases, a consistent pattern of increasing ice thinning rates with decreasing absolute elevation remains detectable.
We compare elevation changes obtained with CryoSat-2 with those obtained by ICESat over the same drainage basins. Giver the small footprint size (one order of magnitude smaller the CryoSat-2 footprint) and the smaller size of the glacier surface considering during the plane fit, the ICESat measurements tend to be more accurate and less affected by noise, yet we observe the same pattern of dh/dt with elevation. Hence the use of larger plane radius in the CryoSat-2 data does not affect our ability to detect elevation changes at the drainage basin scale. The number of observations with CryoSat-2 increases by one order of magnitude with respect to the ICESat, hence allowing a more complete characterization of glacier changes over the entire region. The highest thinning rates are measured by CryoSat-2 over marine-terminating glaciers, especially for those facing the Barents Sea with peak values of 3-4 m/year below 200 m elevation. Thinning rates over land-terminating glaciers are in general well below the rates measured for tidewater glacier and always below 1 m/year. Total Precipitation was above the mean during summer and autumn. The same pattern of increased precipitation and temperatures was observed by Moholdt et al. [11]. They interpreted this result as evidence of regional climatic mass budget not deviating from values registered in the previous decades. The regional glacier mass balance was positive in 2009 and 2010. These two years were characterized by increased precipitation in almost all seasons and winter and summer temperatures being 0.2 • C below the long-term average in 2010.  [82,83] in order to investigate the presence of a trend in seasonal temperatures during the period 2002-2016. We find a significant positive trend during autumn (2.4 • C ± 1.5 degree/decade) and summer (1.1 • C ± 0.6 deg/decade ).
Zeeberg and Forman [32] observed a strong linear relation between mean summer temperatures and ablation. The analysis of seasonal temperatures is therefore useful to validate the seasonal/interannual variability present in the GRACE/CryoSat-2 mass/elevation change time series. We compare time series of mean summer temperature with annual mass balance from GRACE (see Figure 10). We find that the two time series are significantly anti-correlated (R = −0.69) confirming that summer atmospheric temperatures are an optimal predictor of the annual glacier mass balance variability.

Discussion
We present a 14-year mass balance estimates of the NZEM glaciers extending for the entire GRACE mission. Our GRACE estimates agree with Reager et al. [6] for the time period 2002-2014 and Matsuo and Heki [84] for the time period 2004-2012, both obtained using GRACE data products from the 5th release. We find a large disagreement, however, with Sun et al. [39] who reported a mass loss rate of 1.04 ± 0.25 Gt/year between 2003 and 2014 (versus 7.5 ± 4 Gt/year here). We attribute the discrepancy to differences in data post-processing applied to the GRACE data prior to obtaining regional mass change estimates. Moholdt et al. [11] estimated a mass loss using GRACE ranging between 4.1 and 5.8 Gt/year between October 2003 and October 2009. We find a larger mass loss of 9.9 ± 5 Gt/year. Even though the two estimates agree within their error bounds, estimates by Moholdt et al. [11] were obtained using an older GRACE data release (RL04) and a different GIA correction (ICE-5G, [85]).
Our analysis covers a time period two times longer than the one considered in Moholdt et al. [11] and extends by several years the analysis of Reager et al. [6]. Between 2002 and 2016 the NZEM was the 5th largest contributor to SLR among all ice-covered regions outside the Greenland and Antarctica [7]. Our mass balance assessment combines independent measurements from satellite altimetry and time-variable gravity, which increases the overall confidence in our assessment. We find an interannual variability in the time series of glacier mass change that is similar to that observed on the Svalbard and the Franz Josef Land archipelagoes [7]. In these regions, an increased glacier mass loss between 2010 and 2017 was associated with enhanced ice thinning at the margin of all marine-terminating glaciers [62,86]. This pattern suggests that the NZEM glacier mass balance is largely influenced by shifts in atmospheric and ocean conditions over the entire Barents Sea region. The analysis of temperature and precipitation shows that we may divide the period of study in three main phases. In the first phase between 2002 and 2009, an increase in winter and summer temperatures together with positive anomalies in precipitation yielded a negative mass balance. In the second phase, between 2009 and 2010, below normal summer temperature and increased precipitation yielded a pause in the mass loss. Annual and Winter NAO indexes were predominantly negative during these two years [84,87]. The Winter of 2009/2010 was associated with a record negative NAO index in the entire 187-year record [87]. In the last phase, 2010 to 2017, the record mass loss was associated with atmospheric temperatures above normal in all seasons and NAO returning to strongly positive values. The only three years with positive mass balance (2004, 2010 and 2014) were characterized by colder than normal summers.
We conclude that the time series of ice mass change reflects a long-term process of warming over the Barents Sea region, strongly modulated by alternating phases of positive and negative NAO which either enhance the warming or pause the warming, respectively. During positive NAO, a strengthening in the northward-flowing North Atlantic current increases the advection of warm water into the Barents Sea and an intensification of the polar vortex increases precipitation. Decreasing NAO values sill decrease in the influx of warm AW and decrease atmospheric and ocean temperature. With enhanced positive NAO after 2010, the glacier mass loss of the NZEM has increased markedly.
Importantly, we report greater (2-4 times) thinning of the marine-terminating glaciers compared to land-terminating glaciers, suggesting that ice dynamics plays a strong role in the glacier mass balance of NZEM. Changes in ice dynamics involve increases in ice discharge which increases glaciers' thinning compared to the surrounding areas. Faster glacier flow has been documented in the Arctic and Greenland to be related to ice front conditions. A sufficient condition to de-stabilize the glacier fronts is to increase subaqueous melting in contact with ocean waters [88,89]. Surface thinning may contribute to grounding line retreat. Enhanced melt water production and ocean temperature may increase calving. The melting of ice mélange in the fjords in warm summers or longer periods over the year may flush out icebergs from glacier fronts and facilitate the detachment of more ice blocks. By far, however, it is likely that the dominant process triggering faster flow is an enhanced melt rate of ice in contact with ocean waters [90]. There are few data in this region to confirm the presence of warm, salty, subsurface water of Atlantic origin or to determine whether the intrusion of AW into the fjords has increased over time. Our results call for additional studies of ocean conditions in front of the NZEM glaciers in the Barents Sea, following the work of Enderlin et al. [91] to use iceberg melt rates as proxies for oceanic thermal conditions.
The large glacier elevation changes presented here could be associated with the dynamic activation of surge-type glaciers as observed in other regions facing the Barents Sea, such as Svalbard [15,62,92]. In the case of NZEM, however, Carr et al. [14] documented that only three glaciers were actively surging between 2000 and 2013, while no evidence of surge-like behavior were observed between 2012 and 2014 by Melkonian et al. [13]. Therefore, even though we cannot rule out the presence of surge events during the period of study, we can assume that they have a limited impact on our final estimates.
We calculate that almost 60% of the total ice mass loss takes place over marine-terminating glaciers on the Barents Sea coast. Our results suggest that ice dynamics plays a dominant role in the ongoing mass loss of the NZEM glaciers. This is in contrast with the glaciers in the Canadian Arctic Archipelago and Alaska, which are dominated by land-terminating processes, with only a small fraction controlled by marine-terminating processes.

Conclusions
We present a combined analysis of multi-sensor data, using GRACE time-variable gravity for the entire mission, ICESat altimetry for the entire mission and CryoSat-2 altimetry until present to document the mass loss from NZEM glaciers, RHA. We find that the mass loss has been increasing over the entire period, with two brief pauses, and with a marked increase after 2010. This long-term trend and interannual variability are consistent with a sustained warming of the region part of a long-term trend, modulated by alternative phases of negative and positive NAO, which pause or reinforce the mass loss from NZEM glaciers. With NAO remaining positive after 2010, the mass loss is now increasing faster than in the previous decade. We find excellent agreement between the remote sensing techniques. Our altimetry estimate suggests the present-day mass loss is dominated by marine-terminating glaciers. In addition, the overall ice mass loss from the NZEM glaciers is well above the estimated SMB change. Finally, our results indicate that the ocean plays a major role in the evolution of the NZEM glaciers.
Author Contributions: E.C. and I.V. wrote the paper. E.C. and T.C.S. processed the data. All authors discussed the results and commented the manuscript.

Appendix B. Time Series of Ice Elevation Change
The CryoSat-2 mission has a repeat cycle of 369 days and a sub-cycle of 29 days [63]. The same orbit is repeated with a small offset every month. Such a high repetition frequency allows the system to detect changes in ice elevation at a quasi-monthly temporal resolution. Following Wouters et al. [74] and Noël et al. [75], we evaluate monthly elevation anomalies at each 1 km grid centroid summing the constant trend obtained applying the plane fit to the grid cell residuals. We parameterize the elevation anomalies as a function of absolute elevation. We calculate the mean elevation change time series for each 50 m elevation bin. To obtain time series of regional ice volume change, we multiply single bin mean elevation change time series by the relative area. We obtain the regional glacier mass change time series by multiplying the sum of all the 50 m elevation bin volume change time series with the density of ice. We express the uncertainty associated with the single bin monthly volume anomaly as the 1-σ of the elevation residuals used to calculate the mean elevation anomaly value multiplied for the bin area. Different elevation bins are considered uncorrelated. We therefore sum their uncertainties in quadrature to evaluate the total monthly volume anomaly uncertainty. We include the error introduced by the volume to mass conversion in our ice mass change uncertainty estimate.

Appendix C. Spatial Extrapolation-Local Regression Filter
Before proceeding to the extrapolation phase, we apply a noise reduction procedure. For every grid cell with a valid dh/dt value, we find all the dh/dt estimates located at a distance smaller than or equal to the search radius distance. We use linear regression to find the quadratic polynomial function providing the best parametrization of the relation between the considered elevation change estimates and their absolute elevation. If the dh/dt value of the considered grid point is associated with a residual value larger than 3 times the standard deviation of all the residuals, it is marked as an outlier and excluded from the following extrapolation phase. We further discard anomalous positive elevation change rates with magnitudes above 2 m/year. To validate our selection criterion, we evaluate the maximum accumulation rate in the region using Synoptic Monthly Means of Total Precipitation from ERA-Interim. We cumulate the sum of liquid and solid precipitation between July 2010 and July 2016 and estimate their trend. We find maximum thickening rate of 0.6 m/year, which is well below our selected threshold. We finally reject thinning rates above 10 m per year that are well beyond maximum thinning rates observed over the glacier of NZEM by Melkonian et al. [13] and on the fastest glaciers of the Greenland and Antarctic ice sheets [71].

Appendix D. Uncertainty Analysis
We calculate the total uncertainty characterizing our final glacier mass change considering inaccuracy terms due to :
Sampling Bias Error (due to the non-uniform distribution of the elevation change measurements on the glacier surface) ( bias ); 5.
Error associated with the volume to mass conversion ( ρ ); We evaluate measurement error and extrapolation error following the approach used by Nilsson et al. [72]. We organize the measured elevation changes in bins regularly spaced by a distance equal to the autocorrelation length. Following Nilsson et al. [72], we consider elevation changes located in different bins to be uncorrelated. We calculate the elevation change mean measurement error per unit area ( dh/dt ) as: where N is the number of non-empty uncorrelated bins and σ dh/dt is the mean standard deviation of the elevation changes calculated for each bin. Following Nilsson et al. [72] we calculate the mean extrapolation error per unit area ( ext ) as: where, N is again the number of non-empty uncorrelated bins and σ ext is the mean standard deviation of the extrapolated elevation changes located within each bin. We evaluate the Sampling Bias error ( bias ) by performing a Monte Carlo simulation. We generate a synthetic elevation change map on the 1 km by 1 km grid used for the averaging procedure. We generate 1000 random sampling distributions defined in such a way that a fixed portion of the total glacier area is sampled at each iteration. We evaluate the sampling bias as the standard deviation of the difference between the value obtained in the 1000 iterations and the real mass change value. We estimate an error equal to 0.2 Gt/year and we include it in the error budget. Elevation measurements over glacier surface are isolated employing glacier outlines designed using satellite imagery from the same decade of the two missions [23,24]. Following Moholdt et al. [11], we assume the uncertainty associated with errors in glacier outlines to be within 5% of the total glacier area (A tot ). We calculate the volume error vol as the sum in quadrature of all considered error components multiplied by the total glacier area: We finally include in the total error budget the error associated the volume to mass conversion following Moholdt et al. [11]: where ρ ice is the density of ice,V is the estimated volume change and ρ is calculated as ρ = 1 2 (ρ ice − ρ f irn ).