Spatially and temporally resolved ice loss in High Mountain Asia and the Gulf of Alaska observed by CryoSat-2 swath altimetry between 2010 and 2019

Abstract. Glaciers are currently the largest contributor to sea
level rise after ocean thermal expansion, contributing ∼ 30 % to the sea level budget. Global monitoring of these regions remains a
challenging task since global estimates rely on a variety of observations
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 exploit CryoSat's monthly temporal repeat to
reveal seasonal and multiannual variation in rates of glaciers' thinning at
unprecedented spatial detail. We find that during this period, HMA and GoA
have lost an average of −28.0 ± 3.0 Gt yr−1 (−0.29 ± 0.03 m w.e. yr−1) and −76.3 ± 5.7 Gt yr−1 (−0.89 ± 0.07 m w.e. yr−1), respectively, corresponding to a contribution to sea level rise of 0.078 ± 0.008 mm yr−1 (0.051 ± 0.006 mm yr−1 from exorheic basins) and 0.211 ± 0.016 mm yr−1. The cumulative loss during the 9-year period is equivalent to 4.2 % and 4.3 % of the ice volume, respectively, for HMA and GoA. Glacier thinning is ubiquitous except for in the Karakoram–Kunlun region, which experiences stable or slightly
positive mass balance. In the GoA region, the intensity of thinning varies
spatially and temporally, with acceleration of mass loss from −0.06 ± 0.33 to −1.1 ± 0.06 m yr−1 from 2013, which
correlates with the strength of the Pacific Decadal Oscillation. In HMA ice
loss is sustained until 2015–2016, with a slight decrease in mass loss from
2016, with some evidence of mass gain locally from 2016–2017 onwards.



