Edinburgh Research Explorer Ice loss in High Mountain Asia and the Gulf of Alaska observed by CryoSat-2 swath altimetry between 2010 and 2019

. Glaciers and ice caps are currently the largest non-steric contributor to sea level rise, contributing ~30% to sea level budget. Global monitoring of these regions remains a challenging task since global estimates rely on a variety of observations 10 and models to achieve the required spatial and temporal coverage, and significant differences remain between current estimates. Here we report the first application of a novel approach to retrieve spatially-resolved elevation and mass change from Radar Altimetry over entire mountain glaciers areas. We apply interferometric swath altimetry to CryoSat-2 data acquired between 2010 and 2019 over High Mountain Asia (HMA) and in the Gulf of Alaska (GoA). In addition, we extract monthly time series of elevation change, exploiting CryoSat’s high temporal repeat, to reveal seasonal and multiannual variation in 15 rates of glaciers’ thinning at unprecedented spatial detail. We find that during this period, HMA and GoA have lost an average of –27.9 ± 2.4 Gt yr –1 (–0.29 ± 0.03 m w.e. yr –1 ) and –76.3 ± 5.6 Gt yr –1 (–0.89 ± 0.07 m w.e. yr –1 ) respectively, corresponding to a contribution to sea level rise of 0.048 ± 0.004 mm yr -1 and 0.217 ± 0.015 mm yr -1 . Glacier thinning is ubiquitous except for the Karakoram-Kunlun region experiencing stable or slightly positive mass balance. In the GoA region, the intensity of thinning varies spatially 13 Negative mass trends are also observed in all of the sub-regions in the GoA region, with the largest mass losses in the Coast Ranges and the Saint Elias Mountains. The GoA time series reveal an increased mass loss from 2014, most pronounced in subregions along the south-central and south-east coast (Saint Elias Mountains, Chugach mountains and Coast Ranges) at lower elevations. This mass loss acceleration is linked with the change from a negative to a positive Pacific Decal Oscillation (PDO) in 2014 which resulted in increased temperatures. In general, our time series not only display the sensitivity of glaciers to 380 climatic conditions and changes but also illustrate the limitations of linear models when deriving thickness changes, highlighting the importance of higher temporal resolution to generate robust long-term trends. This is the first study to demonstrate the ability of interferometric radar altimetry to monitor large-scale change in thickness, mass and sea-level contribution of glaciers across regions of extreme topography. This, along with recent work in the Arctic and Patagonia (Foresta et al., 2016, 2018; Gourmelen et al., 2018; Gray et al., 2015) demonstrates the potential of such a 385 system to monitor trends in ice mass on a global scale and with increased temporal resolution.


