Analysis of a GRACE global mascon solution for Gulf of Alaska glaciers

Abstract

We present a high-resolution Gravity Recovery and Climate Experiment (GRACE) mascon solution for Gulf of Alaska (GOA) glaciers and compare this with in situ glaciological, climate and other remote-sensing observations. Our GRACE solution yields a GOA glacier mass balance of −65 ± 11 Gt a⁻¹ for the period December 2003 to December 2010, with summer balances driving the interannual variability. Between October/November 2003 and October 2009 we obtain a mass balance of −61 ± 11 Gt a⁻¹ from GRACE, which compares well with −65 ± 12 Gt a⁻¹ from ICESat based on hypsometric extrapolation of glacier elevation changes. We find that mean summer (June–August) air temperatures derived from both ground and lower-troposphere temperature records were good predictors of GRACE-derived summer mass balances, capturing 59% and 72% of the summer balance variability respectively. Large mass losses during 2009 were likely due to low early melt season surface albedos, measured by the Moderate Resolution Imaging Spectroradiometer (MODIS) and likely associated with the 31 March 2009 eruption of Mount Redoubt, southwestern Alaska. GRACE data compared well with in situ measurements at Wolverine Glacier (maritime Alaska), but poorly with those at Gulkana Glacier (interior Alaska). We conclude that, although GOA mass estimates from GRACE are robust over the entire domain, further constraints on subregional and seasonal estimates are necessary to improve fidelity to ground observations.

Introduction

Runoff from glaciers surrounding the Gulf of Alaska (GOA) impacts rates of global sea level (Gardner and others, 2013) and crustal uplift (Fu and Freymueller, 2012), and is the primary source of fresh water to terrestrial and aquatic ecosystems in the region (Neal and others, 2010). Unlike other geodetic mapping methods that only detect multi-annual changes (e.g. Arendt and others, 2002; Berthier and others, 2010), data from the NASA/German Aerospace Research Center (DLR) Gravity Recovery and Climate Experiment (GRACE) have the potential to resolve temporal variations in glacier mass, providing a unique tool for both annual and sub-annual mass-balance investigations. The primary challenges in working with GRACE data are their coarse spatial resolution and the contamination of the glacier gravity signal from other geophysical sources. To date, GRACE has been used primarily to quantify the cumulative mass balance of GOA glaciers in order to determine the total contribution of this region to rising sea level (Tamisiea and others, 2005; Chen and others, 2006a; Luthcke and others, 2008; Pritchard and others, 2010; Wu and others, 2010; Jacob and others, 2012; Sasgen and others, 2012a). Most studies have focused on the long-term trend in a GRACE time series because it is much less affected by errors in the modeling of terrestrial water storage (TWS). Most of the variability in TWS in the GOA region occurs at sub-annual timescales that match variations in glacier mass, making it more difficult to discriminate between glaciological and other sources of mass variation.

Recent efforts have been made to compare GRACE GOA glacier mass-balance estimates with independent observations. GOA mass-balance estimates derived from a high-resolution mass concentration (mascon) approach (Luthcke and others, 2008) have previously been validated against observations from aircraft altimetry (Arendt and others, 2008) and weather station measurements of air temperature and precipitation (Arendt and others, 2009). Chen and others (2006b) compared a GRACE solution sampled at the location of two US Geological Survey (USGS) benchmark glaciers and found good agreement with field mass-balance observations. Other work (Hill and others, 2011) has pointed toward potential errors in existing GRACE solutions for GOA glaciers, especially at sub-GOA basin scales where mass-balance magnitudes exceed what is expected from geodetic data (Berthier and others, 2010). As studies attempt to extract further spatial resolution from GRACE solutions (e.g. Jacob and others, 2012) there is an even greater need to constrain and validate higher-resolution mass-balance estimates for the GOA and other glacier regions.

In this study we analyze the GOA glacier mascons from a new GRACE global mascon solution covering the period December 2003 to December 2010. This new solution is a subset of the high-resolution global mascon solution developed by the Space Geodesy Laboratory at the NASA Goddard Space Flight Center (Luthcke and others, 2013). We label this solution ‘v12’ following the version history of mascon solutions described in the literature (Luthcke and others, 2006a, 2008, 2013; Rowlands and others, 2010; Sabaka and others, 2010). Our goal is to compare temporal (annual to semi-annual) and spatial variations in the v12 solution with independent data sources to assess the accuracy of the local mass-balance signal. We compare the v12 solution with a series of satellite altimetric, in situ mass-balance and climate datasets that serve both to validate GRACE and to provide insight into the geophysical drivers of seasonal and interannual GOA glacier mass balances during the GRACE observation period.

Data and Methodology

Solution domain

