21st century estimates of mass loss rates from glaciers in the Gulf of Alaska and Canadian Archipelago using a GRACE constrained glacier model

Abstract. Ice mass loss rates from glaciers in the Gulf of Alaska and the Canadian Archipelago are expected to increase through the end of century in response to increasing temperatures. Here, we develop a new glacier model constrained by GRACE gravimetry observations for the period between 2002 and 2017. The high temporal and regional spatial resolution of GRACE mass balance estimates allows us to estimate regional glacier sensitivities to atmospheric changes, and account for higher order of glacier dynamics. We use our regionally constrained models to extrapolate future mass loss under different 5 climate emission scenarios. Generally our 21st century sea level estimates are at the high end compared to other studies. We find that the Gulf of Alaska exhibits the highest mass loss rates between -79 to -112 Gt yr−1 between 2006 and 2100 under different scenarios, and displays the highest sensitivity to the specific scenario (RCP 2.6/4.5/8.5). Our estimates for Arctic Canada South are significantly higher than prior work (-57 to -85 Gt yr−1) and are comparable to projected mass loss rates from the Arctic Canada North (-63 to -101 Gt yr−1). 10

alternate basis on the sphere to spherical harmonics and the functions maximize their energy within a specific region while being orthogonal both on the sphere and region. The localized basis is also sparse, improving signal-to-noise and estimation on timeseries of each coefficient. Slepian functions have been used to obtain the mass loss rates from Greenland (Harig and 90 Simons, 2012), Antarctica (Harig and Simons, 2015), and smaller regions such as Gulf of Alaska, Canadian Archipelago, and Iceland (Harig and Simons, 2016;von Hippel and Harig, 2019). The code to generate the mass loss timeseries are available at Harig and Simons (2015).
Here, we briefly summarize the construction of Slepian functions. When considered over partial sphere regions, R, traditional Spherical Harmonics Y lm are no longer orthogonal; their integral products are no longer delta functions (as is the case over 95 the whole sphere). Instead, the integral products up to a certain bandwidth L form a matrix D, the localization kernel, with significant off-diagonal power, R Y lm Y l m dΩ = D lm,l m . (1) We construct a new set of basis functions, Slepian functions, by taking the eigenvalue decomposition of D as (2) 100 Slepian functions, g, are then the eigenfunctions of localization kernel. The functions are constructed as a combination of spherical harmonic functions. The eigenvalues, λ, vary between zero and one, and describe how well each eigenfunction is concentrated in the region. By considering only large eigenvalues that maximize the concentration ratio as we can create a new truncated basis set of only the most concentrated functions, which are orthogonal both on the whole sphere 105 Ω and over the region R. The Shannon number, where N describes the effective number of functions that can be concentrated within a region for a given bandwidth (L) and region size (A). Practically, this includes functions with eigenvalues λ > 0.5. We then use this truncated set of functions in traditional signal estimation where 110ŝ (r) = N α=1ŝ α g α (r), and g = L lm g lm Y lm .
The local signals in GRACE data are efficiently recovered by transforming spherical harmonics Y lm (r) into Slepian functions has been shown to obtain higher spatial resolution information from the existing monthly Level-2 GRACE data products than other methods, with clear estimates of noise uncertainty Simons, 2012, 2016).

Glacier mass balance model from climate data
The mass balance model is developed from the use of input gridded climate reanalysis products, global glacier inventory from the Randolph Glacier inventory (RGI) 6.0 (RGI et al., 2017) and initial assumptions in volume area scaling. The modeling pro-120 cess is based on a temperature-index degree day calculation, and it broadly includes the following components: (a) computation of ablation and accumulation from ERA-Interim temperature and precipitation for the observational period between 2002 and 2017, (b) inputs of hypsometry such as the glacier outlines, area and elevation for volume and mass balance calculations, (c) glacier assumptions for the initial conditions of area and volume. First, we describe the details of the datasets used in modeling, followed by the computation of regional mass balance. 125 We estimate monthly intervals of glacier mass balance from 2002 to 2017 using the temperature and precipitation products from ERA-Interim (ERAI) reanalysis datasets (Dee et al., 2011;ECWMF., 2012). ERAI is a global gridded data assimilation product, provided by the European Center for Medium-Range Weather Forecasts (ECMWF) at a spatial resolution of 0.7 • from 1979 to 2018. We use ERAI geopotential temperature and ERAI precipitation data products in the calculation of ablation and accumulation, respectively, described further below.