Introduction
Glaciers store less than 1% of the mass (Farinotti et al., 2019) and occupy just over 4% of the area (RGI Consortium, 2017) of global land ice, however their rapid rate of mass loss accounts for almost a third of the global sea level rise,has accounted for almost a third of the global sea level rise during the 21 st century (Gardner et al., 2013;WCRP Global Sea Level Budget Group, (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 65 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 potential of interferometric radar systems to contribute an independent observation of ice trends on a global scale and at high temporal resolution.

Study regions 70
The HMA study area includes the Himalaya, the Tibetan mountain ranges, the Pamir and Tien Shan (regions 13, 14 and 15 of the Randolph Glacier Inventory) and is covered by about 100,000 km 2 of glacier area for about 95,500 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 75 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 80 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 during the satellite records (Bolch et al., 2019;Farinotti et al., 2015;Maurer et al., 2019) which has led to widespread glacier slowdown (Dehecq et al., 2019).

85
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 about 26,500 glaciers (RGI Consortium, 2017). The glacierised areas stretch from sea level up to over 50006000 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 90 in mass changes 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 95 substantial influence on the climate (Wendler et al., 2017) and glacier behaviour (Hodgkins, 2009) within the GoA region. In Field Code Changed 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 20142013-4 change from a negative phase of the PDO to a positive phase, contributing to a substantial increase 100 in temperatures in Alaska from 2014 onwards (Wendler et al., 2017). As a result of atmospheric and oceanic forcings, glaciers in the GoA region have been losing mass during 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 105 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). (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 110 the footprint to ~1.6 km within the beam. Over land-ice, the sensor operates in synthetic aperture interferometric (SARIn) mode, which allows delay-Doppler processing to generate an along-track footprint of ~380 m, while cross-track interferometry is used to extract key information about the position of the footprint centre. In practice however, footprint size will vary depending on properties such as surface slopes, scattering properties, and distance from the Point-of-Closest-Approach (POCA). CryoSat-2 orbits the Earth with a 369-day near-repeat period formed by the successive shift of a 30-day sub-cycle. 115 The satellite has an inclination of 92, offering improved coverage of the polar regions. We process level 1b, baseline C data and the corrected mis-pointing angle for aberration of light (Scagliola et al., 2018) 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 received 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 120 power signal in each waveform that corresponds to the POCA between satellite and the ground. In contrast, swath altimetry exploits the full radar waveform to map a dense swath (~5 km wide) of elevation measurements across the satellite ground track beyond POCA (Foresta et al., 2016Gourmelen et al., 2018;Gray et al., 2013;Hawley et al., 2009) (Foresta et al., 2016Gourmelen et al., 2018;Gray et al., 2013;Hawley et al., 2009) providing one to two orders of magnitude more elevation measurements compared with POCA and improving the sampling of topographic lows (Foresta at al., 2016). Because 125 swath processing does not rely on retracking, it can retrieve elevation measurements also for atypical waveforms with no clearly defined leading edge such as those found over complex terrain and where retracking often fails to identify a reliable POCA. This makes the CryoSat-2 sensor at present the only radar altimeter able to survey small glaciers at high resolution.

Rates of elevation change maps
Rates of elevation change and mass balance are based on ~25 million swath elevation measurements in the GoA region and 130 ~8 million swath elevation measurements in HMA acquired from mid-2010 to mid-2019. The distribution of elevation measurements with altitude departs somewhat from the glaciers' hypsometry. Hypsometric representativeness of samples within spatial units is a key requirement for robust glacier trend estimates. 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 as rate of thickness change is often strongly correlated with altitude. Therefore we derive a subset of the time-dependent elevation dataset, removing the impact of such 135 point density biases by filtering out swath measurements so as to match the glacier hypsometry binned using 100m elevation intervals (e.g. Treichler et al., 2019), and generate elevation change and mass change estimates from the reduced sample ( Figure S5). 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 weWe then follow a similar approach to Foresta et al. (2016) and Gourmelen et al. (2018)Foresta et al. (2016 and Gourmelen et al. (2018), however the lower data density and 140 the complexity of the terrain in the GoA region and in HMA require a slight 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]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 145 Digital Elevation Model (DEM) . We subtract the TanDEM-X 90m DEM (German Aerospace Center [DLR], 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 (TanDEM-X is a composite of acquisitions from different years) and differences in penetration between the reference DEM (X-band) and 150 the swath elevation measurements. The errors related to the reference DEM will result in an increase in spread of the elevDiff measurements and is accounted for in the regression model discussed below.
Time-dependent ratesRates of surface elevation change are then calculated for each 100 x 100 km bin individually based on the elevDiff measurements from mid-2010 to mid-2019. In order to achieve the most robust trends we considered several fitting methods, including ordinary least-square, robust regression (e.g. Kääb et al., 2012Kääb et al., , 2015)(e.g. Kääb et al., 2012Kääb et al., , 2015, 155 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 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 160 Supplementary Information S1.1). We exclude solutions that display extremely large variability across various regression models, considering them as unstable solutions results (see Supplementary Information S1.2). 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. We discard bins that did not fulfil a set of quality criteria based on elevation change uncertainties, temporal completeness, interannual changes and stability of regression results (see Supplementary Information 165 S1.2). The remaining bins covered more than 96% of the total glacierised area in the GoA region, and 88% 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).
The rate of change uncertainty is extracted from the standard error of the regression model. We conservatively use a factor of 170 five for uncertainties on areas without coverage of swath measurements (see Supplementary Information S1.4). While our uncertainty methods follow existing approaches and our error bounds are similar in magnitude to Brun et al. (2017),  and Shean et al. (2020) but lower than GRACE-based estimates, several additional potential sources of errors could impact the results, and methods to assess them, not currently available, should be developed. Radar altimetry has been shown to be sensitive to surface slopes, and in particular to slope in the direction of the satellite's flight path, in regions like HMA 175 and GoA this impact will also be seen in the performance of the onboard tracker as for large slopes the system is expected to "lose lock". 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 events (Nilsson et al., 2015); or at seasonal time-scale (Gray et al., 2019). Over most regions however, it has been shown that surface elevation change from CryoSat over annual and pluri-annual time scale are consistent with in-situ, airborne, and meteorological observations 180 Gray et al., 2015180 Gray et al., , 2019McMillan et al., 2014a;Zheng et al., 2018). Using static glacier masks can also lead to errors in regions of rapid dynamic changes. In general, these limitations are known and efforts are currently underway in the community to improve uncertainty analysis, and develop new glaciers outlines products.

Mass balance and contribution to sea level rise
To obtain volume changes we use the glacierised area of the Randolph Glacier Inventory (RGI 6.0) (RGI Consortium, 2017). 185 We assume the standard bulk density of 850 ± 60 kg/m 3 (Huss, 2013) to convert volume changes to equivalent mass changes.
This assumption is considered appropriate for a wide range of conditions and longer-term trends, however, this factor can differ significantly for shorter term periods (<3 years) (Huss, 2013). 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 190 up. We derive mass balance for the unadjusted (biased) and glacier 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 assume an area of the ocean of 361.8 · 10 6 km 2 and consider total contributions from all glaciers and then only those glaciers within exorheic basins in High Mountain Asia, based on the HydroSHEDS dataset (Lehner et al., 2006).

Time series of surface elevation changes 195
CryoSat-2's monthly repeat cycle provides the opportunity to monitor seasonal as well and multiannual trends of surface 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 also used as an additional check of the dh/dt quality (see 200 Supplementary Information S1.2), whilst we exploit the sub-regional and regional time series to analyse spatio-temporal 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 205 series pre-2011. Although time-Note that as opposed to the linear rates 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.
The time series in this paper end in April 2019, with the latest data from June 2019 due to the 90 days window. series are generally reflecting the actual change in surface elevation, there are a number of limitations that are important to keep in mind 210 when interpreting the results from radar altimetry. For the reasons stated above, scattering properties can induce elevation biases at seasonal time-scale (Gray et al., 2019). In addition, integrating changes over very large regions can lead to spatial heterogeneity in the successive time steps, in particular when the data volume becomes too low. These limitations may explain some of the observed patterns, and in particular the few cases where seasonal variability is larger than what is expected fro m our knowledge of SMB in the regions. 215

Uncertainty assessment
The error budget on mass change has three uncertainty sources, which are assumed to be independent and uncorrelated: uncertainty on time-dependent elevation change ( ℎ ), uncertainty on glacierised area ( ) and uncertainty on mass-volume conversion ( ).
The rate of elevation change uncertainty for each 100 x 100 km bin is based on the standard error of the regression model. We 220 conservatively use a factor of five (Berthier et al., 2014;Brun et al., 2017) for uncertainties on areas without coverage of swath measurements: where g is the proportional coverage of glacierised area at 400-metre postings, is (1 − ) and is the standard error of the regression. To retrieve the uncertainty on extrapolated bins we calculate the differences of all non-extrapolated bins between 225 elevation changes using the plane fit approach and elevation changes using the hypsometric averaging method. The standard deviation of these differences is the uncertainty on elevation change ( ℎ ) for all extrapolated bins. We retrieve an uncertainty on elevation change for extrapolated bins of 0.34 m yr -1 and 0.47 m yr -1 respectively for High Mountain Asia and the Gulf of Alaska. To account for errors due to temporal changes in glacier extents and polygon digitization  we use an error of 10% ( = 0.1 ) on the glacierised area in a bin, even though the reported uncertainty of the RGI is ~8% (Pfeffer 230 et al., 2014). Assuming independence between the two error components ( , ℎ ), volume change uncertainty ( ) of a bin is: where ℎ is the elevation change rate of the respective bin. To generate the region-wide volume uncertainty ( ) we combine all the values (including extrapolated bins) in quadrature. We use a density uncertainty of = 60 kg m -3 , and a 235 density mass conversion of p = 850 kg m -3 (Huss, 2013). The total mass balance uncertainty is: where is the total volume change for the region.

Spatial coverage and elevation sampling 240
Using the theoretical pulse-limited footprint size of CryoSat-2, we derive 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 it is estimated that 40% of HMA glaciers are not sampled due to onboard tracking limitations (Dehecq et al., 2013) we estimate that with an appropriate onboard tracking system, the rate coverage for HMA would be as high as 50%. These values are within the high-end of the 245 range of observational methods (Zemp et al., 2019), whilst generally lower than the coverage provided by high resolution sensors Shean et al., 2020). As expected from the relatively large footprint of radar altimeters, we do observe a positive correlation between spatial coverage and glacier size, we do however observe coverage over all glacier sizes ( Figure   S6).
We observe a bias of the total number of swath measurements towards higher altitudes (e.g. Figure S5), which can be attributed 250 to the onboard tracking tending to favour elevations closest to the satellite. However, a comparison of the glacier hypsometry and the spatial coverage of our data ( Figure S4) shows that we still achieve good coverage at low elevations in both regions.
In addition, we interpolate missing data based on the relationship between elevation and elevation changes and therefore still capture the changes in the lower reaches of the HMA and GoA glaciers.

