Glacier changes in Alaska: can mass-balance models explain GRACE mascon trends?

Abstract

Temperature and precipitation data from three weather stations in the St Elias Mountains of Alaska and northwestern Canada were used to drive one-dimensional (1-D) (elevation-dependent) and 0-D degree-day mass-balance models. Model outputs were optimized against a 10 day resolution time series of mass variability during 2003–07 obtained from Gravity Recovery and Climate Experiment (GRACE) mass concentration (mascon) solutions. The models explained 52–60% of the variance in the GRACE time series. Modelled mass variations matched the phase of the GRACE observations, and all optimized model parameters were within the range of values determined from conventional mass-balance and meteorological observations. We describe a framework for selecting appropriate weather stations and mass-balance models to represent glacier variations of large regions. There is potential for extending these calibrated mass-balance models forwards or backwards in time to construct mass-balance time series outside of the GRACE measurement window.

Introduction

Mass-balance models driven by temperature and precipitation data from local weather stations or climate model output are commonly used to estimate glacier contribution to rising sea level (Hock, 2005). Confidence in such models improves when they are calibrated against and agree with mass balances determined from field, satellite or airborne measurements. In many cases, data are available from a small number of benchmark glaciers, but the parameters derived from these glaciers may not represent entire glacier regions. Improved models are necessary for predicting future changes of mountain glaciers, expected to be important contributors to rising sea level in the next 50–100 years (Meier and others, 2007).

Data from the Gravity Recovery and Climate Experiment (GRACE) measure Earth mass variations which, when corrected for non-glacier sources of mass variability, yield estimates of the magnitude and distribution of recent (2003–present) mountain glacier mass trends (Tamisiea and others, 2005; Chen and others, 2006a). Luthcke and others (2008) have developed a high spatial (2 ^˚ × 2 ^˚ ) and temporal (10 day) resolution GRACE solution procedure that reduces GRACE data to temporally varying local mass observations. This procedure has been applied to the glaciers in the Gulf of Alaska region, yielding trends in glacier mass during 2003–07. A subset of GRACE mass concentration (mascon) trends from September 2003 to August 2007 agreed to within 5% of concurrent aircraft altimetry mass-balance estimates in the St Elias Mountains, Alaska, USA, and Yukon, Canada (Arendt and others, 2008). Such validation work illustrates the capability of GRACE mascon solutions to resolve multi-year glacier mass trends of large and dynamic mountain glacier systems.

Degree-day models generally provide robust estimates of glacier mass changes for individual glaciers (Jóhannesson and others, 1995; Braithwaite and Zhang, 1999; Hock, 2003). It has not been determined whether such models, calibrated to individual glaciers, capture regional mass variations of entire glacier systems. Here we investigate techniques for calibrating a regional degree-day mass-balance model against GRACE mascon solutions for glaciers within and surrounding the St Elias Mountains, a region with 39 000 km2 of glacier ice. We design a modelling strategy that accounts for elevation-dependent variations in snowpack evolution and glacier melt. We also test a simplified zero-dimensional (0-D) model that relates positive degree-days and snow accumulation to the regional balance, without accounting for elevation variations of the input parameters. Both models are optimized against the 10 day resolution GRACE mascon time series, using a total of four parameters to adjust temperature, precipitation and melt factors that convert positive degree-days to melt equivalents.

Our goals are to: (1) provide a precise interpretation of the information contained within the GRACE mascon solution; (2) assess the extent to which variations in mass determined from GRACE can be modelled by time variations in temperature and precipitation; and (3) determine the optimum level of model complexity required to reproduce the GRACE mascon series from nearby climate station data. Our analysis will provide a framework for extending the GRACE time series beyond the existing measurement window.

Methods

Grace mascon solution