We calculate gravity variations within a solution domain that encompasses most glaciers in Alaska, USA, as well as glaciers in Yukon and British Columbia, Canada, that are part of the icefields that straddle the US/Canada border (Fig. 1). On its southeastern boundary, our domain ends just south of the southern end of the Stikine Icefield. We exclude those glaciers of the Aleutian Island Arc (1793 km²) and the Brooks Range of northern Alaska (531 km²). The total glacier-covered area of our domain is 82 505 km², calculated from a new glacier inventory (Randolph Glacier Inventory version 2.0 or RGIv2.0) derived almost entirely from modern (~2005–11) satellite imagery (Arendt and others, 2012). We solve for GRACE mascon parameters at every 1 × 1 equal area mascon defined in Luthcke and others (2013). Our GOA region is a subset of the larger GOA domain defined in Luthcke and others (2013) that includes glaciers of the western Canadian Rockies. We exclude those glaciers in this analysis because they have traditionally not been included in earlier GOA mass-balance assessments. As a result, our mass-balance estimates are slightly different from those reported in Luthcke and others (2013). We combine mascons into eight mountain ranges, similar to those defined in Arendt and others (2002).

Fig. 1. GRACE solution domain for GOA glaciers. Boxes with gray outlines delineate 1 × 1 equal area degree mascons, and red outlines show groupings of these mascons into eight mountain ranges: 1. southwestern Alaska Range; 2. central Alaska Range; 3. eastern Alaska Range; 4. Chugach Mountains; 5. Wrangell Mountains; 6. St Elias Mountains; 7. Juneau Icefield; 8. Stikine Icefield. Blue shading showsthe location of Alaska glaciers, derived from a global glacier inventory (Arendt and others, 2012). Black dots show locations of ICESat elevation change measurements. Red dots show locations of USGS benchmark glaciers, and of Mount Redoubt volcano that erupted on 31 March 2009.

GRACE solution

The GRACE mission maps changes in Earth’s gravity field by precisely measuring changes in distance between tandem satellites using a microwave K-band inter-satellite range and range-rate (KBRR) system with GPS data for positioning, star tracker data for orientation and accelerometer data to remove the non-conservative surface forces needed to isolate the gravity signal. GRACE measures changes in gravity from all components of the Earth system. In order to isolate the glacier mass-balance signal, it is necessary to forward-model the time-varying gravity from atmosphere, ocean, TWS and solid Earth variations.

The v12 GRACE solution (Luthcke and others, 2013) calculates surface mass anomalies (mascons) directly from the GRACE KBRR observations taking into account the full noise covariance. We employ a unique method for accelerometer calibration (Luthcke and others, 2006b) and for calculating mass changes as scale factors on the differential set of Stokes coefficients in a spherical harmonic expansion of the geopotential field (Rowlands and others, 2005). The mascons are estimated globally at 1 × 1 arcdeg (~25 000 km²) spatial and 10 day temporal sampling, and we apply anisotropic temporal and spatial constraints between neighboring mascons representing similar surface types (land, ocean and glaciers). These neighbor constraints help to isolate signal and minimize signal leakage in and out of the land ice regions of interest. A detailed error and resolution analysis has shown that the basic mascon spatial resolution within a constraint region is equivalent to a Gaussian spatial smoother with 300 km radius (Luthcke and others, 2013). The signal within the GOA glacier region is isolated using these anisotropic spatial constraints, and the leakage errors estimated in Luthcke and others (2013) are included in our error estimates reported here. We estimate all errors at the 68% confidence (1σ) level.

The v12 mascon solution uses a similar GOA mascon solution domain and similar forward-modeling techniques to those in Luthcke and others (2008), including explicit accounting for post-Little Ice Age isostatic rebound resulting from the collapse of the Glacier Bay ice field (Larsen and others, 2005). A different ocean model is used, namely the Ocean Model for Circulation and Tides (Thomas, 2002), based on the Hamburg Ocean Primitive Equation Model (Drijfhout and others, 1996; Wolff and others, 1997). Over land we continue to use the 0.25° spatial, 3 hour temporal resolution GLDAS (Global Land Data Assimilation System) hydrology model (Rodell and others, 2004) with a 0.25° mask to remove glacierized regions from the TWS corrections, because these are not well modeled in GLDAS. We use a glacial isostatic adjustment (GIA) model based on the ICE-5G deglaciation history (Peltier, 2004), and an incompressible two-layer approximation to the VM2 viscosity profile (Paulson and others, 2007, computed and provided by Geruo A). We also apply geocenter corrections derived from the degree 1 Stokes coefficients determined from Swenson and others (2008). The v12 solution has been iterated three times to minimize residuals in modeled and observed KBRR data, resulting in improved recovery of non-modeled or mismodeled signal.

