Spread of Svalbard Glacier Mass Loss to Barents Sea Margins Revealed by CryoSat‐2

The Norwegian Arctic archipelago of Svalbard is located in the most rapidly warming area of the Arctic, at the interface of Arctic and Atlantic air and ocean masses. The presence of a large number of surge‐type glaciers and the potential for rapid changes in surface mass balance and ice dynamics necessitates regularly updated mass balance assessment. This study uses swath processing of CryoSat‐2 SARIn mode data to obtain glacier elevations for 2011–2017. Individual elevation estimates are collected into 1‐km2 grid cells, and a least squares plane‐fitting technique is used to calculate rates of elevation change, with residuals being used to reveal the temporal pattern. A 7‐year rate of mass change of −16.0 ± 3.0 Gt a−1 is estimated (equivalent to 0.044 mm a−1 of global sea level rise), of which −11.0 Gt a−1 results from the melt and dynamic thinning of nonsurging ice and −5.0 Gt a−1 results from surges. This compares to ‐3.4 Gt a−1 previously estimated using ice, cloud, and land elevation satellite (ICESat) (2003–2008). The west coast remains a major contributor to mass loss from nonsurging ice in the archipelago, with mass loss increasing from areas bordering the Barents Sea. Sea ice concentration and climate reanalysis data sets show ocean and lower atmospheric warming and sea ice decline in this region, likely contributing to enhanced glacier melt and discharge.


Introduction
The Arctic is one of the most rapidly warming areas of the planet (AMAP, 2017b), in part due to the polar amplification of climate change (Serreze & Barry, 2011). This has resulted in pan-Arctic glacier mass loss, dominated by increased surface melting (Box et al., 2018;Gardner et al., 2011;Larsen et al., 2015;Noël et al., 2017). The Norwegian high Arctic archipelago of Svalbard lies at the confluence of cold air masses from the Arctic, and warm, humid air masses from the Atlantic (Hagen et al., 1993). The northernmost branch of the North Atlantic Current, the West Spitsbergen Current, transports heat toward the Arctic Ocean along the west coast of Svalbard, preventing extensive sea ice formation to the west of the archipelago even in winter, and providing an ameliorating influence on atmospheric temperatures (Hagen et al., 1993). A further branch of the North Atlantic Current enters the Barents Sea to the south of Svalbard. This is an Arctic warming "hot spot," exhibiting rapid ocean and atmospheric warming and sea ice decline (Årthun et al., 2012;Lind et al., 2018;Polyakov et al., 2017;Screen & Simmonds, 2010b). This climatic setting leads to strong gradients in temperature and precipitation across the archipelago, and a high sensitivity to heat transport from the south and sea ice conditions to the north (Hagen et al., 1993;Isaksson et al., 2005;Lang et al., 2015;Østby et al., 2017;van Pelt et al., 2019).
Approximately 34,000 km 2 of Svalbard is covered by ice (Nuth et al., 2013;Pfeffer et al., 2014), with around two thirds of this area terminating in ∼800 km of calving cliffs (Błaszczyk et al., 2009;Nuth et al., 2013). The ice cover varies from cirques, valley glaciers and icefields on Spitsbergen, to large ice caps on the eastern islands (Hagen et al., 1993). A large number of these glaciers have been observed to surge (Dowdeswell et al., 1991;Hagen et al., 1993;Jiskoot et al., 1998;Lefauconnier & Hagen, 1991;Sund et al., 2009), or exhibit landforms indicative of past surges (Farnsworth et al., 2016). Estimates for the percentage of Svalbard's glaciers which are surge-type vary from 13% to 90% (Jiskoot et al., 1998;Lefauconnier & Hagen, 1991). Observations of a number of recent surges in Svalbard have suggested stronger coupling to climatic warming than previously appreciated Nuth et al., 2019). Indeed, the surges of Storisstraumen (Basin-3) and Nathorstbreen are the largest observed in Svalbard since the comparably warm 1930s McMillan et al., 2014;Nuth et al., 2019;Sund et al., 2014). The high sensitivity of Svalbard glaciers to climatic shifts, and widespread surge potential necessitates regular mass balance assessment. Here, swath processing Gray et al., 2013Gray et al., , 2015Gray et al., , 2017 of roll-angle-corrected (Scagliola et al., 2018;Gourmelen et al., 2018) CryoSat-2 SARIn mode radar altimeter data is used to quantify the elevation change and mass balance of the glaciers and ice caps of Svalbard for the period 2011 to 2017. These estimates are compared to published estimates for the 2003-2008 ice, cloud, and land elevation satellite (ICESat) period  revealing both the overall change in the mass balance of the archipelago and the geographic distribution of changes. Potential oceanic and atmospheric drivers of changes in the magnitude and distribution of mass balance are investigated using sea ice concentration from satellite microwave radiometry and temperatures from atmospheric reanalysis (sea surface, 2 m, and 850-mb geopotential height) during the ICESat and CryoSat-2 periods.

Data
For the purposes of data processing and analysis, Svalbard was divided into six subregions modified from . Spitsbergen, the largest island, was divided into three subregions: northwest, northeast, and south. Nordaustlandet was divided into two subregions dominated by the two large ice caps, Austfonna and Vestfonna. The sixth subregion encompasses the smaller ice caps on Barentsøya and Edgeøya. The 647-km 2 ice cap covering Kvitøya (Kv in Figure 1d) is outside of the CryoSat-2 SARIn mode mask (Wingham et al., 2006) but, for the purposes of extrapolation, was considered part of the Austfonna subregion. Figure 1d shows this subregional division.
The ice extent within these subregions was defined by the Svalbard Glacier Inventory of Nuth et al. (2013), consistent with Randolph Glacier Inventory 6.0 (RGI6.0) (R. G. I., Consortium, 2017;Pfeffer et al., 2014). The glacier outlines were modified to account for the surge advances of Nathorstbreen in south Spitsbergen (Sund et al., 2014) and Storisstraumen (Basin-3) on Austfonna McMillan et al., 2014) by digitizing the advanced fronts using satellite imagery.

