Grace mascon solution

Estimates of the mass variations of glaciers surrounding the Gulf of Alaska have been obtained from high-resolution GRACE mascon solutions (Reference Luthcke, Arendt, Rowlands, McCarthy and LarsenLuthcke 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 (Reference EkEk and others, 2003; Reference RodellRodell 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 (Reference Raup, Kieffer, Hare and KargelRaup 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 (Reference RowlandsRowlands and others, 2005; Reference Luthcke, Rowlands, F.G., Klosko, Chinn and McCarthyLuthcke 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.

Reference Arendt, Luthcke, Larsen, Abdalati, Krabill and BeedleArendt 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 Reference Arendt, Luthcke, Larsen, Abdalati, Krabill and BeedleArendt 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 Reference Arendt, Luthcke, Larsen, Abdalati, Krabill and BeedleArendt 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. Reference Swenson and WahrSwenson 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 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. Reference Eisen, Harrison and RaymondEisen and others, 2001; Reference RasmussenRasmussen, 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. Reference RasmussenRasmussen, 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. Reference Jóhannesson, Laumann and KennettJó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 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 (Reference De Woul and HockDe 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.

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 (Reference Luthcke, Arendt, Rowlands, McCarthy and LarsenLuthcke 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. Reference Cogley and AdamsCogley and Adams, 1998; Reference RasmussenRasmussen, 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 (Reference Chen, Wilson and TapleyChen and others, 2006b, Reference Chen, Wilson, Tapley, Blankenship and Ivins2007; Reference Luthcke, Arendt, Rowlands, McCarthy and LarsenLuthcke 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 (Reference Arendt, Luthcke, Larsen, Abdalati, Krabill and BeedleArendt 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 (Reference Arendt, Luthcke, Larsen, Abdalati, Krabill and BeedleArendt 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 (Reference Meier and PostMeier and Post, 1987), and may account for a large proportion of the mass loss in a given glacier region (Reference ArendtArendt and others, 2006; Reference Larsen, Motyka, Arendt, Echelmeyer and GeisslerLarsen 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.