We present 10 day mass anomaly estimates together with a 10 day width Gaussian filter applied to the time series. We assume the 1σ errors for each 10 day estimate are the difference between the filtered and unfiltered estimates, which we combine in a root-sum-square fashion to estimate the total error of the GOA mass changes due to noise in the GRACE signal. These errors are combined with leakage and forward-modeling errors described in Luthcke and others (2013) to obtain the total GOA error reported here. We calculate seasonal balances as the difference between yearly maximum and minimum values in the Gaussian-filtered mascon time series, and we identify balance years based on the end-of-summer minima. For mass-balance calculations over the entire period of record we report both an average of the annual balances determined by summation of the summer and winter balances, as well as mass trends determined from a least-squares simultaneous estimate of bias, trend, annual, semi-annual and 161 day cycle (alias period of S2 tide errors).

We present summations of 1 × 1 mascon solutions over each of the eight mountain ranges, even though these mascon aggregates are not formally constrained in the way the total ice, land and ocean constraints are applied in the mascon approach (Luthcke and others, 2013). Our purpose is to ascertain the extent to which we can resolve mass variations at spatial scales smaller than the entire GOA region. Without any formal constraints, we are unable to quantify errors in the GRACE results we present for each of the eight regions. Instead, our comparison of the subregional GRACE estimates to independent climate and mass-balance observations is itself a form of error assessment, from which we can draw conclusions about the spatial resolution of the mascon approach.

ICESat data

We use Release 633 data from the GLA06 altimetry product of NASA’s Geoscience Laser Altimeter System (GLAS) on board the Ice, Cloud and land Elevation Satellite (ICESat; Zwally and others, 2002) to calculate changes in surface elevation of GOA glaciers, and compare these with our GRACE-derived mass balances. We choose ICESat data for this purpose because they cover the entire GOA during 2003–09, thereby fully overlapping with the GRACE observation period. Other elevation datasets derived from airborne altimetry (Johnson and others, 2013) or satellite stereoscopy (Berthier and others, 2010) are also available, but these span discontinuous time periods and do not cover the entire study region. The coverage of ICESat tracks over GOA glaciers is shown in Figure 1.

ICESat data were acquired over 17 repeated campaigns of a 33 day orbit sub-cycle between October 2003 and October 2009. ICESat provides surface elevations of ~70 m diameter laser altimetry footprints, each having a ~170 m along-track spacing. The elevation accuracy is highly dependent on surface slope but has been found to be within 1 m in typical glacier terrain (Moholdt and others, 2010). In addition, there is a recently discovered bias in ICESat elevations resulting from an incorrect selection of the centroid, rather than the Gaussian peak, of the transmit pulse to calculate surface elevation (the GC offset; personal communication from A. Borsa, 2013). We simulate the effect of this time-variable offset using global mean offsets calculated for each observation campaign.

We calculate elevation changes by examining elevation trends and residuals of surface planes fitted to 700 m long segments of near repeat-track data with 50% overlap, following methods in previous studies (Smith and others, 2009; Moholdt and others, 2010, 2012; Gardner and others, 2012, 2013). The method simultaneously solves for cross-track slope effects that are used to correct for repeat-track separation, as well as the true elevation changes used in our analysis. We use the same data-filtering parameters to remove data outliers as detailed in Moholdt and others (2010), except that we increase from 5 m to 10 m the maximum allowable elevation residual with respect to the trend of the plane. This was done to capture more of the steep topography and rapid elevation changes that are characteristic of many GOA glaciers.

A total of 80 000 ICESat laser returns passed our filtering criteria, resulting in 6257 elevation change estimates that each cover a time-span of minimum 2 years within the 2003–09 period. We examined the sampling density within each of the eight mountain ranges, but found that the coverage of ICESat tracks was not sufficient for reliable mass-balance estimation in the smallest GOA regions. However, at the scale of the entire GOA region, we found that ICESat sampling density as a function of elevation matched well with the glacier area–elevation distribution (Fig. 2b), allowing us to proceed with regional volume change calculations. For each 200 m elevation bin we calculated differences between the ICESat sampling and area distribution density, and used that as a weighting factor on the ICESat elevation changes. This corrected elevation change was multiplied by the total glacier surface area within each bin to yield the total volume change in the GOA region. We converted the volume changes to mass changes assuming Sorge’s law (Bader, 1954), and an average density of 900 kg m⁻³. We also calculated annual mass changes between ICESat campaigns in October/ November, which resulted in larger errors but allowed us to investigate links to our GRACE and climate observations.

Fig. 2. (a) Change in elevation determined from 2003–09 ICESat data, averaged within 200 m elevation bins. Error bars are the standard deviation of the mean change per elevation bin, divided by the square root of the number of uncorrelated observations. (b) Histograms of GOA glacier surface area (gray bars) and numbers of ICESat laser returns (orange line) used to estimate GOA glacier mass balance. The x-axis labels the midpoint of 200 m elevation bins over which the histograms and elevation changes are calculated.

