Increased Ice Thinning over Svalbard Measured by ICESat/ICESat-2 Laser Altimetry

: A decade-long pronounced increase in temperatures in the Arctic, especially in the Barents Sea region, resulted in a global warming hotspot over Svalbard. Associated changes in the cryosphere are the consequence and lead to a demand for monitoring of the glacier changes. This study uses spaceborne laser altimetry data from the ICESat and ICESat-2 missions to obtain ice elevation and mass change rates between 2003–2008 and 2019. Elevation changes are derived at orbit crossover locations throughout the study area, and regional volume and mass changes are estimated using a hypsometric approach. A Svalbard-wide annual elevation change rate of − 0.30 ± 0.15 m yr − 1 was found, which corresponds to a mass loss rate of − 12.40 ± 4.28 Gt yr − 1 . Compared to the ICESat period (2003–2009), thinning has increased over most regions, including the highest negative rates along the west coast and areas bordering the Barents Sea. The overall negative regime is expected to be linked to Arctic warming in the last decades and associated changes in glacier climatic mass balance. Further, observed increased thinning rates and pronounced changes at the eastern side of Svalbard since the ICESat period are found to correlate with atmospheric and oceanic warming in the respective regions.


Introduction
Observing variances in global cryosphere characteristics remains an important factor in a warming climate. Consequently, several satellite-mounted laser and radar sensors have been used to continuously measure surface elevation change of glacier and ice caps in the last decades (e.g., [1][2][3][4][5]). While radar systems typically exhibit a large footprint of several meters to kilometers and can penetrate into ice, snow, and firn (e.g., [6]), the advantages of laser altimetry lie in a small footprint of cm to meters and in the measurement of precise surface heights with little subsurface penetration. Glacier topography investigations with laser altimeters were initially mainly conducted through aircraft campaigns (e.g., [7,8]) until NASA launched its first spaceborne satellite laser altimetry mission ICESat (Ice, Cloud, and Land Elevation Satellite), which operated between 2003-2009. The mission data has been successfully used for glacier mass change estimates in non-polar (e.g., [9]) and polar (e.g., [10]) regions. In 2018, its successor, ICESat-2, was launched equipped with an improved laser altimetry system, which is expected to produce height measurements with enhanced accuracy and coverage (e.g., [11]). The first data analyses have shown its potential for cryosphere applications, such as mapping of grounding lines in Antarctica [12] or elevation and mass change assessments for the Greenland and Antarctic ice sheets [11]. However, for smaller glaciated regions like Svalbard, an evaluation of the potential of combining data from both missions through a crossover approach has not been published yet.
In terms of Arctic amplification, the Barents Sea area represents an inter-Arctic warming hotspot [13,14]. Because it is located at the western border of this fast-warming region and at an area of retreating Arctic sea ice, the ice caps and glaciers of Svalbard are exposed to the fastest warming on Earth [15]. An accompanied significant reduction in glacier area and mass balance through the 21st century is expected [16], which makes an investigation of recent elevation and mass change patterns indispensable. Simulations have already shown negative surface mass balance and mass loss since the 1980s at variable rates [15,17,18], which is supported by several remote-sensing based negative elevation and mass change estimates over the last decades [19][20][21]. Therefore, the aim of this study is to evaluate the potential of ICESat and ICESat-2 data in Svalbard to determine elevation and mass change trends in the period of 2003-2008 to 2019. The derived ICESat/ICESat-2 change trends are compared to recent CryoSat-2 results between 2011-2017 [19], and observed alterations relative to 2003-2008 intra-ICESat acquisitions [20] are discussed. The assessment of the used technique to derive ICESat/ICESat-2 results is provided through an in-depth analysis of the performance over Svalbard. In support of recent radar altimetry measurements [19], this use of precise laser acquisitions by ICESat/ICESat-2, with its unique properties, is expected to deliver detailed insights of elevation and mass change patterns in one of the most critical environments regarding climate change.

Study Area
Svalbard is an Arctic archipelago located north of Norway between 75-82 • N ( Figure 1). It is surrounded by the northern Atlantic Ocean to the south and the Arctic Ocean to the north. The western parts of the six major islands are mostly influenced by the warm West Spitsbergen Current, whereas the eastern part is impacted by clearly cooler waters from the Arctic Ocean [22]. Due to the location at the junction of northern cold and dry polar air masses and southwestern warm and humid air masses from the Atlantic currents, meteorological conditions are spatially and temporally versatile [20]. Overall, the archipelago is characterized by strong climate gradients, including milder and more humid conditions in the south and west, and colder and drier conditions in the northeast [23].  [20]. NW is northwest Spitsbergen, NE is northeast Spitsbergen, SS is South Spitsbergen, VF is Vestfonna, AF is Austfonna, BE is Barentsøya and Edgeøya, and KV is Kvitøya. Note that northern parts of SS are included here, but not in [20]. Red points indicate crossover locations of ICESat and ICESat-2 measurements used in this study. Background © Bing Satellite.
The ice masses in Svalbard consist mainly of ice caps, which cover most of the land mass of Nordaustlandet and Kvitøya, and extensive ice fields on the interiors of Spitsbergen and Barentsøya/Edgeøya. Svalbard covers a total area of~60,000 km 2 [23], from which 57% or~34,000 km 2 are covered with glaciers [24]. This corresponds to~10% of the glaciers in the Arctic outside the Greenland ice sheet [23] and represents 6% of the worldwide Glacier coverage outside Antarctica and Greenland [20]. Recent estimates about total ice volume agree on~6200 km 3 or a sea-level equivalent of 1.5 cm [25,26]. Of all glaciers, 15% by number [23] and up to 60% by area [27] are tidewater glaciers terminating into fjord or ocean water [23]. Overall, glaciers are mainly polythermal, thus consisting of cold and temperate ice [23]. They have typically shown slow velocities of 2-10 m yr −1 [28] and a negative surface-mass balance and mass loss at variable rates since the 1980s [15,17,18], supported by several remote-sensing based negative elevation and mass change estimates over the last decades [19][20][21]. Contrasting other Arctic glaciers, Svalbard's ice masses exhibit a low maximum of area-elevation distribution at~450 m a.s.l., compared to ice caps of, for example, Greenland, Arctic Canada, and Iceland, which have corresponding altitudes between~800-1400 m a.s.l. [15]. With a proportion of 60%, most of the ice is located under this rather low hypsometric peak. Additionally, they are very vulnerable to changes in atmospheric temperatures because the equilibrium-line-altitude fluctuates around that peak [15].
Surging behavior is common over the archipelago. Most reported surging glaciers are located in southern Spitsbergen (~45%), while the rest is equally distributed over the other regions [29]. The overall proportion of surging glaciers is uncertain and estimates to range from 13% to 90% of all glaciers [29,30]. Overall, 345 surge-type glaciers have been identified until 2013 [31]. The quiescence phase in Svalbard typically lasts between 30-500 [32] years and is followed by an active (surging) phase ranging between 3-10 years [30,32]. This is much longer than for surge-type glaciers in for example, North America, Iceland, and the Pamirs, where active phases only last 1 to 2 years [30].