Elevation changes and mass balance in High Mountain Asia 255
The total HMA mass balance between 2010 and 2019 was -28.0 ± 3.0 Gt yr -1 (-0.29 ± 0.03 m w.e. yr -1 ), or -18.3 ± 2.3 Gt yr -1 (-0.38 ± 0.03 m w.e. yr -1 ) when including only exorheic basins. This mass loss corresponds to 0.078 ± 0.008 mm yr -1 SLE, or 0.051 ± 0.006 mm yr -1 when including only exorheic basins. 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 260 change across the two regions between the unbiased and biased results, which indicates that we still achieve a stable result with the reduced sample. Our maps of surface elevation change show a heterogeneous pattern in the Himalayan range, with a cluster of slightly positive/near balance trends in the Kunlun and Karakoram ranges (Figure 2), the so-called "Karakoram anomaly" Hewitt, 2005). Another striking feature is the gradient from moderate thinning in Spiti-Lahaul and western Himalaya (-0.25 ± 0.09 m w.e. yr -1 ) to increasingly negative surface elevation changes along the central (-0.43 265 ± 0.14 m w.e. yr -1 ) and eastern (-0.56 ± 0.16 m w.e. yr -1 ) Himalayan mountain range, with the Nyainqêntanglha mountains and Hengduan Shan (-0.98 ± 0.22 m w.e. yr -1 ) showing the highest negative trends.
We display the altitudinal distribution of elevation changes in Figure 6 and a comparison with Brun et al. (2017) in Figure S7.
While some variability exists along the profiles, in particular over regions and elevation range containing fewer glaciers that can reflect a less robust solution and or spatial variability in glacier response, trends between elevation and ice thickness change 270 are clearly visible. 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, S7, Table S4) 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, S7, Table S4), which have been reported previously and been related to debris thickness (Bisset et al., 2020;Brun et al., 2017). 275 We show temporal variability of surface elevation change for the whole HMA region (Figure 4), the RGI second order regions ( Figure 5, S1) and the regions by Brun et al. (2017) (Figure S2). The monthly time series show sustained multiannual trends across almost all of the subregions until 2015-6, and decreased loss or even mass gain from 2016/2017 onwards ( Figure 5, S2), which is also reflected in the full HMA time series (Figure 4).), and consistent with previous observation (Ciracì et al., 2020).

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. 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.7 Gt 285 yr -1 (-0.89 ± 0.07 m w.e. yr -1 ), contributing -0.211 ± 0.016 mm yr -1 to global sea level rise. Surface elevation change maps ( Figure 3) display an expected pattern with clearly visible highermore negative trends towards lower elevations close to the coast. Note that some of the lower rates observed in the St Elias Mountain are likely the result of the presence of accumulation areas of large glaciers e.g. Hubbard and Bering glaciers in these particular grid cells. We present sub-regional estimates aggregated on the RGI 6.0 second order regions in Table 2. The largest mass loss is seen in the Northern Coast Ranges (-1.08 290 ± 0.09 m w.e. yr -1 ; -24.8 ± 2.1 Gt yr -1 ) and Saint Elias Mountains (-1.03 ± 0.10 m w.e. yr -1 ; -34.1 ± 3.4 Gt yr -1 ), 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 (-0.41 ± 0.05 m w.e. yr -1 ), which is also in agreement with other studies (Berthier et al., 2010;Luthcke et al., 2008). 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 295 Ranges.
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 300 in autumn. In contrast, 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). A very noticeable feature within the full GoA time series is the steepening (more negative)acceleration of inter-annual elevation trendsthinning from 2013-4 onwards ( Figure 4). We record an acceleration of thinning from -0.06 ± 0.33 m yr -1 (before 2014January 2011 to January 2013) to -1.1 ± 0.06 m yr -1 (2014)(2015)(2016)(2017)(2018).January 2013 to January 2019). We observe this almost consistently across all the five sub-regions, but this is 305 most pronounced in the Saint Elias Mountains, the Western Chugach mountains and Coast Ranges.