We calculated errors by dividing the standard deviation of elevation change by the square root of the number of uncorrelated observations, assuming a correlation length of 5 km (Moholdt and others, 2012). We assume a fractional area uncertainty of ±10% (Arendt and others, 2006) and a density uncertainty of ±100 kg m⁻³ (Moholdt and others, 2012). The total mass-balance uncertainty was then calculated from the quadrature sum of these error components. Sensitivity tests revealed that our error budget was dominated by the large scatter in ICESat elevation changes, with uncertainties in density and area having relatively small effects on the total error. Also, the uncertainty of the annual mass balance for 2009 was higher than for the other years due to the failure of the last laser instrument one-third of the way into the October 2009 observation campaign.

Satellite snow-cover and albedo data

We use the Moderate Resolution Imaging Spectoradiometer (MODIS)/Terra Snow Cover Monthly L3 Global MOD10CM V005 dataset (Hall and others, 2011) to investigate the spatial extent of snow cover on glacier-free terrain during the GRACE observation period. The MODIS snow-cover data are distributed at 0.05° spatial resolution of the Climate Modeling Grid. Daily maps are determined from a normalized difference snow index and averaged for the monthly product. The daily MODIS snow map accuracy has been reported at 93%, with larger errors occurring over complex terrain, and where snow cover is thin (Hall and Riggs, 2007). The largest source of errors likely results from clouds and mixed pixels (Painter and others, 2009), which tend to confuse the algorithm, generally resulting in an underestimation of snow cover.

We sample the snow-cover grids at locations within our GOA solution domain exclusive of the ice-covered regions. Our purpose is to assess interannual patterns in snow accumulation that may not be captured in ground station observations or reanalysis products, and to quantify potential snow water equivalent errors in the GLDAS TWS model. We note that the 0.25° spatial, 3 hour temporal resolution version of the GLDAS model used here ingests the daily MODIS snow-cover product and adds an arbitrary snow water equivalent of 0.01 kg m⁻² to locations at which MODIS gridcells show >40% snow cover (Rodell and Houser, 2004). Here we assess the sensitivity of TWS calculations to the 40% threshold value by counting the number of MODIS gridcells classified above and below 40% snow cover.

We examine change in the spectral albedo of the glacier surfaces using the MODIS albedo product MCD43C4 for the years 2000–10. The MCD43C4 albedo is determined from the angular integration of a parametric model of the bidirectional reflectance distribution function (BRDF) that has been fit to 16 days of MODIS multi-angular reflectance observations in seven spectral bands. The product is available at 8 day intervals at a 0.05° horizontal resolution. We choose the coarse-resolution constant degree product over the higher-resolution 1000 m product (MCD43B3) for ease of processing, and because our GRACE analysis occurs on a low-resolution grid of mascons. The MODIS albedo product was found to agree with ground observations of snow albedo over the Greenland ice sheet within measurement error (Stroeve and others, 2005). However, particular attention must be given to the quality flags (QF) that accompany the dataset (Schaaf and others, 2011). Here we only use albedo estimated from fully inverted BRDFs (QF = 1). To ensure that we primarily capture changes in glacier albedo and not changes in terrestrial snow extent, we mask out all MODIS pixels with <95% glacier coverage as determined from the RGIv2.0.

Climate data

Here we build upon previous efforts to simulate seasonal variability of Alaska glaciers (Rasmussen and Conway, 2004; Arendt and others, 2009; Rasmussen and others, 2011), including an assessment of the relative influence of land surface physics on air temperature patterns, by comparing surface and 700 mbar air temperature records. For the upper air analysis we calculate the mean of two daily 700 mbar temperatures, one from the US National Centers for Environmental Prediction (NCEP) Reanalysis-2 and the second from the European Centre for Medium-Range Weather Forecasts (ECMWF) ERA-Interim reanalysis products. For each balance year in the GRACE record we calculate glacier area-weighted summer (June–August) departures from the 2004–10 mean for each of the eight mountain ranges. For the surface products we include both monthly temperature and precipitation fields derived from 1961–2009 station data using a total of 322 and 261 weather stations for temperature and precipitation respectively (note that 2010 data were not available). This surface product was calculated from absolute (for temperature) and proportional (for precipitation) anomalies of the station data from 1971–2000 climate normals (personal communication from D. Hill and S. Calos, 2011). These scattered anomalies were then interpolated onto a regular 2 km × 2 km grid covering all of Alaska. This approach assumes that the temporal derivative of the anomaly field is more spatially coherent than that of the normal field. The climate normals were taken from the Parameter-elevation Regressions on Independent Slopes Model (PRISM; Daly and others, 1994). PRISM uses a digital elevation model to correct for orographic effects on precipitation distribution, using regression equations that incorporate elevation and aspect for each gridcell. We assume precipitation is a suitable proxy for snowfall, given the difficulties in determining a solid precipitation temperature threshold when using monthly data. Summer average (June–August) temperature and winter (September–May) total precipitation were calculated to compare with summer and winter mass balances determined from GRACE. We calculated temperature and precipitation anomalies as the difference of these average/totalized values from their mean values over the period of record.