ICESat and ICESat-2
ICESat was launched on 12 January 2003 and ended its operation after seven years [33]. The payload included the 'Geoscience Laser Altimeter System' (GLAS), which contained three 1064 nm laser altimeters with a laser pointing angle determination system and a 1064/532 nm cloud and aerosol LIDAR (light detection and ranging) [34]. While operating, each laser pulsed at 40 Hz and illuminated a spot of about 65 m in diameter on the earth's surface, whereas individual spots are separated by 172 m in along-track direction [3]. GLAS measurements were found to have a <3 cm shot-to-shot precision over ice sheets, which significantly improved the accuracy of elevation measurements over such surfaces [34]. On 29 March 2003 Laser 1 failed and caused an adaption to several month-long measurement campaigns [3,34]. This reduced the GLAS operation time by 73% [34]. Hence, during the whole mission, the adopted 33 to 56-day campaign plan was used to extend the ICESat mission lifetime [35].
Its successor, ICESat-2, was successfully launched on 15th September 2018 [36] and carries the Advanced Topographic Laser Altimetry System (ATLAS) instrument. ATLAS uses a photon-counting LIDAR instrument to measure the time a photon needs from ATLAS to the earth and back, the pointing vector at the time a photon is transmitted, and the position of ICESat-2 in space at the time a photon is recorded [37]. ATLAS transmits green (532 nm) laser pulses at 10 kHz and splits each laser pulse by a diffractive optical element to generate six individual beams in three pairs. The beams in a pair have different transmit energies (energy ration of~1:4) and are separated by 90 m in the across-track direction. Individual pairs are separated by~3.3 km in the across-track direction, and strong and weak beams are separated by~2.5 km in the along-track direction [37]. Through firing at 10 kHz from a nominal~500 km orbit height, a laser pulse is transmitted everỹ 0.7 m in along-track direction for each beam [37]. This ensures minimal measurement gaps and provides a higher fidelity of topography, even over rough and heterogeneous glacier and sea surfaces [38]. With this high sampling rate, a narrow footprint of~14.5 m, and a near global coverage (+/−88 • latitude), the mission is designed to repeatedly measure the ice sheets by revisiting 1387 reference ground tracks with a 10 m horizontal accuracy and 0.03 m vertical precision every 91 days [11].