Estimates of the mass variations of glaciers surrounding the Gulf of Alaska have been obtained from high-resolution GRACE mascon solutions (Luthcke and others, 2008). A series of background models were used to remove gravity signals due to glacial isostatic adjustments, Earth/ocean tides and hydrospheric/atmospheric variations. Gravity signals due to terrestrial water storage on land surfaces were removed using Global Land Data Assimilation System (GLDAS)/Noah 0.25 ^˚ , 3 hour resolution data (Ek and others, 2003; Rodell and others, 2004). Terrestrial water storage corrections were set to zero for glacier areas, which were delineated using a 0.25 ^˚ resolution binary glacier mask calculated from Digital Chart of the World glacier outlines (Raup and others, 2000). After background model corrections have been applied, we assume that the residual gravity signals represent the glacier component of mass variations.

The exceptional spatial (2 ^˚ , approximately 49 000 km2) and temporal (10 day) resolution of our GRACE solution procedure is attained by preserving the gravity information contained within the GRACE intersatellite range-rate measurements, and by parameterizing local mass variations as mascons (Rowlands and others, 2005; Luthcke and others, 2006). A set of 12 mascons were defined for the glacier regions, and the solution error for each mascon was 1–2Gt a ⁻ ₁. Each solution provides an estimate of the time-averaged mass within the glacier system, referenced to the mass from the background models described above, over each 10 day interval.

Arendt and others (2008) selected a subset of high-resolution GRACE solutions from Luthcke and others (2008) representing mass variations in the St EliasMountains (Fig. 1). They validated a subset of the 4 year trend in glacier mass from GRACE using independent aircraft laser altimetry measurements. The close correspondence between these independent regional estimates (−20.6 ± 3.0Gt a ⁻ ₁ vs −21.2 ± 3.8Gt a ⁻ ₁ for GRACE and altimetry, respectively) suggests GRACE accurately represents multi-year glacier mass trends in the St Elias Mountains.

Fig. 1. Location of glaciers (shown in black) of the Gulf of Alaska region and (b) St Elias Mountains, Alaska USA, and Yukon, Canada. Black boxes outline the 2 ^° mascons defined in Luthcke and others (2008). Cordova, McCarthy and Yakutat mark the locations of three weather stations used to drive mass-balance models. Mascons 6, 7 and 10 contain 19, 42 and 17% glacier ice, respectively, relative to the total area of the mascon.