130
The RGI contains glacier information for 563 glaciers from the Gulf of Alaska North, 4000 glaciers in the Canada North and 7413 glaciers in the Arctic Canada South. The glacier outlines are provided at a spatial resolution of 0.5 • x 0.5 • . To maintain the consistency between the RGI regional outlines and GRACE regions of interest, we include only the glaciers defined by RGI regional boundaries. For Gulf of Alaska, we consider glaciers from the Gulf of Alaska North region due to the regional difference in mass loss 135 rates compared to the Gulf of Alaska South region (Arendt et al., 2002). RGI provides hypsometric details of glacier area at an elevation interval of 50 m. From this information, we compute the area elevation distribution of glaciers at every 50 m grid spacing. We compute the regional mass balance based on elevation bins and individual glaciers defined by the RGI.
We use the degree day approach to model the present and future mass balance for transient conditions of temperature and precipitation for glaciers in the Gulf of Alaska and Canadian Archipelago. For every glacierized cell of RGI at 0.5 • spatial 140 resolution, monthly rates of ablation and accumulation are calculated based on the relationship between local meteorological condition and glacier hypsometry Hock, 2006, 2011). The mass balance (M ) of a glacier is given as the sum of ablation (A), accumulation (C) and refreezing (R) as the regional mass balance (Wahr et al., 2016).
We calculate ablation for each glacier, A gl , from ERAI temperature at degree day (DDF , mm.w.e.d −1 C −1 ), which is typically the number of melt days above the threshold temperature T o , to differentiate the solid and liquid precipitation, as Surface temperature or temperature at the glacier, T gl , is calculated from downscaling ERAI geopotential temperature based 150 on the lapse rate, lp, where h is the mean elevation of an individual glacier, and ∆h is the average elevation of glaciers in a region, Accumulation is based on the precipitation, which is the amount of snowfall considered as ice. Like temperature, the precipitation, P gl , is downscaled from convective precipitation P (t) using a precipitation gradient d prec , and average elevation of 155 glaciers in the region ∆p, There is a slight difference in the temperature and precipitation lapse rates because the precipitation is not accurately represented at high altitudes, especially from the climate reanalyses data products (Behrangi et al., 2016). Accumulation from glaciers C gl (t) is then calculated as, where K o is the precipitation correction factor that is used to account for snowfall from all the glaciers.
While modeling the mass balance of individual glaciers in the RGI, the accumulation and ablation rates are computed by matching each glacier to the nearest ERA-Interim grid cell. We ignore the effects of tidewater calving from Gulf of Alaska and Canadian Archipelago since they contribute less to regional mass balance (Larsen et al., 2015). The mass balance is calculated 165 in two ways, based on elevation bins within the region of interest, and by considering each individual glacier. We find that the computation of mass balance from glaciers is useful for model tuning, contrary to the elevation binned mass balance as in previous glacier modeling studies (Huss and Hock, 2015). We also found that the results from these two methods does not vary much.

170
Our glacier model derived using ERAI dataset is highly sensitive to the primary parameters that include the degree-day factor (DDF ), threshold temperature (T o ), precipitation correction factor (K o ) and less sensitive to the temperature lapse rate (∆h), precipitation lapse rate (∆p), precipitation gradient (d prec ), and environmental lapse rate (lp) ( Table 5). In order to constrain our model to an observational record, we use the monthly timeseries of regional mass balance from GRACE (Section 2.1) to compare against our monthly modeled glacier mass balance. We limit our model space to solve for three primary parameters, 175 DDF, T o and K o (bold in Table 5), and fix the remaining parameters to values used in previous studies (e.g. Wahr et al., 2016).
The primary parameters have greater sensitivity for controlling the modeled mass balance. The model is run for individual glaciers within our region of interest, and it is integrated for all glaciers to obtain total mass balance. We then use a grid search optimization using the parameter ranges listed in Equation (11) to solve for the combination that best minimizes the residual sum of squares between modeled mass balance and GRACE observed mass balance, calculated as One key difference between our model and prior work is that we calculate the optimal model parameters (DDF , T o , and K o ) in each region separately because we have separate regional GRACE mass balance timeseries.
The modeled mass balance explains nearly 87% to 88% of variance from the GRACE mass balance (depending on region, see figure A3) using three parameters mentioned in Table 5. Even though the GRACE regional mass balance would not be 185 expected to represent the changes in any specific individual glacier, we compare our modeled mass balance with several individual glaciers that have direct observations in the WGMS in Sections 3 and 4.3.