Materials and Methods
The extraction of change signals from ICESat laser altimetry is typically performed by many different techniques (e.g., [20,[39][40][41][42]). In this study, orbit crossover locations of the ICESat and ICESat-2 missions were used to compare height values measured by both systems at those positions. Some processing steps were thereby performed using NASA's cryosphere altimetry processing toolkit 'captoolkit' by F. Paolo, J. Nilsson, A. Gardner, and T. Sutterly (https://github.com/fspaolo/captoolkit (access on 1 April 2020)). This tool provides a collection of editable Python scripts for many purposes regarding laser and radar altimetry. The crossover method is additionally tested in comparison to common intra-ICESat repeat-track analysis over a test ground track in Vestfonna, which is provided in the Supplementary Materials.

Data Preparation
For ICESat, the GLAS/ICESat L1B Global Elevation Data GLAH06 Release 634 of the National Snow & Ice Data Center (NSIDC) was used between 2003 and 2008. Some corrections were applied to the GLAH06 data, including a detector saturation correction [43] and an inter-laser bias correction for ICESat's laser 2 and 3 [44]. Together with those major corrections, data culling, and general corrections for ICESat data were carried out as follows: • Sorted out data not meeting the on-product variable criteria: 'elev_use_flg == 0', 'sat_corr_flg < = 2', 'sigma_att_flag == 0', and 'i_numPk == 1'; • Applied the saturation correction; • Converted from the Topex/Poseidon ellipsoidal height to WGS84 ellipsoidal height; • Applied the inter-laser correction: subtracting 1.7 cm from Laser 2 and adding 1.1 cm to Laser 3 according to [11] due to the same methodological approach.
Data preparation was achieved for all ICESat high energy spring campaigns (L1A, L2B, L3B, L3E, L3H, L3J), which makes a total of six campaigns from 2003 to 2008. Only spring campaigns were chosen to minimize seasonal effects.
ICESat-2 data came from the ATLAS/ICESat-2 L3A Land Ice Height ATL06 Version 3 dataset by NSIDC. ATL06 data were produced with a special segmentation technique: along-track geolocated ATL03 photon data from each beam were divided into 40 m long overlapping (20 m) segments, followed by fitting photon heights with a linear model as a function of along-track distance [11]. Data were finally corrected for instrument specific deviations, such as first-photon bias and transmit-pulse shape, which produced the final ATL06 product. This includes latitude, longitude, and height, with respect to the WGS84 ellipsoid for each segment [11]. ATL06 data from 14th January 2019 to 14th April 2019 were used to ensure coverage of the complete ICESat spring campaign months while using all measurements of a full 91-day repeat cycle. Data were filtered similar to [11], as follows:

•
Removed data flagged by the on-product ATL06 quality summary value; • Removed data where adjacent segment height differences were >2 m and therefore exhibited high along-track height variability; • Removed data where surface height was >10,000 m due to atmospheric scattering.

ICESat-1/2 Crossover Analysis
Crossover locations were determined between all six ICESat spring campaigns and the 2019 ICESat-2 spring dataset. Due to the 6-beam architecture of ATLAS, crossing of one ground track from ICESat with the tracks of ICESat-2 produced up to six crossovers. To calculate orbit intersections, latitude and longitude values projected in UTM Zone 35N for all footprints of always one ICESat-2 track were inserted into geometric functions with all measurements from always one track of ICESat ( Figure 2). This was repeated until each ICESat-2 track was compared with every ICESat track. The algorithm was implemented by the authors of the 'captoolkit' and produced the exact location of intersections between all used ICESat/ICESat-2 tracks. After a crossover location was determined, it was checked if one of the four nearest footprints exceeded the threshold of 350 m distance to the crossing position. If so, the crossover value was not used in the analysis. Because of the larger distance between individual ICESat footprints, this allowed for one ICESat measurement to be missing, even if the neighboring footprint was close to the orbit intersection. Next, the height data of the four measurements were linearly interpolated to the crossing position and both interpolated elevation values were subtracted from each other. dh/dt values were finally calculated by dividing the absolute elevation change values through the timespan between the two acquisitions (between~11 to~16 years). The calculation of such annual elevation change time series (dh/dt) allows for the comparison of elevation change values between the six ICESat campaigns and the ICESat-2 acquisitions. In total, 8822 crossovers were found between all datasets. The crossing spot is not necessarily located in between the two closest footprints of each system.

Hypsometric Extrapolation
Despite achieving up to six elevation change measurements between one orbit crossing, the spatial extent of the crossovers is still mostly limited to the ICESat reference tracks. Since data coverage is therefore too sparse to determine total change rates through local spatial interpolation, especially in mountainous terrain in NW and NE [20], a hypsometric approach similar to [20,21], was used to estimate ice volume and mass changes. At first, the data were separated into the different regions presented in Figure 1. Next, dh/dt values were split by the surface heights at the crossover locations into 50 m elevation bins for every region. All dh/dt values for each bin were then filtered by excluding all values which exceed the threshold of three times the mean absolute deviation from the median of all values. After filtering, 8021 valid crossover points were left, excluding general outliers and some extreme values through surging events. To derive annual volume changes, 50 m glacier hypsometries were determined by the 20 m resolution 'S0_DTM20_NP-ArcticDEM-Mosaic' Svalbard digital elevation model from the Norwegian Polar Institute (NPI) [45]. The DEM is an updated version of the 'S0_Terrengmodell' from 2014 and includes more recent data [46]. The original DEM was mostly processed using stereo models from aerial photos and from elevation contours, lakes, and coastlines [47]. Data from the 'Randolph Glacier Inventory 6.0' (RGI) [48] served as a mask to clip the DEM to only glaciated regions. Glacier areas were calculated for each region by counting the cells in the model per 50 m bin (e.g., in 100-150 m altitude) and multiplying the number of cells by their median size. A third order polynomial curve was fit to the mean of all dh/dt values for every elevation bin per region ( Figure 3). This is required to fill bins which are lacking dh/dt measurements but have glacier area values.
Extrapolation of measurements to unsurveyed areas with polynomial fits are a common practice in laser (e.g., [20,21]) and radar (e.g., [19]) investigations. In this context, the usage of third order polynomial fits in this study allows for a direct comparison with the hypsometries derived from ICESat by Moholdt et al. [20] and CryoSat-2 by Morris et al. [19]. The third order polynomial fit can produce runaway tails at the data edge, however, as bins which are filled with the polynomial value typically correspond to just one bin with an area proportion ranging between 0.11-0.99% of the total regional area value, their contribution to volume change calculations is neglectable. Storisstraumen is an exception with 8.77% of total area filled with polynomial fits, which is still expected to have just a minor impact on calculations. Using the mean dh/dt value with fits per bin or just polynomial values resulted in no significant difference in final results. For example, the total Svalbard mass loss rate (dM/dt 2 ) calculated with mean dh/dt and fit values was −12.21 ± 4.28 Gt yr −1 , compared to −12.19 ± 4.28 Gt yr −1 , by using just polynomial fit values.
Annual volume change (dV/dt in m 3 ) for each region was finally calculated using the following equation: where n is the number of 50 m bins per region, h is the mean of all dh/dt measurements in m yr −1 , or (if missing) the polynomial fit, and A is the glacier area in m 2 per bin i. Averaged elevation change rates dh/dt were then calculated by dividing the dV/dt value of each region by the corresponding total Glacier area A. Further, mass loss rates were determined by multiplying dV/dt with the density of ice (0.917 g/cm 3 , dM/dt 1 ) and additionally with a density of 0.850 g/cm 3 (dM/dt 2 ) to partly include changes in the firn pack. The conversion factor for dM/dt 2 was chosen based on the suggestions of Huss et al. [49] for periods > 5 years by assuming stable mass balance gradients, the presence of a firn area, and volume changes significantly other than zero. Due to the large spatial extent and extreme dh/dt values through a surging event [50], acquisitions over the Storisstraumen Glacier in Austfonna were treated separately. Values for Austfonna were subsequently estimated without the contribution from the glacier basin, whereas mass loss contribution from the surge is discussed in Section 6.2. The final Svalbard-wide mass change value is the sum of all regional dM/dt 2 estimates and the dM/dt 1 value for Storisstraumen, assuming that the mass loss from this basin was mainly caused by surging ice. The total thinning rate was determined by using all regional filtered crossover values excluding Storisstraumen and corresponding total bin area values. Further, despite the spatial relative limited data, ordinary kriging interpolation was performed with all crossover measurements (including Storisstraumen) to produce an approximate spatial overview of elevation change trends over Svalbard (Figure 4). For this purpose, a grid with a 1 km 2 resolution was roughly cropped to the outline of the RGI, and an experimental variogram based on all crossover values was computed. Next, a spherical model with a range of~30, sill of~1.1, and nugget of~0.05 as parameter settings was fitted to the variogram, and all dh/dt values were finally interpolated over the grid, based on this model.

Error Assessment
Over non-glaciated areas, such as rocky, solid terrain, the surface elevation was not expected to change over time. Hence, dh/dt values at ICESat/ICESat-2 orbit intersections over stable ground were expected to be zero. With this assumption, the crossover-error ∂ cross was estimated by crossing all ICESat campaigns with the ICESat-2 data over ice-free ground. This error assessment approach over non-glaciated terrain is common in geodetic glacier mass balance approaches and was used by [21] over Svalbard with ICESat data. For every region, ∂ cross is determined by calculating the root mean square error (RMSE) of all calculated values (dh/dt at crossover locations over ice-free ground) relative to the assumed elevation change of 0 m yr −1 through the following equation: where n and N are the number of crossovers over ice-free ground for every region, and z i is the dh/dt value in m yr −1 at a crossover location i. After filtering all region-wide values with the threshold defined in Section 4.3, 3870 crossovers over stable ground produced error estimates ranging from 0.09 m yr −1 in Austfonna to 0.21 m yr −1 in northwest Spitsbergen. The Svalbard-wide error is 0.15 m yr −1 by calculating the area weighted mean value of all regional error estimates. Due to the very sparse non-glaciated area in Kvitøya, no crossovers on stable ground exist, and the error for Austfonna is applied in further processing steps. For Storisstraumen, the error of Austfonna was likewise applied. A previous Svalbard-wide mean crossover error of 0.20 m yr −1 was calculated within individual ICESat campaigns over ice where again only small elevation changes were expected [20]. The ICESat/ICESat-2 crossover error was consequently smaller in all regions except northwest Spitsbergen and reduced for the whole of Svalbard compared to this intra-ICESat crossover error. Regional and Svalbard-wide volume error estimates were determined by multiplying the respective ∂ cross values with the corresponding glaciated area values, and mass change errors were further calculated by applying the conversion factors for dM/dt 1 and dM/dt 2 .  Table 1, regional and Svalbard-wide elevation, volume, and mass change values are presented. The estimated Svalbard-wide dh/dt value is −0.30 ± 0.15 m yr −1 , including highest rates at South Spitsbergen with −0.80 ± 0.18 m yr −1 and lowest elevation change at Austfonna and northeast Spitsbergen with −0.07 ± 0.09 m yr −1 and −0.07 ± 0.16 m yr −1 , respectively (Table 1). At low altitudes (<300 m), frontal thinning is observable over all regions in the magnitude of~0.5-5 m yr −1 . Vestfonna on the other hand, experienced only slight low elevation thinning. At higher altitudes, dh/dt values tend to get more positive, but with variations, from diminished thinning at, for example, South Spitsbergen, to thickening in northeast Spitsbergen and Austfonna up to~0.5 m yr −1 . Despite the elevation decrease at the margins, this thickening trend is clearly evident in the interior Glacier and ice cap areas of Austfonna and northeast Spitsbergen. By excluding Kvitøya, a Southwest-northeast gradient is visible, specified by pronounced thinning in northwest Spitsbergen, South Spitsbergen, Barentsøya/Edgeøya, and less elevation decrease in northeast Spitsbergen, Austfonna, and Vestfonna. Overall, all regions show negative total dh/dt values, although error consideration could lead to estimates of balance for Austfonna, Vestfonna, and northeast Spitsbergen. Extreme thinning rates are especially notable (up to 10 m yr −1 ) over Storisstraumen in Austfonna, and are related to a surging event during the observation period [50].

Regional Elevation Change
Northwest Spitsbergen has an overall high annual thinning rate of 0.63 ± 0.21 m yr −1 , in addition to the significant thinning rate of~2 m yr −1 in low altitudes (Figure 3). Only slight thickening <0.5 m yr −1 is measured at higher altitudes >700 m in the interior, where Glacier area is small. The highest negative dh/dt values can be found at the calving fronts of marine terminating glaciers, such as in the Trollheimen region in the south or Albert-I-Land in the north (Figure 4).
Northeast Spitsbergen exhibits much less elevation decrease and additional thickening at elevations >350 m up to~0.5 m yr −1 . The overall negative dh/dt estimate of −0.07 ± 0.16 m yr −1 for the whole region results from thinning rates as high as~5 m yr −1 in combination with relatively large glacier area distribution at low altitudes. In contrast, thickening in higher altitudes is taking place at a much smaller magnitude. Negative values are spatially mostly present at the margins, and highest decrease rates can be found e.g., around the marine-terminating Oslobreen and Kongsfonna/Hachstetterbreen Glacier in the east.
With 0.80 ± 0.18 m yr −1 , South Spitsbergen experienced the highest thinning rates of the whole archipelago, including an elevation decrease of up to~−4 m yr −1 at low altitudes and overall thinning at almost all elevations. Exceptions of minor thickening are only measured at some elevations >700 m. Hotspots of decrease are found in the Wedel Jarlsberg Land region in the south, especially for example, at the Nathorstbreen Glacier system and around the Strongbreen in the east.
At the Barentsøya/Edgeøya islands, measured thinning of up to~3 m yr −1 is found at elevations high as 400 m. Glacier area is mostly distributed up to this threshold, which results in a high regional overall decrease value of −0.54 ± 0.13 m yr −1 . Slight thickening in higher altitudes is mainly found in the interior of both islands.
Next to northeast Spitsbergen, Austfonna experiences the smallest overall thinning rate of −0.07 ± 0.09 m yr −1 , but with thinning up to almost 3 m yr −1 at elevations under 350 m and thickening of up to~0.5 m yr −1 above. Positive dh/dt values can be found in the interior of the ice cap, while thinning is limited to the margins, especially in the east and south. As stated before, extreme dh/dt values of up to −10 m yr −1 are measured over Storisstraumen, thus distorting the overall pattern of thickening in the interior and thinning at the edges.
Vestfonna also shows a small elevation decrease rate of −0.09 ± 0.12 m yr −1 , including the least losses or minor thickening at elevations with the largest glacier area. In the north, thickening takes places in contrast to slight thinning up to~0.5 m yr −1 in the south.
The island Kvitøya in the northeast exhibits thinning over the complete ice cap, pronounced at the margins, which is specified by elevation decrease rates of up to~3 m yr −1 at low altitudes. As the only region without any measured thickening trend, Kvitøya represents the area with the second highest thinning rate of 0.74 ± 0.09 m yr −1 .

Mass Change
Average thinning over the archipelago is taking place with rates up to~2.5 m yr −1 at elevations <450 m, where most of the glaciated area in Svalbard can be found ( Figure 3). Glacier area decreases significantly at higher altitudes in addition to relative moderate thickening rates up to~0.5 m yr −1 . This leads to an overall negative glacier elevation decrease rate for all of Svalbard. Patterns of mass loss are reflected by elevation change in relation to glacier area. South Spitsbergen has the highest mass loss rate of −3.91 ± 0.87 Gt yr −1 , followed by northwest Spitsbergen with −3.35 ± 1.11 Gt yr −1 , and Barentsøya/Edgeøya with −1.06 ± 0.25 Gt yr −1 (dM/dt 2 ). With −0.47 ± 1.11 Gt yr −1 and −0.40 ± 0.54 Gt yr −1 , northeast Spitsbergen and Austfonna experience smaller loss rates due to lower overall thinning rates. Kvitøya shows a similar loss rate of −0.40 ± 0.05 Gt yr −1 despite significant higher thinning rates, a result of the small Glacier area. Vestfonna has the smallest mass loss rate of −0.19 ± 0.24 Gt yr −1 because thinning rates are overall low. In addition to the non-surging mass loss rate of Austfonna, a loss rate of −2.62 ± 0.10 Gt yr −1 (dM/dt 1 ) is obtained for the surge of Storisstraumen. Eventually, the overall estimated mass loss rate in the period 2003-2008 to 2019 for all of Svalbard is −12.40 ± 4.28 Gt yr −1 .

Elevation and Mass Change
Different estimates of Svalbard-wide elevation and mass change rates are needed to compare and verify change signals derived from the ICESat/ICESat-2 crossover method. A previous study [20] used two different ICESat repeat-track methods to determine 2003 -2008 dh/dt values, while another work [21] compared ICESat measurements from 2003-2007 with older topographic maps and DEMs in varying periods between 1965-1990. Additionally, newer investigations over the archipelago were carried out by using CryoSat-2 radar altimetry between 2011-2017 [19]. Comparison of regional values can be difficult, because measurement timespans, spatial coverage, calculation of total change values (e.g., polynomial fits or mean values), and the handling of surges may differ and thus distort elevation and mass change rates. Nevertheless, general regional and Svalbard-wide change patterns between the different acquisition periods can be compared sufficiently. A detailed regional analysis of mass change patterns between the different investigations during the last decades and, additionally, regionally corrected and converted mass change rates from the previous studies [20,21], has already been provided in [19]. For a better overview, the Svalbard-wide regional mass change estimates [19][20][21] in relation to results of this work are shown in Table 2.  [19]. 2 Values are those of non-surging and surging ice combined, AF is just non-surging ice and Storisstraumen is surging ice. KV is included in AF. 3 Value is without Storisstraumen.
By representing more recent results, outcomings from the ICESat/ICESat-2 comparison largely supports the findings of the CryoSat-2 measurements [19]. Similar to the radar acquisitions, a pronounced thinning trend at the margins of the archipelago and decreasing change rates in higher altitudes are observable. While areas of thickening are small and restricted to the high portions of the interior ice fields in the south and northwestern regions (SS, BE, NW), large inland areas of the ice fields in northeast Spitsbergen and the ice caps of Austfonna and Vestfonna are either thickening or in balance. Further, thinning in the northeastern regions is mostly limited to lower altitudes <350 m. The general pattern of elevation change is interrupted by more drastic change rates resulting from surges, which were mostly incompletely measured by the crossover attempt. Nevertheless, the general south-northwest to northeast gradient of elevation change over Svalbard is most consistently observable in the CryoSat-2 [19] and ICESat/ICESat-2 results. ICESat/ICESat-2 acquisitions additionally show higher thinning rates at elevations <100 m in northwest Spitsbergen, northeast Spitsbergen, and Austfonna, which is expected to result from the smaller ICESat/ICESat-2 footprint compared to the CryoSat-2 radar sensor. In this context, smaller elevation decrease rates at altitudes 0-50 m compared to 50-100 m, as observed in northeast Spitsbergen and Barentsøya/Edgeøya, are assumed to result from partially floating glacier tongues of calving glaciers and the frontal recession of land terminating glaciers.
Compared to the intra-ICESat acquisitions in 2003-2008 [20], the ICESat/ICESat-2 measurements indicate more negative elevation and mass change rates overall, except for Vestfonna (Table 2). This is specified by a generally more negative trend in almost all regions, in combination with increased thinning rates at low elevations (Figure 3). Acceler-ated elevation decrease at the margins is especially observable in northeast Spitsbergen, South Spitsbergen, Barentsøya/Edgeøya, and Kvitøya, which eventually leads to an increased thinning rate of~1 m yr −1 at elevations <100 m for all of Svalbard, relative to the intra-ICESat measurements [20]. Furthermore, a switch from thickening to thinning of~1 m yr −1 in east South Spitsbergen (Figure 4) and decreasing thickening trends in Austfonna are observable. The CryoSat-2 results likewise show this switch to pronounced thinning in the last decade, compared to the ICESat period of 2003-2009 [19]. However, the ICESat/ICESat-2 hypsometry for VF resembles that of the ICESat ones [20], whereas the CryoSat-2 measurements show distinct negative elevation change rates at elevations <300 m [19]. The authors of the CryoSat-2 investigation state that the ICESat hypsometry for VF was affected by a raised surface of the western Franklinbreen glacier during a surge [19]. Therefore, the contrasting characteristic of the low-elevation trend in VF can be explained through the inclusion of the ICESat period in the ICESat/ICESat-2 comparison, the positive dh/dt trend of crossover measurements due to their location at the marine terminating front of that catchment, and the significant differences in surface elevation change over the Franklinbreen glacier in the period 2000-2017 relative to the 2011-2017 CryoSat-2 measurements [51]. Nevertheless, the laser altimetry acquisitions and the findings of the recent CryoSat-2 research [19] underline the switch to increased thinning and mass loss since the ICESat measurements in 2003-2008 [20]. In this context, differences in total regional change rates derived by ICESat/ICESat-2 and CryoSat-2 [19] are expected to stem from the different acquisition timespans, the smaller footprint of ICESat and ICESat-2 relative to the CryoSat-2 radar sensor and corresponding coverage variances, and the different presence and handling of surges during the investigation periods.
By looking at Svalbard-wide values, mass change rates for the complete archipelago increased from −3.40 ± 1.60 Gt yr −1 from 2003-2008 to −16.00 ± 3.00 Gt yr −1 in 2011-2017 [19]. This compares to the mass loss rate of -12.40 ± 4.28 Gt yr −1 of the ICESat/ICESat-2 approach from 2003-2008 to 2019. Considering the larger timespan and inclusion of the ICESat period with smaller mass loss rates in the ICESat/ICESat-2 comparison, the slightly reduced loss rate of the latter agrees well with the change rate derived by CryoSat-2 [19].  [52] additionally underlines the switch from less mass loss during the ICESat period to increased loss rates in the last decade.

Contribution of Surges
Overall, an estimation of the contribution from surging ice through ICESat/ICESat-2 is difficult because limited distribution of orbit intersections prevents sufficient coverage of most individual glacier systems. Mass loss from surging events in Spitsbergen (NE, NW, and SS) is expected to be small due to the relatively small area of surging ice in comparison to non-surging ice and the redistributive nature of surges [19]. However, change patterns from two significant surges of the Stonebreen in Edgeøya [53] and the Nathorstbreen in South Spitsbergen [54] may be underrepresented or distort the final change rates of those regions. The contribution to mass loss from those surges reported by Morris et al. [19] is −0. The overall negative elevation and mass loss trend over Svalbard can be directly linked to Arctic warming and associated changes in glacier climatic mass balance (CMB). A general warming of 3-5 • C between 1971-2017 has been observed over the archipelago, with the largest increase in winter and the smallest in summer [16]. In a warming climate, melting enhances over the glaciers surface. The polythermal characteristic of Svalbard's glaciers provide a significant retention capacity to refreeze large amounts of meltwater in porous snow and firn [23]. Climatic simulations indicate that the average summer melt of Svalbard's ice caps exceeded annual total precipitation by 25% before a modest warming trend in the mid-1980s triggered a net mass loss. This underlines the importance of meltwater refreezing in the firn to sustain the ice caps [15]. Different simulations and model, covering the last decades consistently show a reduced firn area and refreezing capacities in combination with elevated equilibrium-line-altitudes (ELA), due to atmospheric warming [15,17,18]. Recent simulations report a fast ablation zone expansion from 27% to 44%, due to an ELA upward movement of~100 m to the hypsometric peak of~450 m in combination with a reduced firn refreezing capacity from 54% in 1958-1984 to 41% in 1985-2018, which resulted in enhanced meltwater runoff by 55% since the mid-1980s [15]. This increased melt could also have accelerated the triggering effect of surges in Svalbard [50].
Contrary to the general negative trend, an observed and modeled mass loss pause between 2005-2012 [15] is believed to be caused by changes in atmospheric circulations [55]. The mean 1979-2005 circulation over Svalbard was a west-southwesterly flow, which changed after 2005 due to more frequent North Atlantic Oscillation negative phases in summer. Therefore, the post-2005 induced northwesterly flow over Svalbard brought colder air over the archipelago, which partly counteracted Arctic warming and explains the mass loss stop between 2005-2012 [55]. This mass loss pause was again replaced by reinforced loss rates due to intensified ice discharge after 2012 [15], resulting, for example, from the surge of Storisstraumen [50]. Accompanied record melt and summer temperatures in 2013 are found to be caused by renewed changes in atmospheric circulations from the northwesterly to again more westerly flows over Svalbard [55]. Overall, these circumstances are also expected to be responsible for the overall lower mass loss trend measured during the ICESat period 2003-2008 [20].

Climatic and Oceanic Drivers of Change
In addition to the general evaluation of changes in climatic mass balance, a direct investigation regarding regional climatic and oceanic changes in the ICESat/ICESat-2 period is needed to explain observed change patterns. Morris et al. [19] used microwave radiometry and ERA5 climate reanalysis data over the Barents Sea region to investigate possible links with increased mass loss during the ICESat 2003-2008 and CryoSat-2 2011-2017 periods, which covers the most part of the ICESat/ICESat-2 timespan. A general sea surface temperature (SST) increasing trend over the region is specified by about 1 • C warming, along the western coast of Spitsbergen, and about 1.5-2 • C in Storfjorden between Spitsbergen and Barentsøya/Edgeøya. 2 m air temperature followed a similar pattern at a smaller magnitude of +0.75 • C and +1 • C, respectively. The changes in sea and air temperatures are accompanied by a decline in sea ice, especially around South Spitsbergen, Barentsøya/Edgeøya, Hinlopen Strait, and the south of Austfonna, in addition to little summertime sea ice in both periods, along the west coast of Spitsbergen [19]. This pronounced change pattern at the eastern side is reflected in measured elevation change rates by the ICESat/ICESat-2 investigation, as increased mass loss comes mainly from the eastern parts of Svalbard. This is specified by changes from slight mass gain in eastern South Spitsbergen measured by ICESat in 2003-2008 [20] to recent mass loss, which shows similarities with the amplified atmospheric and SST warming and the sea ice decline over Storfjorden. Observed increased thinning and mass loss in Barentsøya/Edgeøya and Olav-V-Land (around the Negribreen glacier) in south northeast Spitsbergen likewise reflect the climatic and oceanic changes in that region. Similar conditions can be seen at the southern parts of Austfonna and Vestfonna, which border the Hinlopen Strait, and Kvitøya, where mass loss and thinning is expected to be caused by warming in the atmosphere and sea. Compared to the ICESat measurements from 2003-2008 [20], the ICESat/ICESat-2 results consequently indicate a spread of pronounced thinning and mass loss from the west coast of Spitsbergen into the western Barents Sea regions, thus agreeing with recent findings of Morris et al. [19]. A pronounced regional warming maximum has been observed in the northern Barents and Kara Seas between 75-85 • N and 39-90 • E [14]. Diminished sea-ice import and a corresponding loss in freshwater since the mid-2000s led to a weakened ocean stratification and therefore to amplified vertical mixing and increased upward fluxes of heat and salt. This prevents sea-ice formation and increases ocean heat content, which could soon transform the northern Barents Sea from a cold and stratified Arctic climate regime to a warm and well-mixed Atlantic-dominated climate regime [56]. Next to Svalbard, other regions bordering the Barents Sea responded similarly to this warming trend. Franz-Josef-Land has experienced an acceleration of mass loss since the late 2000s, including significant thinning at the southwest parts, which is linked to prevailing warmer conditions [57]. Glaciers along the Barents Sea coast in Novaya Zemlya show pronounced mass loss rates in combination with high retreat rates at marine-terminating glaciers since 2000 [58]. Furthermore, recent research reveals a doubling of glacier mass loss on the Russian archipelagos compared to measurements before 2010 [59]. Regional warming and correlated amplified thinning at the eastern parts of Svalbard is consequently also reflected in changes of other glaciated areas bordering the Barents Sea region.

Crossover Technique Performance
The crossing of ICESat and ICESat-2 data results in 8021 verified dh/dt values. A previous investigation [20] used 329 crossovers over ice within individual ICESat campaigns in Svalbard to estimate an intra-ICESat crossover error. By assuming that this number of measurements resulted from all ICESat campaigns and all points were used, the increase in crossovers from ICESat/ICESat-2 is over 24 times higher by only using the seven ICESat spring campaigns compared to the number of intra-ICESat crossovers. Therefore, the enhanced coverage and sampling rate of ICESat-2 in combination with the 6-beam architecture resulted in a drastic increase in measurements. Although this represents a significant improvement, the crossover locations are still mostly limited to the rather spare ICESat reference tracks (Figure 1). This prohibits extensive coverage, which makes hypsometric extrapolations to derive total change values over larger regions necessary. Furthermore, the six beams in combination with the offset of repeating ICESat campaign tracks result in many 'point clouds', which can inherit many crossing locations in close proximity. Consequently, the large number of crossovers is mostly found in groups along the ICESat reference tracks. Besides the necessity for some sort of extrapolation, this limited coverage is especially problematic over mountainous, complex terrain, such as in South Spitsbergen. On a regional scale, the limited spatial coverage of measurements and their grouping in clusters leads to an unequal distribution, where certain individual glacier basins and altitudes experience extensive observation, while others have few or no values at all. This is especially problematic over surging basins, where crossing positions may only cover parts of the basin and therefore not represent the complete change pattern. Typical measurement coverage and the grouping of values in 'point clouds' is shown in Figure 5 by using the example of the Nathorstbreen glacier system. Crossover locations are spatially sparse and limited to the middle-upper parts of the basin, which prevents an observation of the complete regional basin-wide change pattern. Similar conditions and resulting problems were found e.g., over the Franklinbreen glacier in Vestfonna, as discussed in Section 6.1. Earlier research with ICESat data over glaciers likewise reported similar problems regarding spatial coverage [9,20]. In addition to limited measurement coverage, the~172 m spacing of neighboring ICESat footprints and the maximal 350 m spacing to the crossover locations may be too large in mountainous, complex terrain, including many small glaciers, and therefore distort actual change trends. Setting lower thresholds would however significantly reduce measurements in those areas, leaving many regions unobserved. The values over the ice caps of Austfonna, Vestfonna, Kvitøya, and the interior ice fields of northwest and northeast Spitsbergen are consequently expected to be more accurate than those over alpine, rough terrain, especially in South Spitsbergen. Next to complex topography, biases can be additionally induced due to the 2003-2008 timespan of ICESat measurements, which can contain changing patterns themselves. Neighboring ICESat footprints over areas with large changes (e.g., surges in 2003-2008) may exhibit significantly different elevation rates, thus resulting in overall heterogeneous dh/dt values when compared with recent ICESat-2 data from 2019.
Beside these shortcomings, crossovers over stable ground showed little deviations, thus indicating a generally high accuracy of derived dh/dt values. This is specified by a Svalbard-wide elevation change error reduction of 25% compared to intra-ICESat crossovers [20], although comparison could be biased, due to the differences of used terrain (stable ground/ice). Improvements through the inclusion of ICESat-2 data compared to intra-ICESat analysis is likewise supported by findings of recent research with the ICESat/ICESat-2 crossover method in Greenland and Antarctica [11]. Investigation of the performance of the six ICESat-2 beams in comparison to the ICESat beam through intra-beam crossovers within a timespan of less than 30 days showed that ICESat crossover errors are higher and increase by a factor of 2.5 over a range of cross track slopes, while ICESat-2 beams comparably perform with lower overall error values [11]. Additionally, a new "density-dimension algorithm" transforms ATLAS acquisitions to ice height retrievals with a spacing of 0.7 m and automatically adopts measurements to variable background conditions [60]. This enables the determination of detailed surface characteristics, such as those on heavily crevassed, rapidly accelerating surge-type glaciers [60]. Technological and methodological improvements, such as those are expected to produce highly accurate elevation change measurements in further crossover attempts by including ICESat-2 data, even in complex terrain.
Eventually, surface height acquisitions along the ICESat reference tracks have shown to be sufficiently able to estimate region-wide elevation change values when extrapolated over larger regions (e.g., [20,21]). In this context, and by considering the different time spans and regional classification deviations, the ICESat/ICESat-2 change signals are in good agreement with the CryoSat-2 estimates [19]. This proves the suitability of also using the crossover technique over a smaller glaciated region like Svalbard instead of the large ice sheets of Greenland and Antarctica, if some restrictions, as described above, are considered.

Conclusions
This study uses laser altimetry data by ICESat from spring 2003-2008 and new measurements by ICESat-2 from spring 2019 to investigate annual elevation and mass change trends over Svalbard in this period. Orbit crossover locations between both missions produce 8021 change signals, which are regionally extrapolated with a hypsometric approach, including polynomial fits. The elevation change between 2003-2008 to 2019 exhibits a negative regime over all regions, specified by a total Svalbard elevation decrease rate of −0.30 ± 0.15 m yr −1 (−12.40 ± 4.28 Gt yr −1 ). A pronounced spatial thinning trend is observable, including the highest rates at the southwest regions (NW, SS, BE) and less elevation decrease rates at the northwest regions (NE, VF, AF). Hotspots of thinning are found at the margins of the archipelago and at some marine-terminating and surging glaciers. Areas of thickening are mostly limited to northern Vestfonna, the interior and northern parts of Austfonna, and the interior ice fields of northeast Spitsbergen. Findings of the ICESat/ICESat-2 comparison mostly agree with recent Svalbard-wide CryoSat-2 measurements in 2011-2017 [19] and likewise show increased elevation decrease rates in almost all regions, as well as a spread of pronounced thinning from the west coast to areas bordering the Barents-Sea compared to intra-ICESat acquisitions in 2003-2008 [20]. The overall negative elevation change regime in Svalbard is expected to be caused by Arctic warming in the last decades and corresponding changes in glacier climatic mass balance, especially through a decreased refreezing capacity since the mid-1980s [15,17,18]. Further, amplified thinning and mass loss since the ICESat period is assumed to correlate with a rise in atmospheric and sea-surfaces temperatures around the archipelago, especially along the west coast, Storfjorden, and Hinlopen Strait [19]. This trend is likewise observable in other glaciated regions bordering the Barents Sea [57][58][59] and is probably caused by recent oceanic transitions in this part of the Arctic [56].
The ICESat/ICESat-2 crossover approach suffers from technical shortcomings of the first ICESat mission and reveals some limitations regarding spatial coverage, especially in mountainous, complex terrain. The performance nevertheless exhibits improvements regarding accuracy and sample size compared to intra-ICESat crossover attempts [20]. Moreover, derived ICESat/ICESat-2 results largely agree with recent research using other data sources [19,52]. This proves the suitability of the crossover attempt for accurately estimating region-wide glacier mass change throughout Svalbard. Since this study represents a first attempt to compare data from ICESat and ICESat-2 over a smaller glaciated region with the crossover technique, further investigation of such small-scale areas is needed to extend knowledge regarding limitations and performance of this approach with data from both missions. This is especially important in non-polar glaciated regions, where spatial coverage of ICESat is even smaller. Overall, future use of intra-ICESat-2 data is expected to drastically improve coverage and accuracy of ice change measurements. This could be especially beneficial for small-scale research objects like, for example, the investigation of smaller surging basins and their contribution to mass loss. Recent research already showed promising results by including ICESat-2 data, and future investigations are expected to deliver highly detailed insights into a changing cryosphere.  Figure S3. dh/dt rates from the repeat-track method along the 2008 ICESat ground track in Vestfonna. Background © Bing Satellite; Figure S4. dh/dt rates from the crossover method along the 2008 ICESat ground track in Vestfonna. Background © Bing Satellite.
Author Contributions: L.S. processed the data, created the graphs, and wrote the paper. T.S. supervised the study and contributed to processing, data evaluation, and writing. M.H.B. initiated and led the study. All authors have read and agreed to the published version of the manuscript.
Funding: This work was partially supported by the FAU Emerging Fields Seed Grant TAPE.
Institutional Review Board Statement: Not applicable.

Informed Consent Statement: Not applicable.
Data Availability Statement: Regional filtered crossover values for glaciated and ice-free regions are available at AWI PANGAEA [link generation in progress].