Based on these results, we are most confident in the GRACE solutions for glaciers of the St Elias Mountains, and therefore focus our present analysis on this region. In contrast to the study by Arendt and others (2008), where subsetting of the GRACE solution was necessary for comparison with available aircraft altimetry data, here we include all 39 000 km₂ of ice in mascons 6, 7 and 10. Glacier areas in these mascons were determined from manual digitization of Landsat 7 Enhanced Thematic Mapper Plus (ETM+) imagery acquired during 1999–2001 (US National Snow and Ice Data Center/World Data Center for Glaciology, http://nsidc.org/glims/). Detail on the calculation of glacier hypsometries can be found in Arendt and others (2008). GRACE time series for mascons 6, 7 and 10 were added for each 10 day interval to obtain the mass variations in our study region in units of Gt (Fig. 2). Mass values were arbitrarily set to zero at the start of each GRACE solution.

Fig. 2. Sum of GRACE high-resolution solutions (mascons 6, 7 and 10) for the period September 2003–August 2007, GRACE observations (black dots, left axis) representing 10 day averages of the cumulative glacier mass balance, with error bars indicating 1σ uncertainties. The smoothed fit to the GRACE observations using a Gaussian filter with a 10 day window (black line), and trend in mass (grey line) recovered from simultaneous estimation of bias, trend, annual and semi-annual sinusoid, are shown. Diamonds (right axis) are the period-to-period differences of the smoothed GRACE observations, representing the 10 day glacier mass balance.

Temporal averaging is necessary to align mass-balance estimates from models or ground measurements with GRACE observations (e.g. Swenson and Wahr, 2006). Here we express the GRACE time series in glaciological terminology, and describe the necessary transformations between mass-balance model output and GRACE. We consider the ensemble of glaciers in mascons 6, 7 and 10 as a single glacier system. Our glacier mass-balance models described below estimate the daily mass balance B of this system in units of Gt. The summation of B from the beginning of the time series to time t (days after the start of observations) is the cumulative mass balance:

(1)

Unlike conventional mass-balance measurements that sample mass variations at discrete points in time, each GRACE measurement (black dots in Fig. 2) represents a time average of the cumulative mass balance during consecutive 10 day periods:

(2)

where k is the index associated with each GRACE observation. For the purpose of plotting, we locate GRACE observations at the midpoint of each 10 day interval, at time t = (10k − 4.5).

Although we are primarily interested in determining the trend in B_t ( ^_cml ) for sea-level studies, we fit our model instead to the period-to-period change in time-averaged mass. For simplicity we call this the 10 day mass balance B ₁₀, even though it is not the change in mass over the 10 day period but rather the difference between two time-averaged mass values:

(3)

in units of Gt, shown as the diamonds in Figure 2. By fitting our model to B ₁₀ we attenuate the variance in the low-frequency domain of the cumulative series, and ensure stationary mean and variance through each season. This results in more robust statistics of model fit.

Model Descriptions

Model 1: elevation-dependent degree-day model

This 1-D model assumes that spatial variations in mass balance can be approximated as a linear function of elevation. We define B as the daily mass balance b at a specific elevation z, integrated across the area distribution function a(z):

(4)

where Z is the elevation range of the area distribution. We model the glacier surface mass balance b using a temperature-index (positive degree-day)model (Hock, 2003):

(5)

(6)

(7)

where T is the daily average air temperature ( ^˚ C), P is the daily total precipitation (rain and snow, mw.e.), DDF_snow/ice is the degree-day factor for snow/ice (m ^˚ C ⁻ 1 d ^{− ₁} ), z is elevation (m), Δt = 1 day, and the subscript ‘met’ refers to values measured at the weather station. Values of T _met and P _met are adjusted for elevation using constant temperature and precipitation lapse rates ΓT ( ^˚ Cm ⁻ 1), ΓP (% m ^{− ₁} ) and a correction factor for precipitation k to account for differences in precipitation regime between the weather station and the glacier regions. δ determines the threshold between positive temperatures for melt and negative temperatures for accumulation of solid precipitation:

(8)

We define a(z) in 10 m elevation bins, and solve Equation (4) at the median elevation of each bin.

Model 2: zero-dimensional model

We introduce a simplified approach to modelling B that estimates the balance at a single elevation. Our approach is similar to linear regression models that relate meteorological parameters to mass-balance observations (e.g. Eisen and others, 2001; Rasmussen, 2004). We are motivated to test this simpler model for two reasons. (1) The GRACE mas-con time series contains no information on the elevation variability of mass change, so that our elevation-dependent model 1 is not well constrained. (2) There exists an elevation on each glacier at which the mass balance at that elevation is equal to the mass balance B averaged over the entire glacier (e.g. Rasmussen, 2004).

We assume a linear function of balance with elevation, in which case the point balance is equal to the average balance at the glacier’s mean elevation z̄. We therefore use Equations (6) and (7) to adjust temperature and precipitation values to the elevation of z z̄ = 1600 m, the mean elevation of all glaciers in the St Elias Mountains, using 0.1%m ⁻ 1 for ΓP. The value of ΓT is determined through least-squares optimization. We calculate B as a linear function of daily positive degree-days (PDDz̄) and daily total snow precipitation (SNOWz̄):

(9)

Input data and model optimization

All simulations use daily average temperature and total (rain and snow) precipitation data from three US National Oceanic and Atmospheric Administration (NOAA) climate stations, during the period 1 August 2003 to 31 August 2007 (Fig. 1). Optimization of both models is achieved using a non-linear least-squares fitting algorithm (Levenberg–Marquardt). Each model has four adjustable parameters: DDF_snow, DDF_ice, k and ΓT for model 1 and c1, c2, c3 and ΓT for model 2. Previous studies used a precipitation lapse rate of ΓP = 0.05–0.1%m ^{− 1} (e.g. Jóhannesson and others, 1995); we choose a fixed value of ΓP = 0.1%m ⁻ 1 because we expect steep precipitation gradients to occur in this primarily maritime region. We calculate B with each model and then determine time-averaged mass and 10 day mass balances B ₁₀ using Equations (1–3).

Next, the best fit is determined between modelled and GRACE observations of B ₁₀. Finally, we calculate the modelled cumulative mass balance using Equation (1) for comparison with GRACE. From this, an assessment is made of each model’s ability to simulate secular trends in glacier contribution to sea-level change.

Model Results

Model 1

The elevation-dependent degree-day model captured 52–58% of the variance in B10 determined from the GRACE time series, with root-mean-square (RMS) errors (normalized by the variance in the GRACE observations) ranging from 0.64 to 0.71 Gt, depending on the choice of weather station (Table 1).

Table 1. Optimized parameters obtained using temperature and precipitation data from three weather stations to drive an elevation-dependent degree-day model (model 1). Model output is calibrated against a 10 day resolution time series of mass variations from GRACE mascon solutions. Root-mean-square error (RMSE) is normalized by the variance in the GRACE observations

Of the three weather stations, Yakutat had the highest r ^₂ and lowest normalized RMS error. Examination of the observed and predicted B ₁₀ time series shows that both the McCarthy and Cordova simulations underestimated mass gains during the accumulation seasons (Fig. 3a). Yakutat and Cordova, both located near the coast, had k values from 0.21 to 0.32, indicating measured precipitation values were too high by a factor of about 4. These stations receive a considerable amount of rainfall, and it is possible that localized rain events at the weather station did not occur as snow at high elevations on the glaciers. The model may also have required a reduction in precipitation because of the large amount of glacier ice at high elevations, where a linear lapse rate in precipitation created too much snow. Degree-day factors ranged from 0.0009 to 0.0026m ^˚ C ^{− 1} d ⁻ 1 for snow and 0.0041 to 0.0065 m ^˚ C ^{− 1} d ⁻ 1 for ice, all within the range values listed in previous studies (Hock, 2003).

Fig. 3. Model 1 and (c, d) model 2 estimates of 10 day resolution mass variations of glaciers in the St Elias Mountains determined from GRACE. The model simulations are calibrated to the 10 day mass balance. In the cumulative mass-balance series, every third data point is plotted for clarity.

We calculated B _⁽cml) and compared total mass losses at the end of each time series. Simulations using Yakutat data agreed best with GRACE observations, with cumulative mass losses that were 1.2 times greater than GRACE observations (Fig. 3b). Mass losses simulated using Cordova and McCarthy data were approximately double those predicted by GRACE, due to the inability of these stations to reproduce the necessary snow accumulation.