Introduction
Mountain glaciers store less than 1% of the global ice mass (Farinotti et al., 2019) however their rapid rate of mass loss accounts for almost a third of the global sea level rise, the largest sea level rise (SLR) contribution from land-ice Shepherd et al., 2020;Zemp et al., 2019). Besides representing an icon for climate change (Bojinski et al., 2014) and 25 impacting global sea level rise, the retreat and thinning of mountain glaciers also affects local communities (Immerzeel et al., 2010). Glacier retreat introduces substantial changes in seasonal water availability, which can have major societal impacts downstream, such as endangering water security for populations relying on surface water, or introducing extreme flooding (Guido et al., 2016;Ragettli et al., 2016). https://doi.org/10.5194/tc-2020-176 Preprint. Discussion started: 27 July 2020 c Author(s) 2020. CC BY 4.0 License.
Despite substantial advances with geodetic remote sensing methods, enhancing the spatial resolution and coverage of ice loss estimates (e.g. Brun et al., 2017), there is currently no demonstrated operational system that can routinely and consistently 40 monitor glaciers worldwide, especially in rugged mountainous terrain and with the necessary temporal resolution. Seasonal changes in surface elevation change are key to improving our understanding of the climatic processes that drive surface mass balance in mountain glaciers. Moreover, seasonal melt cycles of mountain glaciers have downstream impacts on livelihoods exposed to glacial risk, such as geohazards (glacier lake outbursts, glacier lake expansion, flooding) and water and food security (Huss and Hock, 2018;Pritchard, 2019). 45 Prior to CryoSat-2, radar altimetry has traditionally been limited to regions of moderate topography such as ice sheets. The launch of a dedicated radar altimetry ice mission, CryoSat-2, improvement in the ability to accurately map the ground position of the radar echoes, and the full use of the returned waveform via swath processing Gray et al., 2013;Hawley et al., 2009), has seen a near-global expansion of its application to monitoring ice mass changes beyond the two large ice sheets (Foresta et al., 2016;Gourmelen et al., 2018;Gray et al., 2015;McMillan et al., 2014b). Over regions of more 50 extreme surface topography however, such as those found in mountain glacier areas, the use of radar altimetry has been prohibited by the large pulse-limited footprint. While CryoSat's sharper footprint has led to a few promising studies over mountain glaciers (Dehecq et al., 2013;Foresta et al., 2018;Trantow and Herzfeld, 2016), there remain two limiting factors for all radar altimeters for these elevation ranges, the limited range window (240 m for CryoSat) and the closed-loop onboard tracking used to position the altimeter's range window; these 2 factors lead to a decrease in usable data in these environments 55 (Dehecq et al., 2013).
The emphasis in this study is to demonstrate the ability of interferometric radar altimetry to monitor regional mass changes in challenging rugged terrain, despite the abovementioned limitations. For this demonstration, we chose High Mountain Asia (HMA) and the Gulf of Alaska (GoA), two regions with complex terrain which have not been previously monitored with radar altimetry. We use CryoSat-2 swath altimetry to derive elevation and mass changes of mountain glaciers from 2010 to 2019, in 60 addition, we exploit the repeat cycle of CryoSat-2 to generate time series (30 days steps) at sub-regional level, giving new insights into seasonal and interannual changes within the two regions. With this study, we ultimately aim to demonstrate the https://doi.org/10.5194/tc-2020-176 Preprint. Discussion started: 27 July 2020 c Author(s) 2020. CC BY 4.0 License. potential of interferometric radar systems to contribute an independent observation of ice trends on a global scale and at high temporal resolution.

65
The HMA study area includes the Himalayan Mountain Range, the Tibetan Plateau and the Tien Shan mountains (regions 13, 14 and 15 of the Randolph Glacier Inventory) and is covered by about 100,000 km 2 of glacier area for a total of 95,540 glaciers (RGI Consortium, 2017). Climatic conditions in HMA are characterised by two main atmospheric circulation systems which impact the distribution of glaciers and glaciological changes; the Westerlies and the Indian monsoon ( Figure 1). The Westerlies dominate regions in the northwest (Pamir regions, Kunlun Shan, Tien Shan and western Himalayan mountain range) and are 70 responsible for a large fraction of the precipitation deposited particularly during the winter months (Bolch et al., 2012;Li et al., 2015;Yao et al., 2012). The Indian summer monsoon mainly influences glaciers in southern sub-regions (central and eastern Himalayan mountains, Karakoram, Nyainqêntanglha mountains), with decreasing precipitation northward (Bolch et al., 2012;Yao et al., 2012). In contrast to the monsoonal and westerly regimes, the inner Tibetan Plateau is mainly dominated by dry continental climatic conditions. Various studies have found precipitation increases in the Pamir regions and decreases 75 in the central/eastern Himalayan range, affected by changes in the two atmospheric systems, namely the strengthened Westerlies and the weakening Indian Monsoon (Treichler et al., 2019;Yao et al., 2012). As a result of atmospheric forcing, the vast majority of glaciers in the HMA region have been losing mass since the beginning of the satellite records (Bolch et al., 2019;Maurer et al., 2019) which has led to widespread glacier slowdown (Dehecq et al., 2019).