Uncertainty
While our uncertainty methods follow existing approaches and our error bounds are similar in magnitude to Brun et al. (2017), Kääb et al. (2012) and Shean et al. (2020) but lower than GRACE-based estimates, several additional potential sources of 310 errors could impact the results. Radar altimetry, and delay-doppler radar in particular, has been shown to be sensitive to surface slopes, and in particular to slope in the direction of the satellite's flight path, in regions like HMA and GoA this impact will also be seen in the performance of the onboard tracker as for large slopes the system is expected to "lose lock" . While we do observe a decreased coverage compared to other, less mountainous, glaciated regions, we also demonstrated here that measurements do cover the entire elevation range of glaciers in the HMA and GoA regions allowing us to match the glaciers' 315 hypsometry. We also do not observe significant coverage bias in function of glacier orientation with respect to the satellite's track path. The spatial coverage is such that we demonstrably resolve spatial, altitudinal, and temporal evolution of glacier elevation.
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 events (Nilsson et al., 2015); or at seasonal time-scale (Gray 320 et al., 2019). Over most regions however, it has been shown that surface elevation change from CryoSat over annual and pluriannual time scale are consistent with in-situ, airborne, and meteorological observations Gray et al., 2015Gray et al., , 2019McMillan et al., 2014a;Zheng et al., 2018). Using static glacier masks can also lead to errors in regions of rapid dynamic changes. In general, these limitations are known and efforts are currently underway in the community to improve uncertainty analysis, and develop new glaciers outlines products (RAGMAC, 2019).