Field mass balance

We compare GRACE measurements to field observations of glacier mass balance acquired by the USGS at Gulkana and Wolverine Glaciers (Van Beusekom and others, 2010). Biannual observations are carried out at each of three ‘index sites’ that are used to represent different zones of each glacier. The glacier geometry and extent are updated every ~5 years. The glaciers are located in two distinctly different climate settings (Fig. 1). Gulkana Glacier had an area of 15 km² in 2009 and is located on the south flank of the eastern Alaska Range in a continental climate. Wolverine Glacier was 17 km² in 2009 and is located in the Kenai Mountains in a wet, maritime climate. Measurements on both glaciers have been extrapolated to produce estimates of total GOA mass losses (Meier, 1984; Meier and Dyurgerov, 2002) that are consistent with independent geodetic estimates (Cox and March, 2004; Harrison and others, 2009).

Results

GRACE-derived mass balance

Based on our least-squares estimate of the linear trend, the average annual mass balance (B _a) of GOA glaciers was −65 ± 11 Gt a⁻¹ during December 2003 to December 2010 (Fig. 3). When we average only the 2004–10 B _a values we obtain −71 ± 11 Gt a⁻¹. The difference results from the slightly different time period and the biases introduced in the linear fitting due to exceptionally negative summer balances at the start and end of the time series. We report both values because most GRACE studies calculate trends using the least-squares fitting method.

Fig. 3. Cumulative mass balance of GOA glaciers determined from the high-resolution GRACE v12 mascon solution (Luthcke and others, 2013). Ten-day estimates (blue dots), and the same data smoothed using a Gaussian filter with a 10 day averaging window (green curve); and trend (red line) recovered from simultaneous estimation of bias, trend, annual, semi-annual and 161 day cycle. GRACE errors for each 10 day estimate are assumed to be the difference between the filtered (green line) and unfiltered (blue dot) data.

The range in B _a in Alaska between 2004 and 2010 was 174 Gt, approximately double the mass loss rate estimated from a linear fit to the entire time series (Table 1). B _a became progressively less negative from 2004 to 2008, with near-balance conditions occurring in 2008. This was followed in 2009 by the most negative B _a in the 2004–10 period. We examined departures of summer (B _s) and winter (B _w) balances from the 2004–10 mean, calculated as the ratio of each year’s seasonal balance to the mean seasonal balance over the period of the GRACE record. Our sign convention is such that a negative departure in each season indicates a more negative mass balance. B _s departures (−83 to 91 Gt a⁻¹) were larger than B _w departures (−25 to 22 Gt a⁻¹), showing that variations in B _s accounted for most of the variability in B _a (Fig. 4). The slightly positive 2008 and large negative 2009 B _a values were caused by strong positive and negative departures of B _s in the respective years.

Table 1. Annual mass balances (B _a) of GOA glaciers between 2004 and 2010 balance years in units of mass (Gt a⁻¹), and in specific units (kg m⁻² a⁻¹) calculated by dividing by the total ice area in the mascon solution domain (82 505 km²). Balance years begin during the fall of the previous calendar year. B _a from GRACE is calculated as the difference between successive annual minima in Figure 3. B _a from ICESat is derived from average elevation changes between fall campaigns in October/November each year

Fig. 4. Departures (difference from 2004–10 mean) of GOA summer (orange) and winter (gray) mass balances. Seasonal balances are calculated as differences between successive minima and maxima in the GOA mascon time series (Fig. 3). Negative departures for both summer and winter balances indicate a contribution toward more negative annual balances, as compared to the 7 year mean.

We observed large spatial variability in the patterns of glacier mass loss between different balance years (Fig. 5). Notable mass gains occurred for glaciers in the central and eastern Alaska ranges (mountain ranges 2 and 3) during the 2008 balance year, the Stikine Icefield (mountain range 8) during 2007, and glaciers of the Alaska and Kenai peninsulas (mountain ranges 1 and 4) during 2008 (Fig. 1). Large mass losses occurred through all regions during the 2004 and 2009 balance years due to exceptionally negative B _s. During the entire period, the St Elias, Glacier Bay and Juneau icefields regions lost the greatest amount of mass relative to other regions. The amplitudes in glacier mass increased with decreased latitude, and were largest in southeast Alaska.