80
The GoA region, which we define to encompass the mountain range stretching along the Gulf of Alaska to British Columbia (Region 1 of the Randolph Glacier Inventory 6.0, excluding Northern Alaska) is covered by approximately 86,000 km 2 glacier area for a total of 26,490 glaciers (RGI Consortium, 2017). The glacierised areas stretch from sea level up to over 5000 m. a.s.l., representing a large variety of different glacier types. 67% of the glacier area are in land-terminating glaciers, 13% and 20% are marine-terminating and lake-terminating respectively ( Figure 1). Large glacier-to-glacier variations in mass changes 85 have been reported, which are assumed to be driven by climate variability and heterogeneity of glacier elevation ranges (Larsen et al., 2015). The coastal regions along the Alaskan Gulf experience a maritime climate, with the maximum precipitation occurring on the southern slopes of the Coastal Range (Wendler et al., 2017). These mountain ranges act as barriers for the moist air from the Pacific Ocean resulting in rain shadow, i.e. more continental climate, on their leeward side (Le Bris et al., 2011;Wendler et al., 2017). The Pacific Decadal Oscillation (PDO) is another factor which exercises substantial influence on 90 the climate (Wendler et al., 2017) and glacier behaviour (Hodgkins, 2009) within the GoA region. In general, the positive phase of the PDO relates to higher temperatures and more precipitation (Fleming and Whitfield, 2010), whilst a cooling and decrease of precipitation are observed during its negative phase (Papineau, 2001). However, the effects on precipitation especially are spatially heterogenous (Fleming and Whitfield, 2010). Our study period of 2010 to 2019 contains the 2014 change from a negative phase of the PDO to a positive phase, contributing to a substantial increase in temperatures in Alaska 95 from 2014 onwards (Wendler et al., 2017). As a result of atmospheric and oceanic forcings, glaciers in the GoA region have https://doi.org/10.5194/tc-2020-176 Preprint. Discussion started: 27 July 2020 c Author(s) 2020. CC BY 4.0 License. been losing mass since the beginning of the satellite records (Arendt et al., 2002;Berthier et al., 2010;Wouters et al., 2019;Zemp et al., 2019).

Data and Methods
In this section, we give a short overview on the data and methods used in this study. More details can be found in the Supporting 100 Information.