Temporal variability
The seasonal and annual time series variability reflects the influence of atmospheric circulations and precipitation seasonality in High Mountain Asia on ice thickness change. Sub-regions dominated by winter accumulation (generally westerly regimes), such as the Hindu Kush, Western Himalaya and the Pamir region (Pohl et al., 2015;Yao et al., 2012), show the typical seasonal pattern with mass accumulation during winter/early spring and mass losses in the summer/autumn months ( Figure 5). 345 Contrarily, sub-regions such as Central Himalaya, Eastern Himalaya and Hengduan Shan show a more heterogeneous seasonal pattern. The elevation change time series of these three sub-regions displayindicate that the annual cycle has two peaksmaxima, with a first peakmaximum in winter and a second and smaller peak in summer ( Figure 5, S1). Receiving summer-accumulation through the Indian monsoon these sub-regions generally have a precipitation maximum in July/August, however they are also defined by a high variability of precipitation regimes (Maussion et al., 2014) and a high temperature range (Sakai and Fujita, 350 2017), resulting in glaciers with varying types over very short distances (Maussion et al., 2014) (Maussion et al., 2014). The impact of this variability becomes evident when compared to the more periodic seasonal patterns of the Hindu Kush, Western Himalayas and Pamir time series. This also stands in contrast with the inner Tibetan Plateau, dominated by a more continental climate, which displayswhere glaciers exhibit almost no intra-annual cycle.
In general, the heterogeneity of the time series reflects the sensitivity of mountain glaciers to meteorological patterns and 355 changes and emphasises that glaciers in High Mountain Asia cannot be considered as one entity with uniform temporal variability and sensitivity to changes.