Future mass and volume rates
We use the optimized mass balance model developed using ERA Interim to predict future mass and volume loss rates until 2100.
Climate model data from CESM-LE , HadGEM-ES , and MRI-CGCM3 (Yukimoto 190 et al., 2012) was used for transient conditions of temperature and precipitation under several different emission scenarios. Table 2 shows the resolution and details of climate model data. Prior to future projections, we adjusted the bias in temperature T i (t) and precipitation P i (t) between ERA-Interim and each climate model using a delta approach (Radić and Hock, 2006) for the baseline period between 2002 and 2017 as where i = 1, 2, ...12 for each month and c refer to the model/emission scenario combination. The spatial resolution of grid cells near the poles was accounted in the computation of mean temperature and precipitation within our region of interest.
We found that the temperature pattern from both the regions (Gulf of Alaska and Canadian archipelago) closely modeled the climate model data, whereas the precipitation pattern had a slight bias after the correction. In order to incorporate glacier area changes to 2100, we follow the volume-area scaling method (Radić and Hock, 2011;Wahr et al., 2016). We use the initial conditions of glacier area and volume-area scaling (power law) relationship, to determine the glacier area and volume evolution over time (Bahr et al., 1997;Radić and Hock, 2010;Wahr et al., 2016).
Here, volume V relates to S the glacier area, where c is the constant in units of length raised to the power (3-2γ ), and γ is the 205 dimensionaless scaling component based on theortical assumptions by Bahr et al. (1997). For glaciers in Gulf of Alaska, we consider γ = 1.36 and use γ = 1.34 for ice cap as in Arctic Canada North and South. The value of c = 0.03 to 0.026 under the density assumption of ρ ice (Bahr et al., 1997(Bahr et al., , 2015. Following this, we use the relationship between volume and mass based on the density of ice (ρ ice = 917 kg m −3 ) to calculate volume rates as where V gl is the volume of glaciers and M gl (t) is the total mass balance of glaciers computed from equation 6. Here, our model does not incorporate the time series evolution of glacier area changes, since our model constrained by GRACE observations has secular and seasonal trends in mass balance. Therefore, we incorporate certain assumptions in the initial state of glacier geometry. In order to test the effects of change in mass and volume loss based on steady state of glacier area and hypothetical 215 conditions of mass loss, we perform a synthetic test and remove mass from glaciated area below < 1500 m or < 600 -800 m (Gulf of Alaska and Canadian Archipelago) to understand the effects of glacier hypsometry.
In each model run for the period from 2002 to 2017, the model is tuned with three sensitivity parameters such that we obtain maximum variance between GRACE and modelled mass balance. For extrapolation of volume and mass loss rates, the main three parameters are kept constant while tuning the other sensitivity parameters such as lapse rate, precipitation factor and 220 precipitation lapse rate (Ref table 5). For initial conditions of volume, we considered 1.6 x10 4 km 3 for Gulf of Alaska North, 0.97x10 4 km 3 for Arctic Canada South, and 3.2 x10 4 km 3 for Arctic Canada North ( (Wahr et al., 2016)). The new rates of mass balance ∆M gl are obtained for glaciers using the relationship between initial mass M gl 0 and change in mass ∆M gl (t) as, The sea-level estimate (SLE) or change in sea-level is computed following where we assume the area of ocean as 361 x10 6 km 2 .  Table 2, our model is able to explain about 87 to 88% of variance from GRACE observation, depending on region (See figure A3).
The values from correlation and RMSE indicates that our model was able to explain the cumulative mass change obtained from GRACE data. In Figure 1d (Table 5).
While we optimized our model against the regional mass balance, we also compared our model output against direct observations of mass balance from glaciers or ice caps in the Gulf of Alaska and Arctic Canada North ( Figure A4). Gulkana and

240
Wolverine glacier in the Gulf of Alaska region have direct mass balance observations that overlap the GRACE observational period ( Figure A2). Gulkana glacier appears to be losing more mass than the Wolverine glacier, and our model explains 45% of the variance with significant p-value. The Wolverine glacier model explained less variance with insignificant confidence levels ( Figure A5). In the Arctic Canada North, we found that direct observations from Meighen and White glacier explained significant relationship with variance up to 84% and 45%, however, the observations from Devon ice cap was insignificant with 245 45% correlation ( Figure A6). In the above comparison of direct observations, the three main sensitivity parameters were kept constant as mentioned in Table 2 (highlighted). Optimization of direct observations was achieved mainly from tuning the lapse rate in temperature and precipitation, which nearly represents the average elevation of glacier. This experiment indicated that our regional model was able to represent local observations.  Our observations of volume and mass loss rates derived from the optimization of GRACE dataset are consistent with sev-260 eral recent studies and in accordance with the IPCC (2019). We adopted a unique approach in modeling mass balance rates compared to the existing methods that are highly dependent on direct observations for extrapolating regional and global mass balance. We discuss the results from the perspectives of (i) model comparison with previous studies, (ii) model performance and sensitivity analysis, (iii) model representation with the local observations, and (iv) uncertainties in projected estimates of SLE and error propagation.

Model comparison with previous studies
Even though the modeling framework adopted in our study is not directly comparable to the existing glacier models, the results are in accordance with the published sea-level estimates in the IPCC special report (IPCC, 2019). In the glacier model intercomparsion project (GlacierMIP), six glacier models have been compared in terms of their modeling capabilities and uncertainties in the future projection. Like the existing glacier models, we used a degree day method to model the mass balance 270 from gridded temperature and precipitation. In the previous studies, either individual glaciers from the RGI inventory or the total glacier area grouped in elevation bins were calibrated with WGMS direct observations for all regions together (Marzeion et al., 2012;Giesen and Oerlemans, 2013;Hirabayashi et al., 2013;Radić et al., 2014). Based on the model calibration with direct measurements, the volume and mass loss rates were extrapolated for future projections. One of the drawback from this approach was large uncertainties in the regions of sparse direct observations or high topographic regions such as the High

275
Mountain Asia or Andes. Some of the results from these models could not be compared directly due to the differences in glacier inventory version and meteorological inputs from either ERA interim or CRU datasets.
This problem was circumvented by the use of standardized regional mass balance by Gardner et al. (2013) to calibrate the model (Huss and Hock, 2015). In this way, every glacier in the RGI inventory is tuned to match the regional mass balance between 2003 and 2009. This method does not attempt to model the local observations from individual glaciers but does solve 280 the uncertainties from regions deficient in direct observations. Model calibration also included glacier thickness measurements for all glaciers to incorporate dynamics in the future projections of volume and mass loss, instead of volume area scaling. Huss and Hock (2015) also included calibration for grounded ice that is already displaced into ocean from calving fronts. It was estimated that error from displacements due to calving fronts contributed about 10-14% to the global sea-level.
One of the advantages from our model is that it does not require calibration from direct mass balance observations, which 285 are spatially and temporally sparse, with annual field observations. Instead, we calibrate our model using GRACE observations with higher spatial and temporal resolution, which are known for estimating regional mass balance very accurately (e.g. Harig and Simons, 2012). The regionally calibrated parameters (DDF , T o and K o ) are unique for Alaska, Arctic Canada North and South, and they are capable of representing direct observations from individual glaciers ( Figure A4). In this way, we address the bias due to sampling issue from direct observations.

290
In the existing models, calibration did not account for dynamics from glaciers which could reflect the response time of glacier leading to cumulative bias for longer timescales. Only the model by Marzeion et al. (2012) included relaxation time in glaciers for the present and future rates of volume and mass loss. Our model constrained by GRACE intrinsically accounts for higher order dynamics, evident from the seasonal signal at monthly intervals (Fig 1). It is also one of the reason that we obtain South, and 9-12% for Arctic Canada North (Table 3). Our observations are somewhat consistent with Huss and Hock (2015) for Arctic Canada South, where the volume loss is -37 ± 14 % for RCP 4.5, while our observations indicated a volume loss of 300 ∼20%. We do not attempt to compare the initial state of glacier thickness from other studies due to different versions of RGI inventory and assumptions in volume-area scaling.
We tested the sensitivity of our model with three experiments that influences the mass and volume rates. In the dynamic adjustments of glaciers for mass loss, we assume that the glacier will lose mass near the terminus, therefore we remove 10% of mass from the terminus region ( Figure 3). We also tested with a temperature increase of 1 • K and precipitation change of 305 +10%. Our model is highly sensitive to temperature changes with mass loss and SLE increase by 11-28 Gt yr −1 and 6-11 mm, compared to the normal conditions (Table A2). Comparatively, our model was less sensitive to changes in precipitation.
We compare our results to the projection estimates from Glacier model intercomparsion (GlacierMIP) project that incorporates six published global glacier models (Hock et al., 2019) ( Figure 5). Some of the previous models did not attempt mass balance projection for RCP 8.5 or RCP 2.6, therefore we use RCP 4.5 as an intramodel comparsion and later, we compare Gt yr −1 ) from 1996 to 2015 (Noël et al., 2018). This indicates that the mass balance response has increased over the last decade in both Arctic Canada North and South in contrast to the inference from glacier models that reported smaller estimates of SLE for Arctic Canada South . We think that global models initialized from direct observations and inventories does not capture the sensitivity from small glaciers and ice caps, as in Arctic Canada South. There is an uncertainty of 18-330 30%, specifically due to sensitivity from small glaciers (Huss and Hock, 2015;Hock et al., 2019). This leads to the erroneous conclusion that large glacierized areas in the Arctic tend to lose more mass than small glaciers in the Arctic Canada South.

Model performance and sensitivity analysis
We consider a different approach to determine the glacier volume and mass loss rates incorporating higher order of glacier dynamics without the need for inputs from direct observations or glacier thickness measurements (Radić et al., 2014;Huss 335 and Hock, 2015). In the Figure 1d and f, the modelled rates of mass balance agree well for Gulf of Alaska and Arctic Canada North, however there is a discrepancy in the mass loss rates for Arctic Canada South. The modelling parameters indicated in Table 5 achieved good correlation and model convergence for all three regions, however there is slightly a large discrepancy in mass loss rates between GRACE and model for Arctic Canada South. We could have acheived lower discrepancy between the GRACE and modelled time series for ACS, but this will involve compromising some of the model parameters for convergence.

340
This could be unrealistic, for example, our model cannot have precipitation lapse rate > 1500 m as it does not represent the ELA of glaciers in ACS. Also, we found that using these unrealistic values in the model parameters significantly affected in the future projections of mass loss rates and SLE. We agree that our estimates of GRACE mass loss rates were slightly lower at -27 Gtyr −1 , compared to recent estimate of -32 ± 8 Gtyr −1 (Wouters et al., 2019). We also note that the seasonality and inter-annual variability present in GRACE data is not well represented in our model mass balance time series (Figures 1 d-f).

345
GRACE is well known for capturing the seasonality in mass balance, whereas the modeled mass loss rates from reanalysis products often have large biases occurring from poor representation of precipitation in cold regions, either at high latitude or high elevation (Behrangi et al., 2016).
Previously, GRACE observations has been used in modeling sub-annual mass balance from glaciers in the Gulf of Alaska with high degree of fidelity in modeling local observations Luthcke et al., 2013). Monthly GRACE 350 observations track the seasonality from glaciers, up to a variance of 75% and 56% from summer and winter balances respectively . We use this rationale to model the mass balance from individual glaciers, constrained by GRACE observations.
While we employed seven sensitivity parameters to model the volume and mass loss rates from three regions (See Table   5) we optimize our model based on three main parameters (DDF , T o , and K o ) to constrain our model against GRACE mass 355 balance timeseries (Wahr et al., 2016). This method eliminates the need for extrapolation of direct observations for regional mass balance and SLE as in Radić et al. (2014). Further, we find that the degree day factor (DDF ) and threshold temperature (T o ) varies between the three region. The DDF for glaciers in Gulf of Alaska is 2.3 mm d −1 C −1 , whereas for Arctic Canada North and South, the DDF is slightly higher with values of 3 mm d −1 C −1 and 4 mm d −1 C −1 , respectively. This is due to higher mass balance sensitivity of glacier with temperature in the Arctic regions (Hock, 2003;De Woul and Hock, 2005;360 Braithwaite and Raper, 2007). Similarly, the threshold temperature T o varies depending on the average temperature of the region, which is larger for Arctic regions compared to the Gulf of Alaska. By optimizing these three parameters, our modeled mass balance was able to explain ∼ 90% of the variance from GRACE observations. Similar to Radić et al. (2014), we found that our model is sensitive to environmental lapse rate (lp) and temperature lapse rate (∆h). The values for environmental lapse are based on Braithwaite (2008). We introduced and −0.24 and −0.19 for Arctic Canada South (Baffin). Their smaller value in bias is due to extrapolation of regional mass 370 balance from a sample of ∼ 20-34 glaciers from these regions and a different version of RGI inventory. In contrast, GRACE regional mass balance includes observations from all glaciers with proper representation of summer and winter balance without the need for scaling from direct observations.