CryoSat-2
A satellite radar altimeter emits a pulse that illuminates a subsatellite footprint with a diameter of several kilometers. In an idealized situation, the backscatter returned to the satellite initially increases as the pulse returns from the point of closest approach on the surface (POCA) and continues to increase as an increasing area of the surface surrounding the POCA contributes to the backscatter. The backscatter decreases slowly following its peak as a decreasing ring of the surface is illuminated (Massom & Lubin, 2006). Plotting backscatter against time yields a waveform with a steep leading edge and a gradual trailing edge, though this can be modified by complex topography and subsurface returns. The time delay to the POCA is found by "retracking" the waveform. This involves the application of an algorithm that selects a point on the waveform representative of the delay time to the POCA. Retrackers developed for land ice applications tend to focus on the leading edge of the waveform, employing a backscatter threshold (Helm et al., 2014;Wouters et al., 2015) or locating the steepest slope (Gray et al., 2013;Nilsson et al., 2016). Over a flat surface, the point of closest approach will be at satellite nadir. However, the POCA will be displaced upslope where the subsatellite terrain is sloping. Knowledge of the subsatellite topography is necessary to convert the delay time to a geolocated height estimate. This has been shown to introduce considerable biases when compared to surface global navigation satelllite system (GNSS) transects and airborne and spaceborne laser transects (Bamber et al., 1998;Brenner et al., 2007Brenner et al., , 1983Hurkmans et al., 2012;Schröder et al., 2017).
The European Space Agency's CryoSat-2 is a radar altimetry satellite designed for cryospheric monitoring, especially to overcome the shortcomings of conventional radar altimeters over the sloping terrain characteristic of glaciers, ice caps, and the margins of the ice sheets (Drinkwater et al., 2004;Wingham et al., 2006). It carries a single science instrument, the Synthetic Aperture Interferometric Radar Altimeter (SIRAL), a pulse-limited Ku band (13.6 GHz) radar altimeter with dual receiving antennae placed across track. SIRAL can operate in low resolution mode (LRM), synthetic aperture radar mode (SAR) and interferometric SAR mode (SARIn) (Wingham et al., 2006). In SARIn mode the Doppler shift of backscatter from fore and aft the satellite is used to divide the footprint into ∼300-m along-track strips, and the intensity and phase of the backscatter is measured by the two across-track receiving antennae. The phase shift between the two antennae can then be used to locate the scatterer in the across-track direction, allowing precise geolocation of the POCA. This mode is used over glaciers, ice caps, and the margins of the ice sheets.
By their nature, POCA points cluster at high points in the terrain. For an ice cap this will be the topographic divides between drainage basins, potentially biasing elevation change assessments toward the more stable slow-flowing regions at the expense of the often rapidly changing margins (Gray et al., 2015). In more mountainous terrain this may be mountain peaks and ridges or the higher parts of glaciers and icefields, decreasing the number of on-ice points and biasing measurements to the interior of icefields while undersampling the lower reaches of glaciers . In areas where CryoSat-2 is operated in SARIn mode and the across-track slope lies within a suitable range, it is possible to derive a swath of geolocated elevation estimates from delay times greater than that of the POCA Gray et al., 2013). In the case of an across-track slope, the POCA will be displaced upslope of satellite nadir. Beyond the POCA the illuminated footprint will include a weakly illuminated area upslope and a more strongly illuminated area downslope, beneath the satellite and within the main beam. Where the across-track slope lies within the range ∼0.5°to ∼2.0°the interferometric phase and delay time can be used to derive an across-track swath of geolocated elevation estimates. This condition is met across large portions of ice caps and, to a lesser extent, icefields. Under ideal conditions swath processing can provide more than an order of magnitude increase in the number of geolocated elevation estimates over POCA alone and provide estimates for downslope areas in which there are fewer POCA estimates Gray et al., 2013).

ICESat
ICESat was operational between 2003 and 2009. It carried a single science instrument, the Geosciences Laser Altimeter System (GLAS), along with instrumentation for precise pointing angle determination. The GLAS instrument consisted of three lasers designed to measure land elevation and the vertical structure of clouds and aerosols in the atmosphere (Zwally et al., 2002). The lasers illuminated a ∼65-m diameter footprint every 172 m along track, deriving the mean elevation within the footprint from the two-way travel time between the centroids of the emitted pulse and a Gaussian fitted to the returned waveform (Zwally et al., 2002). The small footprint reduces the effect of surface slopes, a significant improvement over the previous satellite altimeters which used radar. The first laser failed in March 2003 after around a month of operation, resulting in a switch from continuous operation to seasonal campaigns (Schutz et al., 2005). The remaining lasers were operated for around 30 days, two or three times a year until 2009. This prolonged the mission lifetime and allowed observation of both seasonal and interannual ice elevation changes.

GRACE Mascons
Satellite altimetry is an indirect method for estimating glacier mass change as a density assumption is required for the conversion of measured volume change to mass change. The Gravity Recovery and Climate Experiment (GRACE) satellites, operational between 2002 and 2017, directly measured changes in mass from the Earth's gravity field. Other geophysical mass signals can be modeled and removed to isolate gravity changes relating to the transfer of water from storage as land ice to the ocean. Processed and corrected GRACE data are available among other products as mass concentrations ("mascons") from the Goddard Space Flight Center (Luthcke et al., 2013) and Jet Propulsion Laboratory (Watkins et al., 2015;Wiese et al., 2018Wiese et al., , 2016. Mascons have a coarse spatial resolution, but an approximately monthly temporal resolution. They are thus useful for independent comparison to the mass change time series derived from CryoSat-2. GSFC Mascons 6081 to 6093 and 13856, and JPL Mascons 19, 20, 37, and 38 were downloaded and combined to produce a satellite gravimetric mass change time-series for the Svalbard archipelago. The linear rates (−19.1 Gt a −1 ) and temporal patterns of mass change are consistent between the two processing chains ( Figure S7 in the supporting information). The GSFC solution was chosen for comparison with the altimetric results as the JPL solution was found to exhibit a signal likely a result of S 2 tidal aliasing (Chen et al., 2009;A. Gardner, personal communication, January 23, 2020; Text S2 and Figure S7).