Comparison of regional mass balance with previous work
A comparison of mass balance results in the literature indicates that, while all the studies agree on the general trend in mass 360 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.
Besides the differences in data and methodology, a part of these disagreements can be explained by the time periods.  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.34.2.3
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 380 the Pamir regions. We used the regions by Brun et al. (2017), the RGI 6.0 second order regions and the HiMAP regions (Bolch et al., 2019) to compare our results with other estimates (Figure 7, Figure S8, Figure S9, Table S1, Table S2). 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) (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 395 in monsoonal regimes which could account for all of these changes and attribute the temperature sensitivity of glaciers in 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 (Kääb 400 et al., 2015), moderate (this study), slight mass losses Shean et al., 2020), and even mass gains (Gardelle et 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. 405 Contrasting estimates have also been published for the Pamir and Pamir Alay mountains (Hissar Alay), where high (Kääb et al., 2015), moderate (this study; Ciracì et al., 2020;Gardner et al., 2013), slight mass losses Shean et al., 2020), and even mass gains (Gardelle et al., 2013) have been reported. Part of the discrepancy can be attributed to time variability in mass loss  and driven by fluctuation in winter precipitation (Smith and Bookhagen, 2018).
CryoSat time series indeed suggest near-balance between 2011 and 2015 and increased mass loss from 2015 onwards, which 410 could account for the higher mass loss estimates in comparison to the DEM differencing studies covering the last two decades Gardelle et al., 2013;Shean et al., 2020).
The spatial thinning pattern in the Kunlun-Karakoram area (Figure 2) confirms the suggestion of previous studies Gardner et al., 2013;Kääb et al., 2015) Gardner et al., 2013;Kääb et al., 2015), that the so-called "Karakoram anomaly" Hewitt, 2005) stretches up to West Kunlun Shan, which hasis considered now 415 become the centre of the anomaly. We record less mass gain in Kunlun (+0.01 ± 0.05 m w.e. yr -1 ; +0.06 ± 0.05 m w.e. yr -1 in the Western part of Kunlun) than previous studiesindicating that the Karakoram anomaly might not persist long-term (Farinotti et al., 2020;Rounce et al., 2020). This observation is also reflected in the elevation change profile of the Kunlun regions, where Brun et al. (2017) find constant thickening at almost all elevation during the survey time period of 2000 to 2016, whilst we record thinning at lower elevations (see Figure S7). These findings suggest a shift towards negative mass 420 balance at lower elevations in the Kunlun region in comparison to the previous decade. 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 short-term trend, however, it displays the limitation of all mentioned studies (including this study) when deriving linear trend s in a region like High Mountain Asia with large inter-annual climate variability and associated glacier changes. 425 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.