Model 2

Our simplified 0-D model captured 54–60% of the variance in the GRACE mascon time series for the three weather stations (Table 2; Fig. 3c). Normalized RMS errors ranged from 0.63 to 0.67 Gt. Simulations using data from Cordova had a slightly better fit than Yakutat; however, all three stations produced mass variations whose amplitude and phase agreed well with the GRACE observations (Fig. 3d). Values of Γ _T ranged from –0.0034 to –0.0049 ^˚ Cm ^{− 1}, about half the magnitude of Γ _T calculated in model 1.

Table 2. Optimized parameters obtained using temperature and precipitation data from three weather stations to drive a simplified 0-D model (model 2). Model output is calibrated against a 10 day resolution time series of mass balance from GRACE mascon solutions. Root-mean-square error (RMSE) is normalized by the variance in the GRACE observations

The coefficients c ₂ and c ₃ describe the balance sensitivity to changes in positive degree-days and snow precipitation at the mean elevation. These parameters are similar to mass-balance sensitivities to temperature and precipitation (S _t and S _p respectively) used to predict glacier response to long-term climate changes. For comparison with published sensitivities, we perturbed the temperature and precipitation time series for the Yakutat weather station by 1 ^˚ C and 10% respectively, and determined the average annual change in PDD_z and SNOW_z. We multiplied changes in PDDz and SNOW_z by c ₂ and c ₃ respectively, and divided by the total area of glaciers in our solution domain to obtain the area-averaged change in balance. These calculations yielded St = −0.45m ^˚ Ca ^{− 1} and S _p = 0.039ma ^{− 1} (10%) ^{− 1}. These are on the low end of published values S _t = −0.65 and –0.84m ^˚ C ^{− 1} a ^{− 1} and S _p = 0.04 and 0.23 ma ^{− 1} (10%) ^{− 1} for Gulkana and Wolverine Glaciers, respectively (De Woul and Hock, 2005).