Fig. 5. Spatial distribution of summer, winter and annual mass balance (orange, gray bars and red lines) in Alaska from the GRACE v12 mascon solution (Luthcke and others, 2013), for each balance year spanning 2004–10. Mountain ranges (numbered 1–8) are labeled in Figure 1. Summer, winter and annual balances are calculated as in Table 1. Blue shading shows glacier locations.

Comparison to ICESat-derived mass balance

ICESat elevation changes averaged over the 2003–09 observation period show a strong elevation dependence, with 2–4 m a⁻¹ of thinning at low elevations tapering to near-zero changes at high elevations (Fig. 2a). We multiplied these elevation changes by the area distribution of GOA glaciers to obtain an average B _a of −65 ± 12 Gt a⁻¹ or −0.79 ± 0.15 kg m⁻² a⁻¹ (Table 1). This is 4 Gt a⁻¹ more negative than GRACE data subsampled to the same period as ICESat. GRACE and ICESat also agree to within error bars for individual balance years (B _a) except for 2007 (Table 1; Fig. 6). Our simulations of the ICESat GC offset indicated that the total impact on the GOA mass balance was <1 Gt a⁻¹. Because this was well within our calculated ICESat error estimates, we do not include this correction in our analysis.

Fig. 6. Comparison of area-averaged annual mass balances in the GOA region derived from GRACE (orange) and ICESat (gray). Balance years begin in the fall of the previous calendar year, at the time of minimum mass (GRACE), or in October/November (ICESat).

Comparison to field observations

We compare field observations of B _w, B _s and B _a at Gulkana and Wolverine Glaciers to their closest 1 × 1 mascon (Fig. 7). We convert the GRACE mass changes to specific values by dividing by the total ice area in each mascon. Our scatter plots illustrate both the strength of the correlation between the two datasets (i.e. similarity of patterns from one year to the next) as well as the level of agreement in magnitude between GRACE and ground observations (measured by the departure of points from the diagonal one-to-one line). At Gulkana Glacier, B _a from GRACE and ground observations were poorly correlated (r ² = 0.05, p = 0.1), and the absolute values of GRACE B_s and B_w values were larger than ground observations. For Wolverine Glacier, B _a had a relatively high correlation (r ² = 0.75, p = 0.79) between ground and GRACE observations, and the magnitudes of field-derived B _s and B _w were very close to those observed by GRACE. We note that the higher correlation between GRACE- and field-derived B _a values for Wolverine Glacier may have occurred due to compensating errors in the associated B _s and B _w values, as indicated by their lower correlation to the GRACE data.

Fig. 7. Comparison of area-averaged summer (orange squares), winter (gray circles) and annual (red diamonds) mass balances from GRACE and those derived from USGS field observations at (a) Gulkana Glacier and (b) Wolverine Glacier between 2003 and 2009. GRACE data are converted to specific values by dividing the mass change by the total ice area in mascons containing each glacier.

Comparisons to climate observations

Both the 700 mbar reanalysis air temperatures and the interpolated ground surface temperatures captured a large amount of the variability in GRACE-derived B _s data (r ² = 0.72, p = 0.02 and r ² = 0.59, p = 0.07 respectively; Fig. 8; note the inverted secondary y-axis). In nearly all cases, surface and 700 mbar anomalies were positive during years of negative B _s. A notable exception was 2009, a year of strongly negative B _s during which 700 mbar temperatures were slightly above normal, but surface temperatures were slightly below normal. This was in contrast to 2004, when a similarly large negative B _s value was matched by strongly positive surface and 700 mbar anomalies. We note that during 2004 the magnitude of the 700 mbar temperature anomaly was approximately four times larger than in 2009, even though B _s anomalies were similar between the two years. During 2004, 700 mbar air temperature departures from the 2004–10 mean were distributed across the GOA region, while in 2009 the highest departures occurred over the St Elias, Juneau and Stikine icefield regions (Fig. 9). 2008 had the largest positive summer balance departure that was matched by the most negative surface and 700 mbar air temperature departures. The 2008 temperature departures were distributed relatively uniformly across the region (Fig. 9).

Fig. 8. Departures of GRACE GOA summer glacier mass balances from the 2004-10 mean (orange circles, left axis; note that negative summer balance departures imply greater summer mass loss). GOA summer (June–August) surface (black squares) and 700 mbar lower-troposphere (black circles) air temperature departures from the 2004–09 mean (right axis; note inverted scale).

Fig. 9. Departures in summer (June–August) 700 mbar air temperature (mean of NCEP R-2 and ERA-Interim reanalysis) from 2004–10 mean at each of the eight GRACE mascon mountain ranges.