Temporal variability
The increased thinning from 2013-4 onwards (see Figure 4, 9, S3), which we observed across all sub-regions, has also been reported by Wouters et al. (2019) and Ciracì et al. (2020) in GRACE time series covering the whole Alaska region. This change correlates with the change from a negative PDO phase to a positive phase in 2013-4, which resulted in increased temperatures (Wendler et al., 2017). This is in agreement with (Wouters et al. (2019)This is in agreement with Wouters et al. (2019), finding 435 that their interannual mass change variations negatively correlate with the May-September PDO and May-August NAO indices. Our findings further suggest that the sensitivity of glaciers to the 2013-4 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 elevations within these regions. 440

Comparison of total mass balance with previous work
Our total mass budget of -76.3 ± 5.7 Gt yr -1 (-0.89 ± 0.07 m w.e. yr -1 ) agrees with existing estimates, including those using GRACE (−76 ± 4 Gt yr −1 , −72.5 ± 8 Gt yr −1 and −69 ± 11 Gt yr −1 by Sasgen et al., 2012, Ciracì et al., 2020and Luthcke et al., 2013 and ICESat (−65 ± 12 Gt yr −1 by  as well as a study from airborne altimetry (-75 ± 11 Gt yr -1 by Larsen et al., 2015) and a consensus estimate combining glaciological and geodetic observations (-73 ± 17 Gt yr -1 / −0.85 ± 445 0.19 m w.e. yr -1 by Zemp et al., 2019) (Figure 8b). Our result is significantly more negative than two GRACE studies, with estimates of −53 ± 14 Gt yr -1 (Wouters et al., 2019) and −42 ± 6 Gt yr -1 (Jacob et al., 2012). Besides the variations in methodologies and data between these studies, also differences in study area extents, glacier masks and volume to mass conversion factors contribute to the spread of total mass change results. Our estimates correspond to the RGI 1 region (excluding Northern Alaska) to make the results more comparable for future studies. In general, our total mass balance is more 450 negative than most other studies' findings, reflecting the increased thinning rates we show in the sub-regional time series from 20142013-4.

Comparison of sub-regional mass balances with previous work
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 455 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.  Berthier et al. (2010), providing insights into changes of thinning rates since this time period. Our results are consistently more negative, however the general pattern with the lowest changes discovered in Alaska Range and the highest rates taking place in the Coast Ranges is in agreement with Berthier et al. 460 (2010). We see the largest differences along the east coast − particularly in the Saint Elias mountains − which are also the areas where the lowering of mean surface elevations after 2013-4 has been most pronounced (see Figure 9, S3). 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 changes at the lowest altitudes but less steep gradients in sub-regions along the east coast ( Figure 11, Table S2). This is particularly pronounced in the Saint Elias Mountains, where Berthier et al. (2010) show near-465 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).
In the Western Chugach Mountains, Alaska Range and Alaska Peninsula we observe a decrease of thinning rates towards the 470 lowest elevations of these sub-regions, which can be attributed to the effect of debris cover and the temporal evolution of glacier extent during the study period, one of the limitations of using static glacier masks. This characteristic has been obser ved, although more pronounced and across all sub-regions, by Berthier et al. (2010) and Arendt et al. (2002).

Conclusion
We exploit CryoSat-2 interferometric swath processed data from 2010 to 2019, with a total of 33 million elevation 475 observations, to generate new and independent mass balance estimates for two mountain regions, the Gulf of Alaska (GoA) 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 ratesa rate of 28.0 ± 3.0 Gt yr -1 (0.29 ± 0.03 m w.e. yr -1 ), and the GoA 480 region has lost mass at ratesa rate of 76.3 ± 5.7 Gt yr -1 (0.89 ± 0.07 m w.e. yr -1 ), for a sea-level contribution of 0.078 ± 0.008 mm yr -1 (0.051 ± 0.006 mm yr -1 from exorheic basins) and 0.211 ± 0.016 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 485 Himalayan range and the Tien Shan and slightly positive/near balance trends in the Kunlun and Karakoram ranges, known as the "Karakoram anomaly". The monthly time series of this paper reflect the sensitivity of glaciers in HMA to meteorological patterns and changes and emphasises that the temporal variability of glaciers in High Mountain Asia varies spatially. We show to 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 demonstrates the potential of such a system to monitor trends in ice mass on a global scale and with increased 500 temporal resolution. It also demonstrates the ability to monitor monthly change and paves the way to an observation-based quantification of seasonal accumulation and melting processes, a task that will likely require combination with regional climate models, and with other sensors such as IceSat-2 and high-resolution DEMs.   grid. The size of the circles is scaled by the total glacierised area within a cell. Note that our total mass change estimate of -76.3 ± 5.7 Gt yr -1 (-0.89 ± 0.07 m w.e. yr -1 ) only include glaciers from the RGI region 1 (Alaska). Including also the Northern Rocky Mountains and the Mackenzie and Selwyn mountains we retrieve a mass change of -77.7 ± 5.7 Gt yr -1 .