Time-dependent elevation from CryoSat-2 observations
We use observations from the SAR Interferometric Radar Altimeter (SIRAL) onboard the European Space Agency (ESA) CryoSat-2 satellite (Wingham et al., 2006). SIRAL is a beam-forming active microwave radar altimeter with a maximum imaging range of ~15 km on the ground. The sensor emits time-limited Ku-band pulses aimed at reducing the footprint to ~1.6 105 km within the beam. Over land-ice, the sensor operates in synthetic aperture interferometric (SARIn) mode, which allows delay-Doppler processing to increase the along-track resolution to ~380 m, while cross-track interferometry is used to extract key information about the position of the footprint centre. CryoSat-2 orbits the Earth with a 369-day near-repeat period formed by the successive shift of a 30-day sub-cycle. The satellite has an inclination of 92, offering improved coverage of the polar regions. We process level 1b, baseline C data and the corrected mispointing angle for aberration of light (Scagliola et al., 2018) 110 supplied by the ESA ground segment using a swath processing algorithm . Level 1b data is provided as a sequence of radar echoes along the satellite track, which translates into single power, interferometric phase and coherence waveforms for each along-track location. The conventional level 1b data processing method consists in extracting single elevation measurements from the power signal in each waveform that corresponds to the Point of Closest Approach (POCA) between satellite and the ground. In contrast, swath altimetry exploits the full radar waveform to map a dense swath (~5 km 115 wide) of elevation measurements across the satellite ground track beyond POCA (Foresta et al., 2016(Foresta et al., , 2018Gourmelen et al., 2018;Gray et al., 2013;Hawley et al., 2009) providing one to two orders of magnitude more elevation measurements compared with POCA. This makes the CryoSat-2 sensor at present the only radar altimeter able to survey small glaciers and ice caps at high resolution.

2.2
Rates of elevation change maps 120 Rates of elevation change and mass balance are based on ~25 million swath elevation measurements in the GoA region and ~8 million swath elevation measurements in HMA from 2010 to 2019. The altitudinal distribution of the measurements departs somewhat from the glaciers' hypsometry. Given the general correlation between ice thickness change and hypsometry in the region of imbalance, a bias in the altitudinal distribution of observations can lead to a bias in the total rate of thinning when integrated over a larger domain. Therefore we derive a subset of the time-dependent elevation dataset, removing the impact of 125 such point density biases by filtering out swath measurements so as to match the glacier hypsometry (e.g. Treichler et al., https://doi.org/10.5194/tc-2020-176 Preprint. Discussion started: 27 July 2020 c Author(s) 2020. CC BY 4.0 License. 2019), and generate elevation change and mass change estimates from the reduced sample ( Figure S4). We remove data sequentially based on measurement uncertainty. This process reduces our sample size by 15% for the GoA and by 30% for HMA. To generate rates of elevation change we then follow a similar approach to Foresta et al. (2016) and Gourmelen et al. (2018), however the lower data density and the complexity of the terrain in the GoA region and in HMA require a slight 130 adaptation of the methodology. We bin the elevation measurements into regions of 100 x 100 km, sufficiently large to contain the necessary number of measurements in each bin to ensure sufficient robustness and representativity. Due to the increased bin size (the pixel size used by Foresta et al. [2016] and Gourmelen et al. [2018] is 1000 m) and the variation of elevations within each bin, the topographic signature cannot simply be modelled and therefore needs to be removed using an auxiliary Digital Elevation Model (DEM) . We subtract the TanDEM-X 90m DEM (German Aerospace Center [DLR], 135 2018), which has a near-complete coverage and is contemporaneous of the CryoSat-2 observations, from the swath elevation measurements. The remaining elevation differences (hereinafter referred to as elevDiff) are due to time-dependent elevation change that can be related to glacier thickness change as well as errors in the two data sets, temporal heterogeneity and differences in penetration between the reference DEM and the swath elevation measurements.
Time-dependent rates of surface elevation change are calculated for each 100 x 100 km bin individually based on the elevDiff 140 measurements. In order to achieve the most robust trends we considered several fitting methods, including ordinary leastsquare, robust regression (e.g. Kääb et al., 2012Kääb et al., , 2015, weighted regression (e.g. Berthier et al., 2016;Foresta et al., 2018;Gourmelen et al., 2018), random sample consensus (RANSAC) and Theil-Sen estimator (e.g. Shean et al., 2020). We found that the best results were achieved with a weighted regression model of the elevDiff measurements, similar to the methods of Gourmelen et al. (2018). However, whilst their weights are calculated only according to the power attribute, here we assign 145 each observation a weight based on power and coherence, i.e. measurements with high power and low coherence within the sample will have lower weights assigned (see Supplementary Information S1.1). When fitting the model, we iteratively exclude measurements that are more than 3σ from the mean distance to the fitted line, until no more outliers are present (e.g. Foresta et al., 2016Foresta et al., , 2018. The rate of change uncertainty is extracted from the standard error of the regression model. We conservatively use a factor of five for uncertainties on areas without coverage of swath measurements (see Supplementary 150 Information S1.4). After generating time-dependent elevation changes and elevation change uncertainty for each 100 x 100 km bin we discarded bins that did not fulfil a set of quality criteria (see Supplementary Information S1.2). The remaining bins covered more than 96% of the total glacierised area in the GoA region, and ~90% in HMA. To estimate values for the gaps in our dh/dt map we use the altitudinal distribution of elevation change rates on a sub-regional level (Moholdt et al., 2010a(Moholdt et al., , 2010bNilsson et al., 2015), applying the hypsometric averaging methods described in the Supplementary Information (S1.3). 155 It is a well known observation that microwave pulses scatter from the surface as well as the subsurface, which can lead to elevation change bias in regions of historically anomalous melt event (Nilsson et al., 2015); over most regions however, it has been shown that surface elevation change from CryoSat is consistent with in-situ, airborne, and meteorological observations Gray et al., 2015;McMillan et al., 2014a;Zheng et al., 2018) https://doi.org/10.5194/tc-2020-176 Preprint. Discussion started: 27 July 2020 c Author(s) 2020. CC BY 4.0 License.

2.3
Mass balance and contribution to sea level rise 160 To obtain volume changes we use the glacierised area of the Randolph Glacier Inventory (RGI 6.0) (RGI Consortium, 2017).
We assume the standard bulk density of 850 kg/m 3 (Huss, 2013) to convert volume changes to equivalent mass changes. Mass change uncertainty is calculated based on the main three error sources: rates of change uncertainty, uncertainty of the glacierised area and volume-mass conversion uncertainty (see Supplementary Information S1.4). To obtain a region-wide mass balance, mass changes of each individual bin are summed up. We derive mass balance for the unadjusted (biased) and glacier 165 hypsometry adjusted (unbiased) elevation dataset (see Section 2.2), our final mass balance numbers are for the unbiased case.
To generate the contribution to sea level rise (SLR) we exclude large endorheic basins in High Mountain Asia and assume an area of the ocean of 361.8 · 10 6 km 2 .

Time series of seasonal surface elevation changes
CryoSat-2's monthly repeat cycle provides the opportunity to monitor seasonal as well and multiannual trends of surface 170 elevation. We therefore generate time series with a monthly step (30 days) and a 3-month (90 days) moving window, using the median of all the elevDiff observations (residuals from the reference DEM) within a time period with reference to the first month. Time series are generated on bin size level (100 x 100 km), on sub-regional level (using the RGI sub-regions) and for the full study region. The time series of the 100 x 100 km bins are mainly used as an additional check of the dh/dt quality (see Supplementary Information S1.2), whilst we exploit the sub-regional and regional time series to analyse spatio-temporal 175 variability in thickness change across both, the GoA and HMA regions. To generate region-wide time series for HMA and the GoA we use an area-weighted mean of the sub-regional time series. Note that the regional and sub-regional time series displayed in this publication start in January 2011 (with the earliest data from November 2010 using the 90 days window), since we retrieve less swath measurements for the first few months of CryoSat-2's life cycle, impacting the quality of the time series pre-2011. 180 3 Results

Spatial coverage and elevation sampling
Using the pulse-limited footprint size of CryoSat-2, we achieve a total spatial coverage of glaciated regions of 55% in the GoA and 32% in HMA respectively. These values are the combined result of the absence of recorded returns due to orbit separation and onboard-tracking limitation (Dehecq et al., 2013), and of data quality. Given that 40% of HMA glaciers are not measured 185 due to onboard tracking limitations (Dehecq et al., 2013) the rate coverage of the measurable area is about 50%. These values are within the high-end of the range of observational methods (Zemp et al., 2019), whist generally lower than the coverage provided by very-high resolution sensors Shean et al., 2020). Despite the relatively large size footprint of https://doi.org/10.5194/tc-2020-176 Preprint. Discussion started: 27 July 2020 c Author(s) 2020. CC BY 4.0 License. radar altimeters, we observe a coverage of all glacier sizes, with approximately 10% area of very small glaciers covered ( Figure   S5), as expected we do observe a positive correlation between spatial coverage and glacier size. 190 Representative elevation sampling of swath measurements is important, since change in thickness is strongly correlated with elevation and any systematic elevation bias of the sampling has the potential to affect dh/dt trends. Spatial coverage and number of points show a different relationship with hypsometry, which is due to the overlap between adjacent CryoSat footprints.
Comparing the glacier hypsometry with the spatial coverage of our data shows no significant region-wide elevation bias, however, we observe a bias of the total number of swath measurements towards higher altitudes (e.g. Figure S4), which can 195 be attributed to the onboard tracking tending to favour elevations closest to the satellite.
As expected, we retrieve slightly smaller specific mass losses with differences in the order of 0.05 m w.e. yr -1 for the GoA region and 0.02 m w.e. yr -1 for HMA without removing the elevation bias towards higher altitudes. We do not observe any significant changes in the overall spatial pattern of elevation change across the two regions between the unbiased and biased results, which indicates that we still achieve a stable result with the reduced sample. 200

Elevation changes and mass balance in High Mountain Asia
The total unbiased HMA mass balance between 2010 and 2019 was -27.9 ± 2.4 Gt yr -1 (-0.29 ± 0.03 m w.e. yr -1 ) and -18.3 ± 1.6 Gt yr -1 (-0.38 +-0.03 m w.e. yr -1 ) excluding endorheic basins. Our maps of surface elevation change show a heterogenous pattern in the Himalayan range, with a cluster of slightly positive/near balance trends in the Kunlun and Karakoram ranges ( Figure 2). This spatial pattern confirms the suggestion of previous studies Gardner et al., 2013;Kääb et al., 205 2015), that the so-called "Karakoram anomaly" Yao et al., 2012) stretches up to West Kunlun Shan, which has now become the centre of the anomaly. Another striking feature is the gradient from moderate thinning in Spiti-Lahaul and western Himalaya to increasingly negative surface elevation changes along the central and eastern Himalayan mountain range, with the Nyainqêntanglha mountains and Hengduan Shan showing the highest negative trends. In contrast with other studies Shean et al., 2020) we find a heterogenous pattern in the Tibetan Plateau and Eastern 210 Kunlun, with some scattered glaciers displaying higher mass losses.