Model comparison

Model 2 produced a better fit to the GRACE observations for all weather stations. We attribute this to the fact that we lack sufficient information, in particular balance variations with elevation, with which to constrain model 1. For both models, optimization produced parameters that are physically realistic and within the range of values determined from conventional mass-balance estimates. We note that our optimization procedure minimizes variations in amplitude but makes no phase adjustments. It is therefore encouraging that our models generally capture the seasonal patterns in the GRACE observations, such as the timing of melt and accumulation season commencement.

Discussion

Uniqueness of the solutions

Both our model input and the GRACE observations exhibit annual variations with similar phase. This may account for some proportion of the good performance of our optimized models. Although we have shown that our models are able to simulate high-frequency variations in mass (B ₁₀), these time series do have annual trends that could presumably be simulated with an optimized model driven by other annually varying time series. We investigate this issue by testing the extent to which weather station data used in this analysis can be used to predict glacier mass variations of other regions with well-defined annual cycles. In doing so, we can determine whether the good fit obtained in our analysis is unique to glaciers of the St Elias Mountains.

Model 2 was optimized against B10 calculated from each of 12 high-resolution mascon solutions for Gulf of Alaska glaciers (Fig. 1), all of which have distinct annual cycles of similar phase to those in the St Elias Mountains (Luthcke and others, 2008). We compared predicted and GRACE observations of B10 using data from the Yakutat weather station. The five smallest normalized RMSE and largest r ² values (indicating best model performance) were achieved with mascons 6, 7 and 10–12 (Fig. 4). Recall that mascons 6, 7 and 10 were used in the present analysis to represent glaciers of the St Elias Mountains.

Fig. 4. Coefficient of determination (r2, p < 0.001) and (b) RMS error in the fit of model 2 to each of the 12 GRACE mascon solutions for the Gulf of Alaska glaciers (Luthcke and others, 2008) (see Fig. 1). Errors in each model run are normalized by the variance of the associated GRACE mascon time series. The model is driven by temperature and precipitation data from Yakutat weather station. Solutions are plotted as a function of distance between the weather station and the centre of each mascon. Numbers beside points label each of the 1–12 mascons.

The results are plotted as a function of distance from Yakutat weather station. The error generally becomes larger and the r ² values become smaller with distance. This preliminary analysis suggests that data from Yakutat weather station could be used to model glaciers within approximately 500 km radius equally well, but that model performance would decrease outside of that radius. This is well within the length scale of mass-balance correlations (up to 1200 km) found in other studies (e.g. Cogley and Adams, 1998; Rasmussen, 2004).

Annual land mass variations