Comparison with local observations
In the Gulf of Alaska, Gulkana and Wolverine glacier have direct mass balance data that overlap with our observational period 375 from 2002 to 2017. In order to model the local observations, the three main sensitivity parameters were kept constant, while adjusting the other parameters until the model bias is minimized (Table 5). Table A3 Wouters et al., 2019). These studies, like our model indicated that the temperature is a contributing factor for glacier melt in the Gulkana glacier. Our model did not exhibit good correlation for Wolverine glacier due to its proximity near the coast, with mass loss rates controlled by precipitation. The result was similar to the approach where individual mascons from GRACE observation were compared with direct observations .

385
In the Arctic Canada North, our model was able to represent in-situ observations from Devon ice cap (r = 0.42, p = 0.13), Meighen ice cap (r = 0.84, p = 0.001) and White glacier (r = 0.44, p = 0.09) ( Figure A6). Most of the glaciers and ice caps are known to be losing mass rapidly in the 21st century based on a series of in-situ and satellite observations (Bezeau et al., 2013;Noël et al., 2018;Cook et al., 2019). Agreement of modelled mass balance is consistent for Meighen glacier, with slight inconsistency for Devon ice cap owing to different input datasets and modeling techniques (Lenaerts et al., 2013). These 390 observations indicate that our model constrained by GRACE is capable of representing in-situ mass balance.

Conclusions
We present estimates of glacier volume and mass loss from a GRACE constrained glacier model that can be compared with the recent glacier model intercomparison project (GlacierMIP). We have three key findings from our glacier model and they are as follows. (a) We have demonstrated that the regional bias can be significantly low compared to calibration from direct 405 observations, during the observational period. This is due to the use of GRACE solutions for regional mass balance. This assumption holds true for any RGI glacier region, especially with HMA, Andes or Caucasus with complex topography. (b) Our model is able to represent local observations with good level of confidence in the mass balance timeseries. The unique sensitivities values obtained for each region represents the mass balance change from individual glaciers. This is useful for understanding the response of benchmark glaciers under the present and future conditions of warming without the need for 410 conventional field observations.
(c) We find projected mass losses from the Arctic Canada South that are far higher than other studies, and which are on par with losses from the Arctic Canada North. This indicates that Arctic Canada South has greater sensitivity in the recent decade, and our model is able to capture this sensitivity. In order maintain the consistency with global glacier models, we have used