Glacier elevation changes and mass balance in the Gulf of Alaska
In general we find much higher mass losses in the Gulf of Alaska than in High Mountain Asia, in agreement with global studies Zemp et al., 2019). Over an area of ~86,000 km 2 , including all 26,490 glaciers in the RGI region 1 except Northern Alaska, we estimate a total mass balance of -76.3 ± 5.6 Gt yr -1 (-0.89 ± 0.07 m w.e. yr -1 ), contributing -0.199 ± 215 0.015 mm yr -1 to global sea level rise. Surface elevation change maps (Figure 3) display an expected pattern with clearly visible higher negative trends towards lower elevations close to the coast. The largest mass loss is seen in the Coast Mountains and Saint Elias Mountains, especially the Yukutat and Glacier Bay region, which is in line with the spatial patterns of Luthcke et al. (2008) and Luthcke et al. (2013). The lowest thinning rates are observed in the Alaska Range mountains, which is also in agreement with other studies (Berthier et al., 2010;Luthcke et al., 2008).

Contrarily, sub-regions influenced by the summer Indian monsoon (Central Himalaya, Eastern Himalaya and Hengduan Shan),
where glacier accumulation and ablation happen at the same time, show a more heterogenous seasonal pattern. In the time 235 series of these three subregions the annual cycle has two peaks, with the first peak reflecting thickening in winter and the second and smaller peak reflecting the summer accumulation through the influence of the Indian summer monsoon (see Figure   S1). The inner Tibetan Plateau, dominated by a more continental climate, displays almost no intra-annual cycle. In general, the heterogeneity of the time series reflects the sensitivity of mountain glaciers to meteorological patterns and changes.