There was no significant correlation between GRACE-derived B _w and winter (September–May) precipitation determined from PRISM (r ² = 0.16, p = 0.43) and from the reanalysis products (r ² = 0.07, p = 0.57; Fig. 10). The winter balance years 2006 and 2008 had above-normal winter balances together with above-normal surface snow accumulation values, while the opposite occurred for the 2007 and 2009 balance years.

Fig. 10. Departures of GRACE GOA winter mass balances from the 2004–10 mean (gray circles, left axis). GOA winter (September–May) total winter snowfall departure from the 2004–09 mean (right axis). Data are derived from station observations (black squares) and from reanalysis model output (red squares).

Spectral albedo values averaged over May–July for GOA glacier surfaces showed a 50% increase in shortwave absorption during 2009 relative to the average of the remaining years in 2004–10 (Fig. 11). The greatest decrease in albedo occurred in the visible and near-infrared spectrum (0.40–0.9 μm).

Fig. 11. Average May–July MCD43C3 MODIS-derived spectral albedo for all gridcells within a 250 km radius of Mount Redoubt and having a glacier coverage of 95% or greater for the years 2000–10. Areas listed next to sample years indicate the total gridcell area sampled each year.

Snow cover on glacier-free terrain (an annual average of all snow-cover values in the GOA region exclusive of glaciers) had maximum values of 90–96% and minimum values of 7–13% within the 2004–10 time period (Fig. 12). The summer of 2008 had the largest minimum snow cover (13%), corresponding with the least negative B _s discussed above. For most years the average snow cover was <40%, and therefore below the threshold for inclusion in the GLDAS model, during May–July. To assess potential errors in eliminating cells with <40% snow cover, we multiplied the number of cells with <40% snow cover by 0.01 kg m⁻² (Rodell and Houser, 2004), yielding a value of 2.5 Gt. We lack further information to more accurately assess the water equivalent depth of snow remaining, but we provide this estimate to illustrate the potential order-of-magnitude impact of late-season snowpacks on our mass trend retrievals.

Fig. 12. Average percent snow cover on glacier-free terrain in the mascon solution domain determined from the MODIS Monthly L3 Global snow-cover dataset. Red line is the monthly time series, and gray bars mark the magnitude of the summer minimum value.

Discussion

Our mass-balance estimates averaged over the 2004–09 balance years from GRACE and ICESat compare well, providing mutual validation of the magnitude of GOA glacier contributions to rising sea level during this period. Our estimates are more negative than two other GRACE estimates for GOA glaciers over the same time-span (−42 ± 6 Gt a⁻¹ (Jacob and others, 2012) and −54 ± 13 Gt a⁻¹ (Sasgen and others, 2012b)). The relatively strong agreement between GRACE- and ICESat-derived B _a for individual years provides good support for our assessment that interannual variations in GOA GRACE signals are due to fluctuations in glacier mass. The elevation dependence of ICESat elevation changes also matches recent airborne altimetry observations (e.g. Johnson and others, 2013), lending further support to the accuracy of our ICESat analysis for GOA glaciers.

Both surface and 700 mbar air temperatures explain much of the variability in B _s as observed by GRACE. The slightly stronger correlation for 700 mbar suggests lower-troposphere data may be more suitable than surface data for analysis of snow and ice variations, supporting previous studies (Rasmussen and Conway, 2004). This is likely because 700 mbar temperatures are less subject to localized surface effects and are more indicative of regional temperature variations. The relatively weak correlation of B _w with snow accumulation observations likely results from difficulties in quantifying snow cover from precipitation data, and the paucity of ground observations to constrain the reanalysis and surface precipitation products.

We attribute the strongly negative B _s in 2009 to the 31 March 2009 eruption of Mount Redoubt that spread ash over a large area, causing a significant reduction in surface albedo as measured by MODIS. Ash dispersal maps from the eruption show the plume progressing primarily north and east of the eruption center, with ash fall covering an estimated 80 000 km² in the vicinity of Mount Redoubt (Wallace and others, 2013). Trace ash (particle size <0.8 mm) was observed as far away as Fairbanks, Alaska, located 550 km north-northeast of the volcano (Wallace and others, 2013). These observations are consistent with our 2009 field operations, during which several of us observed ash as far away as the central Alaska and St Elias Ranges (~300 and 650 km from Mount Redoubt, respectively). Snow albedo is very sensitive to the presence of highly absorbing atmospheric aerosols such as those expelled during a volcanic event (Conway and others, 1996). We suggest the Redoubt eruption caused an enhancement of surface melt through greater absorption of solar radiation, providing a mechanism for increased GOA glacier mass loss even in the absence of higher air temperatures.