Sea Ice Concentration and Climate Reanalysis
Sea ice concentration and climate reanalysis datasets were used to investigate relevant climate conditions during the ICESat and CryoSat-2 periods. Daily sea ice concentration maps at 3.125-km resolution from satellite microwave radiometry were obtained from the University of Bremen (https://seaice.uni-bremen. de/sea-ice-concentration). These utilize AMSR-E (Advanced Microwave Scanning Radiometer-Earth Observing System instrument onboard the National Aeronautics and Space Administration's (NASA's) Aqua satellite) and AMSR2 (onboard JAXA's Global Change Observation-Water satellite) data and the ARTIST sea ice algorithm (Spreen et al., 2008). Sea surface and atmospheric temperatures were obtained from ECMWF ERA5 climate reanalysis at 0.25°latitude and longitude resolution (https://climate.copernicus.eu/climate-reanalysis). These data were used to derive maps of mean June-July-August (jja) sea ice concentration, sea surface temperature, 2-m temperature, and 850-mb geopotential height temperature for 2004(ICESat) and 2011-2017, and the difference between the two periods.

Digital Elevation Models
A digital elevation model (DEM) is required to select the most likely of the three sets of POCA and swath height estimates resulting from phase ambiguity, and in the extrapolation of elevation change estimates to unsurveyed areas. This study utilizes a 20-m resolution DEM produced by the Norwegian Polar Institute, based on photogrammetric mapping of airborne stereo-imagery from late summer (NPI, 2014. Some unprocessed areas, mainly in south Spitsbergen, were filled with a filtered version (Nuth et al., 2013) of the ASTER GDEM v2 (Fujisada et al., 2012) to replace older DEM source data. The mosaicked Svalbard DEM was validated against CryoSat-2 altimetry data, resulting in a mean difference and standard deviation of 0.4 ± 11.6 m over glacier surfaces. The area that covers the Nathorstbreen basin was replaced by the necessary parts of two tiles from the newer ArcticDEM (Porter et al., 2018) to give a postsurge advance topography that is more representative of the CryoSat-2 period. These tiles are composites of photogrammetric DEMs derived from high-resolution optical satellite imagery acquired between 2011 and 2017 (primarily 2013 to 2016). It was not necessary to update the DEM for the surge of Storisstraumen because CryoSat-2 phase ambiguities could be solved throughout the observation period.

Elevation Change and Extrapolation 3.1.1. Retracker and Swath Processor
This study used the maximum slope retracker and swath processor of (Gray et al., 2013(Gray et al., , 2015(Gray et al., , 2017 to derive POCA and swath heights and locations. All available CryoSat-2 L1b BaselineC data covering the Svalbard archipelago between January 2011 and December 2017 were provided by the European Space Agency. The erroneous roll-angle arrays were updated using the star tracker mispointing angle corrections (Scagliola et al., 2018). For each ∼300 m along track, the data files contain time histories of the power, phase and coherence of the returns from the POCA and surface and subsurface beyond the POCA. Each complex waveform is first smoothed prior to the application of the "retracker." The real and imaginary components of the complex waveform are separately smoothed with a low-pass-filter so that when recombined into power and phase, the reduced phase noise will improve geocoding and reduce height noise. The POCA position is assumed to be the maximum slope of the first significant leading edge of the (power) waveform. The corresponding coherence and phase values are extracted, and POCAs with a coherence of <0.7 are discarded. POCA latitude, longitude, and height above the WGS84 ellipsoid are then calculated from the delay time, satellite positioning and attitude information, and phase. Beyond the POCA the waveform contains backscatter from either side of the POCA. Where there is an across-track slope, the POCA is displaced upslope and there is a backscatter contribution from a region below the satellite, within the main beam, and a contribution from a weakly illuminated, range ambiguous region upslope (see Gray et al., 2013, Figure 1). It is assumed that the differential phase is dominated by returns from the region within the main beam. Phase is unwrapped in both directions starting from a central region of high coherence. The smoothed, unwrapped phase, along with information on the satellite positioning, interferometric baseline and delay time is used to geocode additional height estimates in the across-track swath. Height estimation and geocoding is repeated three times using the original phase and that phase ±2π, and the estimates compared to the reference DEM to select the most likely solution. The resulting swath points for each waveform form a straight line across the ice surface in the cross-track direction. These are converted to a regularly spaced set of geocoded height estimates by straightforward averaging in cross-track segments ∼100 m in length.
The POCA and swath points were not corrected for the penetration of the radar signal into the snowpack. This has previously been shown to introduce a bias of ∼0.5 m on Austfonna compared to ground-based GNSS and airborne laser scanner transects ( Journal of Geophysical Research: Earth Surface significant interannual trend with respect to these reference data is observed, meaning the effect of penetration on the derived dh dt is minimal.

Filtering
A robust multistep filter was applied to remove unreliable points caused by incorrect positioning of the satellite's internal tracking loop. Points were discarded if the tracking loop "lost lock" and missed backscatter from the glacial ice. This effect is most pronounced over the broad region of calving cliffs at the southern margin of Austfonna, where the tracking loop often seems effected by strong, off-nadir returns from open water or sea ice, with the altimeter only regaining lock some distance inland. Most of the resultant (erroneous) points are located off ice and would therefore be removed using the glacier outlines in the next stage of filtering, but a small number are located on ice. The delay time to these points is greater than the delay time to nadir sea level ( Figure S9). To remove them, a threshold of the modal delay time to nadir sea level minus 10 range gates (<5m, to account for variations in sea surface elevation, but remaining below the elevation of marine-terminating glacier fronts) was calculated. Points with delay times greater than this threshold were removed. Mountainous Spitsbergen is also challenging terrain for radar altimery. The POCA may be located off-ice on mountain peaks or ridges, and the large elevation variability over short distances as the satellite passes over mountains and fjords may also result in loss of lock.
In the second stage of filtering echoes from ice covered surfaces were isolated using glacier outlines. A third filter was applied to identify and remove points located on ice for which the tracking loop was not correctly positioned, and part of the leading edge of the waveform was missed. Points were removed if the backscatter value of the first range gate of a waveform was greater than two times the standard deviation from the mean of the first backscatter values of all remaining waveforms from that satellite pass. Finally, any swath points originating from the same waveform as a filtered out POCA were removed. A limited number of additional passes, mostly at the far eastern extreme of the Austfonna ice cap, were also removed after manual examination of waveforms. The ground track of these passes was mostly over the ocean, with glacial ice surfaces at large off-nadir angles and the tracking loop positioning and therefore backscatter of the initial samples varying considerably over short distances resulting in a high standard deviation, and preventing the third stage of filtering identifying and removing them.

Plane Fit Algorithm
Remaining POCA and swath points were grouped into 1-km 2 grid cells. This resolution was found to offer a reasonable compromise between retaining sufficient spatial detail in dh dt mapping and maximising data coverage, as well as limiting topographic complexity within each grid cell. An iterative least squares regression method with outlier (±10 m) removal was used to simultaneously fit a surface plane and to estimate a linear rate of elevation change : where x, y, and t are easting, northing, and time offsets between each point and the mean, c 1 and c 2 are the surface slope gradients in the east and north directions, dh dt is the linear rate of elevation change, and v are the elevation residuals in the regression. Estimates derived from fewer than 10 elevation measurements, fewer than four individual satellite passes, or a temporal span of less than 2 years, were rejected (Wouters et al., 2015). The offsets between each measured elevation and the fit, known as residuals, contain temporal elevation variability as well as errors in the surface fit and elevation measurement. These were used later in the processing chain to add seasonality and interannual variations to the linear rates of elevation change. Figure 1 shows the areal coverage and hypsometric distribution of POCA and swath points, as well as successful dh dt measurements. The POCA point density is generally low, with higher densities confined to prominent ice divides in northeast Spitsbergen, Austfonna, and Vestfonna. The POCA clustering at high points leads to a bias toward higher elevations (Figure 1a, inset), areas likely to exhibit more gradual changes than lower elevations toward the glacier fronts. Swath points (Figure 1b) are always located downslope of POCA points, making them highly complementary datasets, increasing the areal coverage and reducing the bias to higher elevations. The elevation distribution of swath points closely matches the glacier hypsometry ( Figure 1b, inset). Combining POCA and swath points, dh dt estimates were successfully obtained for 58% of the ice area, with a reduced bias toward higher elevations compared to POCA alone ( Figure 1c). Extensive coverage is achieved over areas of relatively simple topography such as the ice caps on the eastern islands (69-84%; Table S1). Coverage in the more mountainous Spitsbergen is less extensive (39-49%; Table S1), and still biased to the interior of the icefields, and hence toward higher elevations. Swath points do not provide as great an increase in coverage in these subregions (Table S1).

Extrapolation
The extrapolation of elevation change estimates to unsurveyed areas requires separation of actively surging ice (where elevation changes can be rapid and dynamically driven) and nonsurging ice (non-surge-type and quiescent surge-type glaciers where elevation changes are more gradual). Surging basins are outlined in black in Figure 1d and were defined by searching for high magnitude elevation changes in the initial CryoSat-2 results and with reference to published observations of surges occurring during the time period of this study (Benn et al., 2019;Dunse et al., 2015;McMillan et al., 2014;Sevestre et al., 2018;Sobota et al., 2016;Sund et al., 2014).
The proportion of surveyed to unsurveyed grid cells vary spatially depending on local topography and across-track slopes. Estimation of regional volume or mass change requires extrapolation to unsurveyed areas, for which a number of methods exist . This study uses a hypsometric polynomial method, which assumes a relationship between elevation and elevation change. Variations on this method are commonly used to extrapolate sparse radar (Foresta et al., 2016Wouters et al., 2015) and laser (Moholdt et al., , 2012Nuth et al., 2010) altimetry measurements to unsurveyed areas. While coverage on the eastern islands is extensive (Figures 2b, 2d, and 2f), coverage on Spitsbergen is more sparse, with greatest coverage in the gradually changing interior (Figures 2a, 2c, and 2e). A parameterization of the

10.1029/2019JF005357
Journal of Geophysical Research: Earth Surface relationship between elevation and elevation change is therefore necessary to avoid underestimation of subregional mass loss.
The glacier inventory was used to extract glacier elevations from the DEM, allowing the calculation of the average ice elevation within each 1 km 2 grid cell. The fractional ice coverage was also calculated. For each subregion separate third order polynomials for surging and nonsurging ice were fitted to characterize the relationship between the elevation of each grid cell and its rate of elevation change. Separate polynomials were fitted for Nathorstbreen and other surging glaciers in south Spitsbergen as they are at different stages in the surge cycle. To avoid unrealistic behavior (large rates of elevation change) in the polynomial fits at data-sparse highest elevations, the small number of elevation change rate estimates and unsurveyed grid cells requiring extrapolation in these intervals were treated as though they were at the centre of the highest 50-m interval containing >1% of subregional ice area. The rate of elevation change in each unsurveyed grid cell was then calculated according to the appropriate polynomial and the average ice elevation within that grid cell without binning the data. Changing the order of the polynomial and the merging of data-sparse high-elevation intervals was found to not significantly influence the result. For example, the estimated archipelago-wide mass loss rate of −16.0 Gt a −1 from third order polynomial compared to −16.4 Gt a −1 from first-order polynomial, with much of the difference due to the Nathorstbreen basin.
CryoSat-2 has a 30-day orbital subcycle, suggesting the possibility to resolve elevation changes at monthly resolution (Ciracì et al., 2018;Noël et al., 2018). A similar hypsometric polynomial approach was used to extrapolate monthly elevation residuals, except that a first-order polynomial was chosen due to the sparser data coverage. Separate surging and nonsurging polynomials were derived describing the relationship between monthly mean elevation residuals and elevation. Monthly maps of smoothed elevation residuals were then calculated from the polynomials and the DEM.

Volume to Mass Conversion
Volume change can be estimated as the elevation change multiplied by the ice area within each grid cell (0-1 km 2 ) as defined by the glacier inventory. Cumulative volume change on a monthly basis was calculated by adding the monthly residual volume to the linear volume change (the linear rate multiplied by the time since January 2011). Conversion of volume change to mass change requires an assumption about density. Huss (2013) presented a comprehensive assessment of density assumptions, recommending a conversion factor of 850 ± 60 kg m −3 for observation periods longer than 5 years, where the mass balance gradient is stable, firn is present, and volume changes are large.
The surge advance of Storisstraumen (Basin-3), Austfonna was separately accounted for by defining presurge and postsurge calving fronts on satellite imagery, and calculating the added volume and mass of the advance area, based on 2017 CryoSat-2 glacier elevations and past bathymetry observations from Solheim (1991).

Error Budget
The error in mass balance derived from CryoSat-2 altimetry is a combination of the error in the dh dt estimates from the plane-fitting algorithm, the extrapolation error, and the volume to mass conversion error. In a linear regression such as Equation 1, the standard deviation of the third parameter dh dt is usually estimated by where D is the design matrix containing the east, north, and time offsets (x,y,t), cov is the covariance matrix, cov (3,3) is the variance of the time parameter, v is the residuals, and N is the number of observations. Since observations along individual CryoSat-2 tracks will be strongly correlated, N is set to the number of separate tracks such that N−3 represents the statistical degrees of freedom in the sense of a standard error. The uncertainty in surveyed volume change was then calculated as the sum of these individual dh dt errors multiplied by the grid cell area (A):

Journal of Geophysical Research: Earth Surface
The extrapolation error (ε _ V extrapolated ) was estimated by randomly removing 50 ± 3% of the dh dt estimates for a given subregion (Text S1) and calculating new hypsometric polynomials and volume change estimates. This procedure was repeated 10,000 times for each subregion, the highest and lowest 5% of estimates discarded (to remove outliers), and the spread of the remaining estimates assumed to be representative of the error in extrapolation. The total volume change uncertainty is then This was then combined with the density conversion uncertainty (ε ρ ¼ ± 60 kg m −3 Huss, 2013) to estimate the mass change uncertainty of a subregion: where dV dt is the volume change rate and ρ is the density conversion factor ( ρ ¼ 850 kg m −3 ). The total mass change uncertainty for the Svalbard archipelago was estimated by combining the mass change uncertainties of the six subregions as root-sum-squares.
For the ICESat results, the volume change errors of  were used as ε _ V sur veyed , and a similar extrapolation error as for CryoSat-2 was calculated. The total errors were then calculated according to Equations 5 and 6.

Comparison of CryoSat-2 and ICESat
Comparing dh dt and mass balance estimates from CryoSat-2 for the period 2011-2017 with estimates from ICESat for the period 2003-2008 reveals the overall change in the mass balance of the archipelago, and allows the mapping of the geographic distribution of changes between the two periods. The geographic distribution of glacier changes can, in turn, be compared to changes in environmental variables to aid in the identification of potential atmospheric and oceanic driving mechanisms.  estimated changes in elevation and volume from ICESat based on a plane fit algorithm and hypsometric polynomial extrapolation applied to GLA06 release 28 data from 2003 to 2008. The archipelago was divided into subregions similar to the ones used here, except that Nordenskiöldland is included in the south Spitsbergen subregion here. To map changes in dh dt each individual ICESat estimate was compared with the CryoSat-2 rate within the 1-km 2 grid cell with which it overlaps. A similar density conversion scheme (850 ± 60 kg m −3 ; Huss, 2013) was used to convert ICESat-derived rates of volume change to mass change.

Elevation and Mass Changes
The topography of the archipelago is the major factor controlling the spatial extent of successful dh dt retrieval.
On mountainous Spitsbergen (Figure 1c and histograms in Figures 2a, 2c, and 2e) the coverage varies between 39% and 49%, with a moderate bias toward higher elevations. The coverage is more extensive (69-84%) over the relatively simple topography of the ice caps of the eastern islands (Figure 1c and histograms in Figures 2b, 2d, and 2f), where swath processing provides a significant increase in the number of elevation estimates (Figure 1b). The apparent sampling bias to higher elevations of Austfonna is a result of the low-altitude ice cap Kvitøya being outside the SARIn mode mask.

Journal of Geophysical Research: Earth Surface
Across the archipelago between 2011 and 2017 there is a general pattern of thinning at low elevations and smaller changes at high elevations (Figure 2), superimposed on a southwest to northeast gradient in thinning ( Figure 3). These general patterns are punctuated by a number of surges. Thinning is prevalent across the two southernmost subregions, south Spitsbergen and Barentsøya-Edgeøya (Figures 2e and 2f), with the highest (nonsurging) rate of thinning (∼3 m a −1 ) occurring at the lowest elevations in south Spitsbergen. Though rates of thinning are lower at higher elevations, they are still considerable (∼1 m a −1 in south Spitsbergen, ∼0.75 m a −1 in Barentsøya-Edgeøya) around the peak in hypsometry at 400 m. Thinning persists to higher elevations in south Spitsbergen, while the highest elevation ice in Barentsøya-Edgeøya is close to balance, though represents only a small percentage of the glacier area in the subregion. In northwest Spitsbergen the rate of thinning transitions from ∼2 m a −1 in the lowest 250 m to ∼0.75 m a −1 at intermediate elevations, with ice above ∼700 m close to balance (Figure 2c).
Rates of marginal thinning in the subregions of the north and east of the archipelago (northeast Spitsbergen, Austfonna, and Vestfonna) were similar to those in the other subregions. The lowest elevation ice in northeast Spitsbergen and Austfonna (Figures 2a and 2d) thinned at rates of around 2 m a −1 , while values were

10.1029/2019JF005357
Journal of Geophysical Research: Earth Surface and 2e) was most pronounced between 100 and 600 m (an elevation interval containing ∼63.5% of the island's ice). The greatest increase occurred in the southernmost subregions, with thinning increasing by ∼1 m a −1 in south Spitsbergen and Barentsøya-Edgeøya (Figures 2e and 2f). Austfonna maintains an area of high-elevation thickening, but the rate has decreased and moved ∼100 m up glacier (to elevations above 400 m) since the ICESat period (Figure 2d). The ICESat polynomial for Vestfonna (Figure 2b) was influenced by low-elevation thickening of the surging glacier Franklinbreen (Pohjola et al., 2011;Schäfer et al., 2014), but is otherwise unchanged.
The distribution of mass loss from the archipelago reflects the pattern of elevation change and the area of ice within each subregion (Figure 4 and Table 1). The greatest sources of mass loss are northwest and southern Spitsbergen (−3.34 ± 1.28 and −4.57 ± 1.49 Gt a −1 ) and the surge of Storisstraumen on Austfonna (approximately −3.56 Gt a −1 ). Despite the similar pattern of elevation change to southern Spitsbergen, mass loss from Barentsøya-Edgeøya was lower (−1.93 ± 0.60 Gt a −1 ), due to the smaller ice area. The rate of mass loss from the northern and eastern subregions was low (excluding the surge of Storisstraumen). Combined, the total rate of mass loss for the archipelago between 2011 and 2017 was −16.0 ± 3.0 Gt a −1 (Figure 5).
The derived time series of glacier mass reveals seasonal and interannual variability, for example, the strong melt season of 2013 (Lang et al., 2015) and the onset of the surge of Storisstraumen between 2012 and 2013 . Comparison with the satellite gravimetric time series of the Svalbard mass balance ( Figure 5) shows a similar pattern of more moderate mass loss through 2011 and 2012, more rapid mass loss since 2014, and extreme mass loss during the summer of 2013. However, the seasonal mass cycle is not directly comparable due to seasonal variations in CryoSat-2 signal penetration and glacier density.

Climate and Sea Ice Conditions
Summer (jja) sea ice concentration from microwave radiometry and temperatures from ERA5 climate reanalysis over the Barents Sea region were used to investigate potential climatic links with the increase in mass loss from Svalbard glaciers between the ICESat (2003)(2004)(2005)(2006)(2007)(2008) and CryoSat-2 (2011CryoSat-2 ( -2017 periods. A decline in sea ice concentration between the two periods is prevalent across the region (Figures 6a-6c). Sea ice decline is evident around southern Spitsbergen, Barentsøya-Edgeøya, Hinlopen Strait, and to the south of Austfonna. There was little summertime sea ice along the west coast of Spitsbergen during either period. Changes in the pack ice to the north were more modest.
According to ERA5 reanalysis data, sea surface temperatures increased across much of the Barents Sea region (Figures 6d-6f), with warming of approximately 1°C along the western coast of Spitsbergen and into Whaler's Bay to the north of Spitsbergen, and warming of 1.5-2°C in Storfjorden between Spitsbergen and Barentsøya-Edgeøya. This is consistent with observations of an increase in subsurface ocean temperatures in the northern Barents Sea since the mid-2000s (Lind et al., 2018;see Figures 1a and 1b). The spatial pattern of 2-m air temperature change is almost identical to the pattern of sea surface temperature change (Figures 6g-6i), though with smaller magnitudes (+0.75°C along western Spitsbergen, +1°C in Storfjorden). In contrast, temperature changes higher in the atmosphere (850-mb level) were less pronounced, with the greatest increases to the north, and lower magnitude changes in the Barents Sea (Figures 6j-6l). This pattern was interrupted in 2013 ( Figure S10) when westerly atmospheric circulation caused an exceptionally warm summer (Lang et al., 2015) with high mass losses throughout Svalbard (Figures 4 and 5) and the Russian Arctic (Ciracì et al., 2018).

Elevation and Mass Change
In general, thinning is greatest at the low-elevation margins, with lower magnitude changes at higher elevation, including the gradual thickening of some (high elevation) ice cap interiors (Figure 3), a pattern that is common throughout the Arctic (Ciracì et al., 2018;Foresta et al., 2016;Gardner et al., 2011). This pattern is superimposed on a southwest to northeast gradient in thinning. Thinning in the south and northwest is rapid (∼2 m a −1 ) and extends a considerable distance inland from the icefield margins, leaving relatively small areas in balance at high elevation (Figures 3 and 2c, 2e, and 2f). Larger areas of the interior of the icefields in northeast Spitsbergen and the ice caps in Nordaustland are in balance or gradually thickening (<0.5 m a −1 ; Figures 3 and 2a, 2b, and d). Thinning of the northeasterly subregions is generally not as rapid and  Figure 1 of  shows strong thinning across Wedel Jarlsberg Land in the west of south Spitsbergen, partly compensated for by areas of thickening in the east of south Spitsbergen. This compares with widespread thinning in south Spitsbergen shown in Figure 3 (and more negative dh dt from the ICESat to the CryoSat-2 periods in Figure S4). This suggests that southeastern Spitsbergen has switched from near-balance conditions to thinning.
The increased mass loss from south Spitsbergen and continued mass loss from northwest Spitsbergen is primarily driven by surface mass balance and discharge of nonsurging ice. Mass losses from surging glaciers on Spitsbergen have been relatively small during the CryoSat-2 period   (Table 1). The gray area shows the uncertainty in the archipelago-wide mass change rate estimate. The cumulative mass change between 2011 and 2016.5 from satellite gravimetry (dashed black) is derived from GSFC mascons (Luthcke et al., 2013).

10.1029/2019JF005357
Journal of Geophysical Research: Earth Surface (2011)(2012)(2013)(2014)(2015)(2016)(2017). This is in part due to the redistributive nature of surges, and their small area compared to nonsurging ice. The apparent small mass loss resulting from surges in southern Spitsbergen, including that of Nathorstbreen-one of the largest surges ever observed on Svalbard (Sund et al., 2014;Sund & Eiken, 2010)-is likely a result of the most active phase of the surge predating the CryoSat-2 period (Nuth et al., 2019;Sund et al., 2014), and potentially an underestimation resulting from the poor data coverage  (Figure 1), Estimated rates of mass loss for Nathorstbreen calculated by applying a conversion factor of 850 ± 60 kg m −3 to the volume change estimates of Nuth et al. (2019) are −0.91 Gt a −1 between 2008 and 2010, and −0.65 Gt a −1 between 2010 and 2014. The latter estimate partially overlaps with the CryoSat-2 period, over which a rate of mass change of −0.54 Gt a −1 is estimated. For Strongbreen, which surged around 2016 (Benn et al., 2019), an insignificant mass gain of 0.16 Gt a −1 is obtained, indicating mass redistribution rather than major losses. The uncertainty in surge-related mass loss is reflected in the high uncertainty in the rate of subregional mass change (Table 1).
Considering the small ice area, mass loss from Barentsøya-Edgeøya is also considerable at −1.93 ± 0.60 Gt a −1 , with approximately equal contributions resulting from the surge, or frontal destabilisation , of Stonebreen and mass loss from other glaciers and ice caps. This compares to a rate of −1.18 ± 0.15 Gt a −1 between 1971 and 2005 , and −0.39 ± 0.30 Gt a −1 between 2003 and 2008 . The specific mass balance of nonsurging ice (−0.62 m.w.e. a −1 ) in the Barentsøya-Edgeøya subregion is comparable to the rates in the south and northwest Spitsbergen subregions ( Table 1).
Rates of mass loss and the specific mass balance of nonsurging ice were more modest from northeast Spitsbergen (−1.08 ± 1.41 Gt a −1 , −0.13 m.w.e. a −1 ), Vestfonna (−0.58 ± 0.48 Gt a −1 , −0.23 m.w.e. a −1 ), and nonsurging ice on Austfonna (−0.94 Gt a −1 , −0.13 m.w.e. a −1 ), possibly reflecting the generally cooler climatic conditions to the north and east (Lang et al., 2015;Østby et al., 2017;van Pelt et al., 2019). The lower specific mass balance on Vestfonna could be related to the lower elevation of the ice cap compared to Austfonna or the icefields of northeast Spitsbergen (Figure 2). The greater rates of thinning during the ICESat period in the lowermost 100 m of Austfonna and northwest Spitsbergen subregions (Figures 2c  and 2d) may be due to ICESat's smaller footprint better resolving elevation changes near calving fronts.
Prior to 2013, ice in northeast Spitsbergen was on average in balance or thickening slightly (Figures 4a and  S3), a continuation of the near-zero mass balance (+0.44 ± 0.59 Gt a −1 ) observed between 2003 and 2008 by . The year 2013 stands out as a year of strong mass loss, with subsequent years characterised by close to balance conditions or mass losses lower than the 1965-2005 average of −2.08 ± 0.21 Gt a −1 observed by Nuth et al. (2010). Thinning is prevalent across the relatively low-elevation Olav V Land (Figures 3 and S1), particularly at the margins of large marine-terminating glaciers (Hinlopenbreen, Hochstetterbreen, Sonklarbreen, & Negribreen), while the high-elevation portions to the west are in balance or thickening.
Across Hinlopen strait to the east, the long-term thickening trend has continued in the interior of Austfonna . Despite this, the mass change of nonsurging ice is negative, suggesting that surface mass balance has decreased or that dynamic thinning of marine-terminating ice in the south and east of the ice cap (McMillan et al., 2014) has reversed the positive mass balance observed between 2003 and 2008 (+0.52 ± 0.37 Gt a −1 , . The surge of Storisstraumen is the largest single source of mass loss (−3.56 Gt a −1 ) from the archipelago and is comparable to the mass loss from the south or northwest Spitsbergen subregions (Figures 4c-4e and Table 1). Nearby Vestfonna has switched from 1.0 ± 0.29 Gt a −1 of mass gain between 1990 and 2005 , to moderate mass loss in 2003-2008 (−0.33 ± 0.24 Gt a −1 ,  and in 2011-2017 (−0.58 ± 0.48 Gt a −1 ). The increased mass loss is dominated by the thinning of Franklinbreen in the northwest of the ice cap (Figures 3 and S4), and the strong melt season in 2013 (Figure 4b).
In total, the geodetic mass balance of Svalbard between 2011 and 2017 was −16.0 ± 3.0 Gt a −1 . Surging glaciers contributed around one third of the mass loss, despite representing only about one eighth of the ice area. The remaining two thirds of the total mass change resulted from surface mass balance and discharge from nonsurging ice. Contributions from surface mass balance and discharge cannot be easily separated, but the widespread pattern of low-elevation thinning is likely the result of increased surface melt, while the rapid thinning of some large marine-terminating glaciers and fast flowing portions of Austfonna suggests a dynamic component.  Figure S7) compared to CryoSat-2 (−16.0 Gt a −1 ; Figure 5) between 2011 and 2017 could be a result of biases in GRACE corrections, the CryoSat-2 volume-to-mass conversion, or the impact of marine-terminating glacier retreat, which is poorly constrained on Svalbard (Błaszczyk et al., 2009). The results are also consistent with Wouters et al. (2019) who derived global glacier mass change time series for the entire GRACE mission (2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016), with estimated mass change rates of −7.2 ± 1.4 and −9.1 ± 4.1 Gt a −1 based on two solutions from spherical harmonics and JPL mascons, respectively. These estimates fall between those of ICESat (−3.4 ± 1.6 Gt a −1 ) and CryoSat-2 (−16.0 ± 3.0 Gt a −1 ) for the partly overlapping periods.
The CryoSat-2 time series (Figures 4 and 5) exhibit pronounced peaks in mass (equivalent to peaks in elevation) during early summer. These are a consequence of the transition from volume scattering within a predominantly dry snowpack, to surface scattering from a melting snowpack (Gray et al., 2015). This is a seasonal analogue for the raising of the radar reflection horizon accompanying the increase in surface melt extent in Greenland Simonsen & Sandberg Sørensen, 2017;Thomas et al., 2008). In between peaks, the radar does not necessarily "see" the build-up of the snowpack over winter, as the last summer surface remains a strong reflector, and volume scattering results in a retracked elevation below that of the snowpack surface. At the onset of surface melt the water content of the upper snowpack increases, masking volume scattering returns. The radar now ″sees″ the snowpack surface, which manifests as a jump in elevation. Through the summer the radar tracks the melting snow and ice surface. The winter season amplitudes thus do not show the correct temporal pattern of seasonal mass gain. Also, given that a value of 850 ± 60 kg m −3 is used for volume-to-mass conversion of the time series, seasonal mass gain through snowfall and mass loss through melt are overestimated. Nevertheless, the time series are useful to examine the interannual variability in mass change (e.g., the strong melt season in 2013 and the onset of the surge of Storisstraumen between 2012 and 2013) and allow the extraction of estimates of annual mass balance. They also indicate a potential to study end-of-winter snowpack depth and melt season surface lowering on a regional scale, which may allow more accurate estimation of seasonal mass changes using more refined density assumptions or field observations.

Links to Climate and Sea Ice Conditions
Comparing altimetry and climate data for the period 2011-2017 (CryoSat-2) with the period 2003 gives an indication of the origin of this increased mass loss. There is spatial correspondence between regions of increased glacier mass loss (Table 1 and Figure S4), areas of greatest sea ice concentration decrease, and areas of sea surface and lower-atmospheric warming ( Figure 6). Northwest and south Spitsbergen are the major contributors to the mass loss from nonsurging ice (Figures 4c and 4e), but while northwest Spitsbergen and Wedel Jarlsberg Land in western south Spitsbergen display only modest increases in mass loss, mass loss from the eastern coast of south Spitsbergen has increased considerably ( Figure S4). This is consistent with the persistent absence of sea ice (Figures 6a and 6b) and high sea surface and lower atmospheric temperatures (Figures 6d, 6e, 6g, and 6h) on the west coast due to the influence of the West Spitsbergen Current, and the sea ice concentration decline, and strong sea surface and lower atmospheric warming in Storfjorden. Sea ice decline and warming also occurred around Barentsøya-Edgeøya, to the south of Austfonna, and in Hinlopen Strait. Additional increases in mass loss (or switches from mass gain to mass loss) have occurred adjacent to these areas on Barentsøya and Edgeøya, along the calving cliffs of the southern and eastern coasts of Austfonna and in Olav V Land in northeast Spitsbergen. In some of these areas there is evidence of rapid thinning at marine-terminating margins and an increase in inland thinning between the ICESat and CryoSat-2 periods. Regions of low magnitude mass change are limited to Vestfonna, northern Austfonna, and the high-elevation portions of the icefields of northeast Spitsbergen. Here, sea ice concentration did not change substantially, and sea surface and 2-m temperature increases were lower than elsewhere around the archipelago.
Sea ice concentration decline, sea surface and lower atmospheric warming is prevalent across the northern Barents Sea and Kara Sea (Figure 6), and glacier thinning and mass loss has been substantial from the archipelagos of the Russian high Arctic. Glacier thinning in Franz Josef Land has increased in recent years, with greatest thinning detected in the southwest (Zheng et al., 2018). Thinning in Novaya Zemlya has been most rapid on marine-terminating glaciers of the Barents Sea coast (Melkonian et al., 2016), while calving front retreat is prevalent on both the Barents Sea and Kara Sea coasts (Carr et al., 2014(Carr et al., , 2017. Furthermore, the temporal patterns of mass loss appear similar across the archipelagos of the Norwegian and Russian Arctic ( Figure 5 here versus Figure 2 of Ciracì et al., 2018).
Taken together, the results presented here demonstrate a spread of Svalbard mass loss from the west coast glaciers influenced by the West Spitsbergen Current to areas bordering the northwest Barents Sea and Storfjorden, where sea ice concentration has declined and sea surface and lower-atmospheric temperatures have increased. This, combined with observations of increased mass loss from Franz Josef Land and Novaya Zemlya (Melkonian et al., 2016;Zheng et al., 2018), suggests a central role for the warming of the Barents Sea region in recent ice loss from the Eurasian Arctic. The Barents and Kara Seas are experiencing the most rapid regional warming in the Arctic (AMAP, 2017a; Screen & Simmonds, 2010a), with the Barents Sea undergoing a process of "Atlantification" (Årthun et al., 2012;Lind et al., 2018;Polyakov et al., 2017). There has been an increase in Atlantic heat transport into the Barents Sea, due to a strengthening and warming of Atlantic water inflow (Årthun et al., 2012). A reduction in sea ice flux into the ocean basin since the mid-2000s has reduced the freshwater content of the upper water column, weakening stratification and enhancing upward heat flux from underlying Atlantic waters. This reduces local sea ice formation and permits heat loss from the ocean to the lower atmosphere (Lind et al., 2018;Polyakov et al., 2017).
Increased ocean temperatures, decreased sea ice concentration, and increased meltwater production driven by atmospheric warming have the potential to increase mass loss from marine-terminating glaciers and contribute to geometric and hydrological changes, which trigger or sustain surges Sevestre et al., 2018). The impact on surface mass balance is complicated by the competing effects of increased melt due to atmospheric warming, and potential increases in precipitation due to the warmer atmosphere and sea surface, and reduced sea ice cover. However, observed and modeled precipitation trends on Svalbard are weak (van Pelt et al., 2019) and potential increases might be counter-balanced by an increased fraction of rain and reduced pore space in the firn to retain liquid water. Elevation changes alone do not allow the partitioning of increased discharge and reduced surface mass balance with certainty, but there is evidence that both contribute to the observed mass loss. Firstly, there is a widespread pattern of low-elevation thinning across both marine-and land-terminating glaciers which indicates an increase in surface melt. A number of larger marine-terminating glaciers exhibit strong frontal thinning which, combined with observations of increased flow speed (Figure 4 of  is evidence of dynamic mass loss. On Austfonna, thinning is associated with fast flowing ice grounded below sea level (McMillan et al., 2014). There is also increasing evidence that the surges are linked to climatic warming either through increased meltwater penetration to the glacier bed  or the effect of thinning and retreat on conditions near glacier termini (Nuth et al., 2019;Sevestre et al., 2018).

Summary and Conclusions
Swath processing of CryoSat-2 allows the estimation of elevations beyond the point of closest approach, leading to an order of magnitude increase in the number of elevation estimates in ideal conditions. Over the relatively simple topography of the ice caps of the eastern islands of Svalbard, the number of swath points obtained is high, leading to a significant increase in areal coverage and a reduction in the bias toward higher elevations. Over the more complex terrain characteristic of Spitsbergen, POCA points cluster in the high-elevation interior of ice fields, and cross-track slopes are often too steep for the successful retrieval of swath points. This, combined with occasional "loss of lock"-periods when the satellite fails to keep the ice surface within the receiving window over undulating terrain-leads to less extensive coverage of POCA points and a smaller increase in coverage from swath points. Still, by combining POCA and swath data, it is feasible to derive long-term rates of elevation change for >40% of the glacier area, allowing robust estimation of glacier mass balance at a subregional scale on Svalbard.
Derived rates of elevation change show an archipelago-wide pattern of thinning near glacier fronts and more balanced conditions in the interior of icefields and ice caps, interspersed by a number of glacier surges with more rapid elevation changes. The overall glacier mass change of Svalbard is estimated to be −16.0 ± 3.0 Gt a −1 for the period 2011-2017, driven principally by the surge of Storisstraumen and widespread glacier thinning in southern and northwestern Spitsbergen. This corresponds to a global sea level rise contribution of 0.044 mm a −1 . Compared to earlier altimetry data from the ICESat mission (2003)(2004)(2005)(2006)(2007)(2008), mass loss has 10.1029/2019JF005357 Journal of Geophysical Research: Earth Surface increased along the southeastern coast of Spitsbergen to Olav V Land, over Barentsøya-Edgeøya, and the southern and eastern parts of Austfonna, indicating a spread of mass loss from the west coast of Spitsbergen into the rapidly warming Barents Sea region. Sea ice concentration and climate reanalysis data show that Storfjorden and the Barents Sea off eastern Svalbard experienced a summertime decline in sea ice concentration, and warming of the sea surface and lower-atmosphere between the ICESat and CryoSat-2 periods. This warming likely increased surface melt, leading to widespread low-elevation thinning, and may have also increased submarine melt leading to dynamic mass loss from larger marine-terminating glaciers. Increased meltwater availability and the thinning and retreat of glacier termini may also have contributed to the large number of surges observed over the CryoSat-2 period. The continued mass loss from the western coast of Spitsbergen and the spread of mass loss to Barents Sea margins means regions of low magnitude mass change are now limited to Vestfonna, northern Austfonna, and the high-elevation portions of northeast Spitsbergen.