Altitudinal distribution 240
We display the altitudinal distribution of elevation changes in Figure 6 and a comparison with Brun et al. (2017) in Figure S6.
In general, we observe decreasing negative trends with increasing altitudes, which is an expected pattern Gardelle et al., 2013). We find the steepest gradient ( Figure 6, S6, Table S3) in the Nyainqêntanglha/Hengduan Shan, which is in line with the findings of Brun et al. (2017). We also observe lower or even inverse gradients in Bhutan/East Himalaya, Spiti-Lahaul/West Himalaya, Karakoram/West-Kunlun and Pamir (Figure 6, S6, Table S3), which have been reported 245 previously and been related to debris thickness (Bisset et al., 2020;Brun et al., 2017). The main differences are observed in the elevation profile of the Kunlun regions, where Brun et al. (2017) find constant thickening at all elevation during the survey time period of 2000 to 2016, whilst we record thinning at lower elevations (see Figure 6). These findings suggest a shift towards negative mass balance for the Kunlun region in comparison to the previous decade.

Comparison of regional mass balance with previous work 250
A comparison of mass balance results in the literature indicates that, while all the studies agree on the general trend in mass loss and spatial variability of mass loss, there is large degree of variability between estimates. While some of the variability can be attributed to the diversity of time-span and region boundaries used, there are also clear differences between observation methods (Figure 8a). Note that here we are only comparing region-wide mass trends with the results closest in space and time to this study, whilst sub-regional differences are discussed in the next section. 255 Our total mass balance of -27.9 ± 2.4 Gt yr -1 (-0.29 ± 0.03 m w.e. yr -1 ) is in good agreement with the -28.8 ± 12 Gt yr -1 by This trend seems to have continued in more recent years, with Ciracì et al. (2020) observing an acceleration in mass loss of 10 265 ± 5 Gt yr -1 per decade for the period of 2002 to 2019, which could explain our more negative mass balance in comparison to Brun et al. (2017Brun et al. ( ) [2000Brun et al. ( to 2016 and Shean et al. (2020Shean et al. ( ) [2000Shean et al. ( to 2018.

4.1.4
Comparison of sub-regional mass balances with previous work Our higher regional mass loss when comparing to the two DEM differencing studies by Brun et al. (2017) and Shean et al. (2020) are mostly down to differences in the South-eastern Himalayaespecially Nyainqêntanglha/Hengduan Shanand in 270 the Pamir regions. We used the regions by Brun et al. (2017) and the RGI 6.0 second order regions to compare our results with other estimates (Figure 7, Figure S7, Table S1). For a full discussion of regional differences between estimates of recent studies refer to Bolch et al. (2019). Our results are in line with general findings by Bolch et al., (2019)  periods 2000-2016 and 2000-2018] fit in with the generally observed acceleration of mass loss in South-East Asia over the past decades (Maurer et al., 2019;Zemp et al., 2019). Some studies suggest the weakening of the Indian summer monsoon as the primary source of increased thinning (Salerno et al., 2015), whilst other studies find no widespread precipitation decrease in monsoonal regimes which could account for all of these changes and attribute the temperature sensitivity of glaciers in 285 monsoon-dominated regions as the main driver (Maurer et al., 2019). In fact, glaciers in Hengduan Shan, Nyainqêntanglha and Eastern Himalaya have been found to exhibit the highest sensitivity towards temperature in the whole HMA region (Sakai and Fujita, 2017).
Contradictory estimates have also been published for the Pamir and Pamir Alay mountains (Hissar Alay), where high , moderate (this study), slight mass losses Shean et al., 2020), and even mass gains (Gardelle et 290 al., 2013) have been reported. Time series by Ciracì et al. (2020) suggest an acceleration of mass losses within the Pamir-Karakoram region since 2011, which could explain our higher mass loss estimates in comparison to the DEM differencing studies Gardelle et al., 2013;Shean et al., 2020). In fact, our time series suggest that the glaciers in Pamir and Hissar Alay were indeed thinning moderately from 2011, with a slowdown from 2016-7 onward.
We record less mass gain in Kunlun (+0.01 ± 0.03 m w.e. yr -1 ; +0.06 ± 0.03 m w.e. yr -1 in the Western part of Kunlun) than 295 previous studiesa potential sign for the widely discussed and predicted disappearance of the Karakoram anomaly. However, our time series suggests increased mass gain from 2016 in Western Kunlun, and also mass gain in the Karakoram (Figure 5).
At the same time, we also observe decreased thinning rates in Inner Tibet and East Kunlun. These changes could be a shortterm trend, however, it displays the limitation of all mentioned studies (including this study) when deriving linear trends in a region like High Mountain Asia with large inter-annual climate variability and associated glacier changes. 300 We generally find better agreements with Shean et al. (2020) the study including an additional 2+ years (2017 and 2018) in comparison to Brun et al. (2017), and thus more closely aligned with our time period. This potentially indicates that a large part of the disagreements could be related to inter-annual variability and survey time period.

4.2
Gulf of Alaska 305

Temporal variability
We display temporal variability of surface elevation change for the whole GoA region (Figure 4), the RGI sub-regions ( Figure   9, S3) and for different elevation bands within sub-regions ( Figure 10). Figure 9 shows negative trends across all the subregions. The four coastal sub-regions -Alaska Pena, Western Chugach Mountains, Saint Elias Mountains and Coast Ranges display a seasonal oscillation, with an annual surface elevation maximum in spring and annual surface elevation minimum 310 in autumn. In contrary, the seasonal cycle of Alaska Range mountains is shifted, with the thickness maximum in winter, which is also somewhat visible in the time series by Luthcke et al. (2008). This change correlates with the change from a negative PDO phase to a positive phase in 2014, which resulted in increased temperatures (Wendler et al., 2017). Our findings suggest that the sensitivity of glaciers to the 2014 temperature change increases towards lower elevations (Figure 10). The fact that the strongest impact is observed in the coastal regions is likely due the higher sensitivity of maritime glaciers in Alaska to temperature change (Gregory and Oerlemans, 1998) and the lower 320 elevations within these regions.

Altitudinal distribution
We observe a clear correlation between surface elevation changes and altitude ( Figure 11, Table S2), with the highest negative trends at low altitudes in the Saint Elias Mountains and Coast Ranges. In comparison to the study by Berthier et al. (2010), which is based on sequential Digital Elevation Models over the time period from 1962 to 2006, we observe similar elevation 325 changes at the lowest altitudes but less steep gradients in sub-regions along the east coast. This is particularly pronounced in the Saint Elias Mountains, where Berthier et al. (2010) show near-balance at around 1000 metre above sea level, whilst our estimates suggest a surface elevation change of -1.5 m yr -1 at the same altitude. These findings indicate a propagation of thinning upstream since 1962 to 2006. In contrary, whilst on the Alaska Peninsula elevation changes have increased at lower altitudes, the limit of the thinning area has stayed the same since the survey time period of Berthier et al. (2010). 330 In the Western Chugach Mountains, Alaska Range and Alaska Peninsula we observe the effect of debris cover on the altitudinal profile, with a decrease of thinning rates towards the lowest elevations of these sub-regions. This characteristic has been observed, although more pronounced and across all sub-regions, by Berthier et al. (2010) and Arendt et al. (2002).
negative than most other studies' findings, reflecting the increased thinning rates we show in the sub-regional time series from

345
Since there is no prevalent sub-region mask used by more recent studies, we cannot directly compare and validate our results on a sub-regional level. Mass balance or surface elevation change estimates that overlap with our time period are either spatially not resolved (e.g. Gardner et al., 2013;Zemp et al., 2019), presented on a glacier to glacier basis (e.g. Larsen et al., 2015) or GRACE mascon extents (e.g. Luthcke et al., 2008Luthcke et al., , 2013. In addition to the comparison with Berthier et al. (2010) we present sub-regional estimates aggregated on the RGI 6.0 second order regions in Table 2 for future reference.

Conclusion
We exploit CryoSat-2 interferometric swath processed data from 2010 to 2019, with a total of 33 million elevation observations, to generate new and independent mass balance estimates for two mountain regions, the Gulf of Alaska (GoA) 360 and High Mountain Asia (HMA). We also generate observations at sub-regional level and extract elevation-dependant thinning rates, revealing contrasting mass loss across sub-regions. Finally, we extract monthly time series of elevation change, exploiting CryoSat's high temporal repeat, to reveal seasonal and multiannual variation in rates of glaciers' thinning. We find that between 2010 and 2019, HMA has lost mass at rates of 27.9 ± 2.4 Gt yr -1 (0.29 ± 0.03 m w.e. yr -1 ), and the GoA region has lost mass at rates of 76.3 ± 5.6 Gt yr -1 (0.89 ± 0.07 m w.e. yr -1 ), for a sea-level contribution of 0.048 ± 0.004 mm yr -1 and 365 0.217 ± 0.015 mm yr -1 respectively for HMA and the GoA. These estimates are broadly consistent with the range of estimates generated by previous studies and highlight the significant discrepancies that remain in the assessments of mass loss for these two regions.
In HMA we find the most negative surface elevation trends in the Nyainqêntanglha mountains, Hengduan Shan, the East Himalayan range and the Tien Shan and slightly positive/near balance trends in the Kunlun and Karakoram ranges, known as 370 the "Karakoram anomaly". This spatial heterogeneity is known to be influenced by the two main prevailing atmospheric circulation systems in the Himalayas; the Indian summer monsoon and the winter Westerlies. Our monthly time series give insights into the impact of these two systems on glacier dynamics, with the sub-regional accumulation and ablation cycles being influenced by the seasonality of the precipitation. We show sustained multiannual trends across almost all of the subregions until 2015-6, and decreased loss or even mass gain from 2016/2017 onwards. in 2014 which resulted in increased temperatures. In general, our time series not only display the sensitivity of glaciers to 380 climatic conditions and changes but also illustrate the limitations of linear models when deriving thickness changes, highlighting the importance of higher temporal resolution to generate robust long-term trends. This is the first study to demonstrate the ability of interferometric radar altimetry to monitor large-scale change in thickness, mass and sea-level contribution of glaciers across regions of extreme topography. This, along with recent work in the Arctic and Patagonia (Foresta et al., 2016(Foresta et al., , 2018Gourmelen et al., 2018;Gray et al., 2015) demonstrates the potential of such a 385 system to monitor trends in ice mass on a global scale and with increased temporal resolution.