GRACE solutions for mountain glacier regions are complicated because the solution domain contains land cover in addition to glacier ice (Chen and others, 2006b, 2007; Luthcke and others, 2008). Whereas studies of multi-year glacier mass trends assume that all snow on land melts in a given season and does not affect the overall mass trend (Arendt and others, 2008), our present study of seasonal variations requires accounting for precipitation on both land and glacier surfaces. We have attempted to correct for terrestrial water storage using the GLDAS/Noah model described above, but snow precipitation data used to constrain GLDAS are sparse in Alaska (particularly in mountainous regions). In addition, our models assume that all meltwater immediately leaves the GRACE solution domain, when in reality some proportion of water will become stored within the glacier drainage system, raise lake levels or percolate into the land subsurface where it will remain stored for some time. Further studies will be necessary to quantify the magnitude of these errors and their effect on our modelling results.

Glacier dynamics

The flow of glaciers complicates our interpretation of the GRACE mascon time series because some mass may be transferred across mascon boundaries. Temperate glacier flow velocities have an annual periodicity that roughly matches annual climate variations, so that the magnitude of mass transfer probably varies on a seasonal basis. In addition, numerous glaciers in the St Elias Mountains are surge-type, several of which were actively surging during the GRACE measurements (Arendt and others, 2008). Surging results in a rapid displacement, but not removal, of glacier mass; these events will therefore not be captured by GRACE unless the surging glacier spans more than one mascon. Many other glaciers in the St Elias Mountains terminate in lakes or the ocean. Lacustrine and tidewater glacier dynamics can result in large mass changes that are not directly associated with trends in climate (Meier and Post, 1987), and may account for a large proportion of the mass loss in a given glacier region (Arendt and others, 2006; Larsen and others, 2007).

These dynamic changes have implications for the interpretation of our model parameters determined from simulations optimized to the GRACE time series. Because the GRACE signal records the effects of both climate and dynamics on glacier mass, and because some mass may flow between mascons, our optimized parameters will be different from those that would have been obtained from conventional surface mass-balance and local meteorological observations. In particular, the degree-day factors in model 1 and the coefficients (similar to mass-balance sensitivities) in model 2 should be interpreted with caution as they may not be suitable for use in conventional modelling studies of surface mass-balance.

Conclusions

Time variations in temperature and precipitation have been used in optimized models to simulate the annual and 10 day variations in the GRACE mascon time series for glaciers of the St Elias Mountains. Model optimization yielded parameters that generally agreed with values from previous studies. Modelled and GRACE observations of melt and accumulation seasonality agreed well despite the fact that the model optimization did not involve any phase adjustments. While this does not constitute an independent validation of the GRACE time series, it is encouraging that such simple models can capture the seasonal variability in glacier mass. We suggest that on a regional scale, local complexities resulting from dynamic changes, topographic shading or other factors average out to produce a mass-balance signal that is well captured by local climate station measurements. This has implications for regional degree-day modelling efforts, and suggests that simple models might perform well at capturing regional mass changes (provided the models are well calibrated).

We did not observe large differences between observed and modelled regional mass changes, indicating either that dynamic glacier losses do not make large contributions to the variability of glacier mass in this region, or that the timescales of these dynamic events have a seasonal component that matches the variability of the climate data. The highly dynamic nature of this region, as indicated by recent aircraft altimetry measurements, suggests the latter.

Our simplified 0-D model (model 2) produced simulations of the GRACE mascon time series almost identical to those of a model that distributed mass changes over all elevations within the region. We suggest that model 2 can be used to calibrate GRACE time series throughout the Gulf of Alaska region against existing climate station data. Investigations can then be carried out to determine whether these calibrated models can be used to extend the GRACE time series backwards or forwards in time. Existing aircraft altimetry measurements from previous studies can be used to further constrain and validate the model. Our model does not account for dynamic feedbacks between surface area, elevation and mass balance, or for the potentially large mass changes associated with tidewater glacier dynamics. We therefore caution against extrapolating model predictions in time before comparison with independent datasets.