The sequence and timing of air temperature and snow accumulation anomalies appear to be important in controlling GOA mass-balance patterns. The relatively large winter accumulation during the 2008 balance year was followed by a cool summer that allowed snow to remain unmelted on glacier-free surfaces. The fact that most of the GOA mascons had <40% snow-cover values during the summer months means that they were not accounted for in the GLDAS model (Rodell and Houser, 2004). In the event that summer snowpacks registering <40% cover in the MODIS product are of sufficient extent or thickness to be sensed by the GRACE satellites, this is a bias in our calculations that should be addressed in future work.

B _s and B _w values predicted by GRACE are much larger than is physically realistic for Gulkana Glacier, and far exceed any in situ seasonal balances measured during its entire record dating back to 1966. These erroneous GRACE amplitudes likely result from some combination of signal leakage across the solution domain (Sabaka and others, 2010) and signal contamination due to TWS variations on glacier-free surfaces. Concerning signal leakage, we note that although our GRACE solution does have spatial and temporal constraints applied to the entire GOA solution domain, no constraints have been applied to either the individual mascons or mountain ranges. These constraints help to minimize leakage of signal between different broad-scale geophysical systems, such as land, glacier and ocean surface types. Without constraints, individual mascons can contain signal that is to some extent smeared across a larger region of the Earth’s surface. For example, Luthcke and others (2013) found that signal from a single mascon can affect other mascons at spatial scales up to 600 km. These problems reflect fundamental limitations of the GRACE satellites that are made worse when applying GRACE to a complex region such as the GOA that has a highly nonuniform distribution of glaciers near a highly dynamic ocean boundary. The fact that our comparisons with Gulkana Glacier, located in a region with sparse glacierization, were much worse than with Wolverine Glacier, located in a heavily glacierized region, suggests that signals from the large ice masses in the St Elias and Chugach Mountains are leaking into mascons in the Alaska Range and other regions.

TWS modeling errors could also account for discrepancies between ground observations and GRACE. Luthcke and others (2008) showed that forward modeling of TWS using the same dataset employed for our v12 solutions (the GLDAS/NOAH product) did not match GRACE observations over much of the glacier-free surfaces of interior Alaska during 2003–08. TWS errors result from a paucity of snow-depth, groundwater and other observations needed to calibrate the GLDAS product. If TWS is mismodeled over glacier-free surfaces located within glacier mascons, this will result in erroneous attribution of terrestrial snow cover, groundwater and other TWS signals to glacier surfaces. Such errors have the potential to be greatest for glacier mascons with low concentrations of glacier ice, wherein a greater proportion of TWS corrections are being applied. We note that for those years when all snow cover melts on glacier-free terrain, these errors only affect the sub-annual components of mass change, and not the annual trend.

Conclusions

GOA glaciers are controlled by numerous environmental forcings that complicate interpretation and modeling of their mass changes, especially over short time periods. During 2003–08, the mass balance of GOA glaciers was strongly correlated to mean summer air temperatures, emphasizing the importance of acquiring air temperature data from high mountain regions in order to improve model predictions of Alaska glacier surface mass balances. At the same time, our 2009 observations highlight potential limitations in models that only include temperature as a driver of glacier ablation. Due to the regional reduction in albedo from the 31 March 2009 eruption of Mount Redoubt, temperature alone was an insufficient proxy for glacier surface mass balance, and other factors (e.g. changes to the optical properties of the snow and ice surface) may have been the dominant controls.

Studies aimed at matching GRACE observations to field observations require consideration of leakage errors and fundamental resolution constraints inherent in satellite gravimetry datasets. We have shown that individual mascon time series combined additively for the purpose of hydrological analysis within mountain ranges generally do not match the magnitude of ground observations, even though they are highly temporally correlated. We attribute this to the smearing of signal within the GOA domain and the mismodeling of TWS variations, both of which are enhanced by the complex spatial distribution of glaciers within the region and the lack of sub-GOA solution constraints.

Snow on ground likely plays an important role in controlling the mass balance of individual glacier catchments, and biases to GRACE-derived mass balances may occur when the summer snowpack on glacier-free terrain does not completely melt in a given year. Further improvements to TWS models and to the acquisition of accurate snow water equivalence data are important to improve recovery and attribution of GRACE mass-balance signals in this region.

We have provided the first independent validation of GRACE glacier mass-balance estimates spanning the entire GOA region, using data from the ICESat mission. The agreement of these two estimates strongly supports our assertion that inaccuracies in sub-annual mass balances, due to problems outlined above, are sufficiently random over the length of the GRACE record that they do not bias our multi-annual estimates. This suggests that our existing GRACE processing constraints used to separate mass change signals occurring within the entire GOA region from surrounding land and ocean surfaces are robust. The relatively strong correlation between B _a values derived from GRACE and ICESat indicates that future satellite altimetric missions, which will have an even higher sampling density than ICESat, will be a valuable tool for annual GOA mass-balance assessments.