Geometry, mass balance and thinning at Eklutna Glacier, Alaska: an altitude-mass-balance feedback with implications for water resources

Abstract

We analyzed glacier surface elevations (1957, 2010 and 2015) and surface mass-balance measurements (2008–2015) on the 30 km² Eklutna Glacier, in the Chugach Mountains of southcentral Alaska. The geodetic mass balances from 1957 to 2010 and 2010 to 2015 are −0.52 ± 0.46 and −0.74 ± 0.10 m w.e. a⁻¹, respectively. The glaciological mass balance of −0.73 m w.e. a⁻¹ from 2010 to 2015 is indistinguishable from the geodetic value. Even after accounting for loss of firn in the accumulation zone, we found most of the mass loss over both time periods was from a broad, low-slope basin that includes much of the accumulation zone of the main branch. Ice-equivalent surface elevation changes in the basin were −1.0 ± 0.8 m a⁻¹ from 1957 to 2010, and −0.6 ± 0.1 m a⁻¹ from 2010 to 2015, shifting the glacier hypsometry downward and resulting in more negative mass balances: an altitude-mass-balance feedback. Net mass loss from Eklutna Glacier accounts for 7 ± 1% of the average inflow to Eklutna Reservoir, which is entirely used for water and power by Anchorage, Alaska's largest city. If the altitude-mass-balance feedback continues, this ‘deglaciation discharge dividend’ is likely to increase over the short-term before it eventually decreases due to diminishing glacier area.

1. INTRODUCTION

The changes in glacier hypsometry caused by negative mass balances can conceptually be separated into terminus retreat (change in area) and thinning (change in surface elevations). The relative proportions of retreat to thinning determine the feedback between the change in hypsometry and future glacier mass balance. The glacier response is ‘stable’ if, in the absence of further climatic change, the annual mass balances trend toward zero and the glacier geometry approaches a new steady state (Harrison and others, 2009). Many mountain glaciers, particularly if they are steep, respond to negative mass balances primarily through terminus retreat, which reduces ablation losses in a negative feedback that stabilizes the now smaller glacier toward neutral mass balances (Harrison and others, 2001; Oerlemans, 2008). Conversely, thinning shifts the hypsometry to lower elevations and increases ablation, causing the mass balances to become more negative (Harrison and others, 2001; Paul, 2010; Huss and others, 2012). This response, termed altitude-mass-balance feedback, can be ‘unstable’ where the positive feedback between surface altitude and mass-balance outweighs terminus retreat and amplifies future mass losses (Böðvarsson, 1955; Trüssel and others, 2013).

The nature of a glacier's geometric response to negative balances has important hydrologic consequences. The loss of glacier mass during warm years with negative balances leads directly to enhance downstream flows: the ‘deglaciation discharge dividend’ (Collins, 2008). The relative importance of terminus retreat versus thinning affects the size and persistence of that deglaciation dividend (Huss and others, 2012), and will ultimately alter the seasonality of downstream flows (Barnett and others, 2005). If the mass balance returns to equilibrium, that deglaciation dividend is lost, and downstream flows may be further diminished by water loss from interception and evaporation on newly deglaciated terrain (Collins, 2008). Predicting the likely timing and magnitude of these changes is critical in situations where runoff provides a valuable water resource.

Runoff from Eklutna Glacier, in Southcentral Alaska, is used for both potable water and hydroelectric power production for Anchorage, Alaska's largest city. Eklutna Glacier includes 74% of the 39 km² total glacierized area within the 307 km² Eklutna basin, yet 45–50% of the total reservoir inflow is runoff from the 64 km² area sub-basin that includes Eklutna Glacier (Larquier, 2011). Larquier (2011) estimated that the larger (101 km²) sub-basin to the east (including 19 small glaciers with a total area of 10 km²) accounts for <40% of the inflow to the reservoir, and that mass loss from Eklutna Glacier proper accounted for 24 and 3% of the reservoir budget during 2009 and 2010, respectively. Despite the relevance of this basin to the Anchorage water supply, the only prior glaciological research beyond casual observations of terminus retreat consists of stream gauge measurements and short duration ablation measurements (partial ablation seasons; 1985–1988; Brabets, 1993).

In this paper, we investigate the ongoing response of Eklutna Glacier to changing climate using surface mass-balance and surface elevation observations. We use those to calculate both glaciological and geodetic mass balances, which we compare. We evaluate how changes in the amount and density of firn, and the timing of our geodetic observations, impact the observed changes in surface elevations. We then calculate ice flux from the mass balance and thinning. Concurrent with mostly negative annual balances, Eklutna Glacier has experienced substantial thinning of a broad high-elevation basin that includes much of the accumulation zone. This is particularly important because much of the glacier surface area is near the present-day equilibrium-line altitude (ELA). In the discussion, we consider the rates and distribution of thinning, and use the ice-flux calculations to demonstrate the implications of those changes for future glacier hypsometry and consequently future mass balances and water supplies.

2. SITE DESCRIPTION

Eklutna Glacier is located in Southcentral Alaska's Chugach Mountains, 50 km northwest of Prince William Sound (Fig. 1). It is 10.2 km long, 29.5 km² in area, and ranges in elevation from 580 to 2100 m (2010 statistics; Fig. 1c). Two tributaries converge 2.7 km above the terminus, but only 6% of the glacier area is located below the convergence. We refer to the longer, larger tributary (56% of the total glacier area) as the main branch, and the shorter, smaller tributary (44% of the area) as the west branch. Below the confluence, a medial moraine delineates the branches, showing that the west branch contribution ends 0.5 km above the terminus. Within the main branch, much of the glacier area is contained in a broad, low-slope ‘upper basin’. This upper basin constitutes 60% of the main branch area yet occupies only 10% of the glacier's total elevation range. Importantly <10% of the branch area is on the steeper flanks (15–45°) of the peaks above this very flat basin. Below a lateral constriction at ~1360 m, the lower glacier is progressively steeper and narrower toward the terminus (Fig. 1c). The west branch, in contrast with the main branch, is smaller, narrower, much steeper in its upper reaches, and extends to higher elevations. The west branch hypsometry is distributed over a broader range of elevations and shows two distinct peaks at 1400 and 1600 m (Fig. 2). We will show that these fundamental differences in hypsometry, slope, aspect and shading give rise to different patterns of both mass balance and thinning in each branch, and we treat them separately for most of our analysis.

Fig. 1. (a) Location of Eklutna Glacier in Southcentral Alaska. (b) Detail of box in (a) showing Eklutna Glacier (1, red), Eklutna Reservoir (2, blue), the watershed (outlined in black), Anchorage (3), and Wolverine Glacier (4, green). (c) Eklutna Glacier topography, 2010, showing the 1957, 2010 and 2015 extents. We refer to the area above 1360 m in the main branch (red line) as the upper basin.

Fig. 2. Eklutna Glacier mass-balance observations (circles) and profiles (lines), 2008–2015. Both panels show the 2010 glacier surface hypsometry in gray with values on the upper horizontal axis. Circle diameters reflect ± 0.2 m estimated measurement uncertainty.

Eklutna Glacier is in a transitional zone between maritime and continental climates (Bieniek and others, 2012). From 2008 to 2015 we measured 4–6 m snowpacks in the accumulation area (above ~1450 m elevation) by the end of winter and high (1–3 m w.e. a⁻¹) mass-balance amplitudes (Meier, 1984). Summer mass-balance rates at the terminus ranged from −8 to −12 m a⁻¹. Air temperatures near the ELA vary by more than 40°C, ranging from −30°C to +10°C. A majority (50–60%) of precipitation usually falls in the four months from July to October (Bieniek and others, 2012).

3. DATA

In this paper, we rely on two main types of measurements: surface mass-balance and glacier surface elevations. Additionally, we use weather data from an on-glacier automatic weather station (AWS), and ice velocities from the ablation stakes.

3.1. Glaciological mass balance

Direct measurements of surface mass balance were made seasonally between 2008 and 2015 at 3–8 sites. At each site (Fig. 1c) we placed an ablation stake, which we maintained over annual time intervals by periodically shortening or extending the stake. Site visits were timed to be close to the beginning and end of the melt season, so that ablation measurements include most of the seasonal amplitude. The observed change in snow and or ice thickness was converted to water equivalence based on density observations. Each spring, accumulation was calculated from depth and density measurements collected in snow pits and cores at these sites.

Measurement sites evolved over the course of this project. Three primary sites were established near the centerline of the main branch in 2008, but were subsequently moved up-glacier in 2009 to better distribute them over the glacier hypsometry (Fig. 1c). Additional sites were added throughout the program as logistics allowed (data at doi.org/10.5066/F7MP51CB/). Stake and pit measurements were complemented by late-summer observations of ELA. No attempt was made to estimate internal or basal processes; hereafter ‘mass balance’ refers only to the ‘surface mass balance’ (Cogley and others, 2011).

During the melt season from 2008 to 2015, we maintained an on-glacier AWS at the main branch measurement site near the long-term ELA (~1400 m; Fig. 1c). This station was composed of a floating pyramid with sensors for air temperature, relative humidity, wind speed and direction, and incoming and reflected solar radiation measured at 2 m above the surface. A sonic ranger mounted on a drilled ablation stake measured surface lowering. These hourly data were used to estimate how the measurements of net ablation were distributed through the melt season.

We used a mapping grade (L1) GPS system to collect and post-process stake positions from which we calculated seasonal and annual surface velocities. Baselines were <40 km, and the resulting uncertainties were <10% of the total velocities (Table S3).

3.2. Geodetic data

We acquired seven glacier surface elevation and boundary datasets, three of which are presented here and used in our primary analyses. The other four complement the primary acquisitions and support the same conclusions, but analyses are complicated by the limited spatial coverage and shorter time intervals, so we present them in the supplemental material for simplicity. The earliest glacier geometry available is a DEM produced from 12 July 1957 stereo imagery published in the US Geological Survey National Elevation Dataset (NED; Gesch and others, 2002). The DEM has a 2 arc-second resolution (~30 × 60 m² post-spacing). It was produced from 1: 63 360 maps with a 100 ft (30.48 m) contour interval. The glacier boundary was copied directly from the map used to create the NED.

The US Geological Survey contracted full-coverage airborne lidar acquisitions over several areas of Alaska in 2010, including Eklutna Glacier, using an Optech Gemini airborne scanning system that produced a point cloud with nominal post-spacing of 1.9 m. The US Geological Survey EROS Data Center filtered point-cloud outliers and gridded the data at 2.5 m (http://lidar.cr.usgs.gov/). We digitized the 2010 glacier boundary from a hillshade rendering of this DEM. These data were collected on 16 September, 5 d prior to the fall mass-balance measurements in 2010.

We also acquired Worldview 3 stereo pairs on 24 August 2015 (under the Nextview license, Digital Globe). Surface elevations and the glacier boundary were derived with photogrammetric methods using the Ames Stereo Pipeline (Shean and others, 2016) and gridded at 2 m. We digitized the 2015 glacier boundary directly from the orthorectified visible image. These images were collected 27 d prior to the fall mass-balance measurements in 2015.

4. ANALYTICAL METHODS

We utilized the data described above to create time series of mass-balance and surface elevation changes as functions of elevation, and used it to assess the recent geometric evolution of Eklutna Glacier, the relationship of that evolution to mass balance, and the implications for water resources. Analyses were performed separately for each branch due to differences described previously.

4.1. Glaciological mass balance

Using observations from stake data, we calculated glaciological mass balances as functions of elevation (balance profiles) for each year on each branch of the glacier (Fig. 2) using the field survey dates (i.e. a floating date system; Cogley and others, 2011). Observations within 50 m vertically were binned together, and we interpolated a profile using a ‘pchip’ function (Fritsch and Carlson, 1980), which rounds discontinuities in the profile (piecewise linear interpolation yielded glacier-wide results within ± 0.01 m w.e. a⁻¹). Outside the elevation range of our observations, we apply the mean profile gradients from all years in the accumulation zone and ablation zone to extrapolate the balance profile above and below the observations, respectively. We utilized the measurements from our lowest site, near the confluence of the two tributaries, for both branches. Prior to commencement of direct measurements in the west branch in summer 2011, we used the ELA position (estimated from late season satellite imagery) as an upper elevation balance measurement. From 2011 forward, we combined the observed west branch ELA position with available direct measurements to fit a linear profile.

Our results are ‘conventional balances’ in the sense of Elsberg and others (2001), meaning that glacier-wide mass balances were calculated over a time-variable hypsometry. The hypsometry was linearly interpolated from 2010 to 2015, and that rate of change was extrapolated to 2008 and 2009.

4.2. Change in glacier thickness

In order to determine the spatial distribution of glacier mass change, the annual ice-equivalent thinning rate (the deficit between the mass balance and ice flow) is required. We started by differencing measured surface elevations from 1957, 2010 and 2015, assuming no uplift, basal erosion or subglacial sediment deposition. Over both intervals (1957–2010 and 2010–2015), we calculated raw surface elevation changes first by co-registering each pair of DEMs (using ‘universal coregistration’; Nuth and Kääb, 2011), and then differencing a bilinear interpolation of the finer DEM from the post-locations of the coarser DEM.

The derived changes in surface elevation are not a direct measurement of ice-equivalent thinning. The interannual variability in mass-balance drives changes in firn extent, thickness, and density, which can be large even over short timescales (Huss, 2013). Elevation data acquired before the end-of-season mass-balance observations do not account for the ablation or the ice flow between the two dates, both of which can be large in a high mass turnover environment.

From 1957 to 2010, we had minimal data to model changes in the firn distribution and density, or late season rates of surface melt and ice flow – processes that are small compared with the long time interval and large thinning signal. Without a firn model, Huss (2013) recommends accounting for a likely loss of firn by assuming a glacier-wide mean value of 850 kg m⁻³ for the observed volume loss. That recommendation was based on a distributed model where densities were elevation dependent, and hence taken to be 900 kg m⁻³ in the ablation zone and <850 kg m⁻³ in the accumulation zone. We followed suit, assuming a density of 900 kg m⁻³ for material lost in the ablation zone (below the observed snowline at 1360 m in the 1957 imagery) and calculating the effective density of material in the accumulation zone required to average a glacier-wide value of 850 kg m⁻³. That effective density works out to 820 kg m⁻³ for the accumulation zone, which yields a more conservative thinning rate within the upper basin.

Over the shorter (2010–2015) interval, we used surface mass-balance observations to model the changes in firn density from year to year, and to estimate the portion of the measured ablation and ice flow after the surface elevation acquisitions. For a glacier-wide geodetic mass balance, model results suggest that the constant density assumption would be adequate for the 5-year interval from 2010 to 2015 (Huss, 2013), but we are interested specifically in the thinning within the accumulation zone, which we would expect to be more sensitive to interannual changes in firn density. To determine firn thicknesses and densities at the times of the 2010 and 2015 DEM acquisitions, we used a simple elevation dependent model of compaction and densification based on the surface mass-balance profiles (Huss, 2013).

For mass balances prior to 2008 (for model spin-up purposes only) we used the mean measured balance profile at Eklutna Glacier adjusted by the annual deviation from the mean profile as estimated by the Wolverine Glacier mass-balance time series (Van Beusekom and others, 2010; O'Neel and others, 2014). Wolverine is smaller (16.7 km²) than Eklutna and ~80 km south but in a similar climate (Fig. 1b). Eklutna mass balances show lower interannual variability (Fig. 3b), but a similar pattern and cumulative mass loss (r ² = 0.81, p = 0.005), implying that Eklutna Glacier has had a similar relationship to regional climate in recent decades.

Fig. 3. Glacier-wide glaciological mass balances 2008–2015. (a) Eklutna Glacier annual balances, with the same colors as Figure (2). (b) Comparison with Wolverine Glacier annual mass balances.

To adjust 2010 and 2015 surface elevation measurements to the end of the melt season when surface mass-balance measurements were made, we calculated late season ablation by partitioning the total summer ablation into two periods based on the positive degree days before and after the surface elevation observations. Then, we used mass continuity to estimate horizontal flux divergence rates, which we used to calculate the late season vertical component of ice flow.

Mass continuity relates changes in glacier thickness, ∂h/∂t, to the mass-balance rate, $\dot b$ , through the horizontal ice-flux divergence, $\nabla _{xy} \cdot \vec q$ , of the glacier flow,

(1) $$\displaystyle{{\partial h} \over {\partial t}} = \dot b - (\nabla _{xy} \cdot \vec q)$$

(Cuffey and Paterson, 2010). Equation (1) assumes that the density is constant through time (i.e. Sorge's Law; Bader, 1954), so we incorporated our modeled changes in the density structure of the firn with the right-hand term in Eqn (2) and solved for the flux divergence,

(2) $$(\nabla _{xy} \cdot \vec q) = \displaystyle{{\dot b} \over {\rho _{\dot b}}} - \displaystyle{{\partial h} \over {\partial t}} - \int_{h_{\rm b}} ^{h_{\rm s}} {\displaystyle{1 \over \rho} \displaystyle{{D\rho} \over {Dt}}{\rm d}z}, $$

where $\rho _{\dot{ b}} $ is the density of the accumulation or ablation, h _s is the glacier surface, h _b is the glacier bed, ρ is the density, and z is the elevation (Cuffey and Paterson, 2010). We used Eqn (2) to solve for the flux divergence over the time period of the surface elevation measurements, from 16 September 2010 to 24 August 2015, using the partitioned ablation and the modeled firn densities to estimate the mass balance and density terms over the same time period. The flux divergence is not necessarily constant through time, so we assumed that seasonal variations in flux divergence scale linearly with seasonal variations in horizontal surface velocities.

The last step for the analysis of the 2010–2015 surface elevation changes was to rearrange Eqn (2) and solve for the ice-equivalent change in surface elevation at the time period of the mass-balance observations, from 22 September 2010 to 19 September 2015. So we utilized the constraint of mass continuity on the relationship between changes in surface elevation and mass balance to maximize the value of the high-quality surface elevation datasets in 2010 and 2015, and then solved for ice-equivalent units over a common time frame.

4.3. Geodetic check

To further check for glacier-wide biases in the mass-balance profiles, we assessed the internal consistency of our results by comparing glaciological and geodetic mass balance over the 2010–2015 interval (e.g. Cox and March, 2004; Huss and others, 2009; Zemp and others, 2013). A geodetic check requires that the calculated changes in glacier thickness and mass balance occur over the same time frame, and that changes in glacier thickness are compensated for changes in density. Our analysis honors these requirements through the process above.

4.4. Ice flux

One way to evaluate the mass redistribution that is represented by the change in ice-equivalent surface elevations is to compare it to the ice flow. To do that we estimated balance and thickness change fluxes, $Q_{\dot{ b}} $ , and Q _∂h/∂t , directly from mass-balance and thickness change measurements, rather than from velocities and ice thicknesses. We then related ice flux, Q, to our mass-balance and thickness change measurements through mass continuity and Eqn (1). The total ice flux through a cross-sectional area is the sum of the balance and thickness change fluxes

(3) $$Q = Q_{\dot b} + Q_{\partial h/\partial t} $$

(Brown and others, 1982). Over the surface of the glacier, Ω, with a range of surface elevations from z _min to z _max, we estimated the balance flux through a cross section at surface elevation, z _i , as the integral of the surface mass balance above the cross section,

(4) $$Q_{\dot b} (z) = \int_{z_i} ^{z_{\max}} {\dot b{\rm d\Omega}} \;. $$

Likewise, we estimated the thickness change flux as the negative of the volume change above the cross section. For a thinning glacier, as is the case at Eklutna, this flux is positive and we refer to it as a thinning flux for clarity, where

(5) $$Q_{\partial h/\partial t} (z) = - \int_{z_i} ^{z_{\max}} {\displaystyle{{\partial h} \over {\partial t}}{\rm d\Omega}} $$

(Rasmussen and Meier, 1982; Nuth and others, 2012).

5. UNCERTAINTIES

Uncertainties in our results stem primarily from measurement errors and from the spatial extrapolation of measured values to unmeasured areas. In this section, we explain our estimates of these uncertainties for glaciological mass-balance and the surface elevation changes. Then we use a sensitivity test to evaluate the potential impacts of sparsely sampled surface mass-balance profiles on our analysis methods.

5.1. Mass balance

At individual mass-balance sites we assumed a nominal measurement uncertainty of 0.2 m w.e. a⁻¹ (Heinrichs and others, 1995; Cox and March, 2004; Beedle and others, 2014) due to mismatches in timing of seasonal mass extrema, measurement errors introduced by surface roughness, bending stakes and snow density variations. To calculate glacier-wide mass balance, the interpolated balance profile between these stakes was extrapolated across the entire glacier. These extrapolation errors are difficult to quantify (Fountain and Vecchia, 1999; Jansson, 1999; Beedle and others, 2014), and even on well-sampled glaciers direct measurements may over- or underestimate the glacier-wide mass balance (e.g. Huss and others, 2009; Zemp and others, 2010).

Several lines of evidence suggest that extrapolation errors on the main branch are likely small, and that extrapolation errors on the west branch could be much larger. First, the upper basin in the main branch is broad with low slope, resulting in homogeneous aspect, slope and shading, which should give rise to similar ablation rates at a given elevation. As we would expect, transient snowlines and final equilibrium lines closely follow elevation contours through the upper basin of the main branch, whereas in the west branch, snowlines are more variable. Second, ground-penetrating radar (500 MHz) measurements from Eklutna Glacier in spring of 2013 show that snow accumulation is primarily a function of elevation with minimal lateral variability – particularly in the main branch. A multi-variable regression between terrain parameters (e.g. elevation, aspect, slope and curvature) and accumulation found that elevation explained most of the glacier-wide variability in accumulation, with a standardized regression coefficient of 0.75 (McGrath and others, 2015). Repeating the same analysis on each branch individually, the elevation coefficient increased to 0.80 on the main branch, supporting the assertion that the distribution of accumulation is largely a function of elevation on the main branch.

Nonetheless, our sparse mass-balance measurements provide little basis for a direct quantitative assessment of extrapolation errors, and we do not report uncertainty in the glacier-wide mass balance. We do assess the glaciological measurements for systematic bias through a geodetic check, presented in Section 6.4. We also use sensitivity tests to evaluate the impacts of mass-balance errors on our analysis, presented in Section 5.3.

5.2. Thinning

We have three surface elevation datasets. The reported or nominal vertical accuracies are ±15 m for the 1957 NED, ±0.3 m for the 2010 LiDAR, and previous studies have found relative accuracies <0.5 m for Ames Stereo Pipeline DEMs of Worldview 3 images (Shean and others, 2016). Rather than relying on unverified reported accuracies, we calculate the uncertainty in the difference between two DEMs directly though analysis of spatial auto-correlation over areas of stable ground (Rolstad and others, 2009; Shean and others, 2016). In any given elevation bin, A, the uncertainty within that bin, σ, is calculated from the spatial distribution of elevation differences over stable ground. We calculated the degree of spatial nonstationarity of variance (expressed as a semivariogram) and assessed the potential for errors as follows:

(6) $$\sigma = \sqrt {c_0 \cdot \displaystyle{{a_0 /\pi} \over A} + \; c_1 \cdot \displaystyle{{a_1 ^2} \over {5 \cdot (A)}}}, $$

where a ₀ is the minimum data spacing, c ₀ is the variance at that minimum spacing (i.e. the nugget), c ₁ is the total variance at distances greater than the range (i.e. the sill), and a ₁ is the range of the sill (Webster and Oliver, 2007; Rolstad and others, 2009).

For the 1957–2010 difference, we infer a larger uncertainty of ±45 m for portions of the accumulation zone with poor contrast on the original photos. This condition has been shown elsewhere to cause large photogrammetric errors (Arendt and others, 2002). We assumed that these potentially substantial measurement errors dominate the uncertainty in the 1957–2010 elevation differences, and adopt them directly as a conservative estimate of the total error within each elevation bin.

5.3. Sensitivity testing

The geodetic check assesses bias in the glaciological measurements; but in our analysis, the geodetic and glaciological mass-balance time series are not completely independent of one another, because the surface mass-balance measurements were used to parameterize the firn model and the adjustments to a common date. Furthermore, the ice-flux analysis is sensitive to the spatial distribution of mass balance, which our geodetic check does not directly test. We address these two issues with sensitivity tests.

The models we used to adjust the calculated thinning rates are dependent, in part, upon the mass-balance record. The models reduce error in our geodetic balances substantially in comparison with the common practice of ignoring these processes altogether (e.g. Fischer, 2011; Das and others, 2014), but introduce the appearance of circularity where errors in the two input datasets could offset each other in the geodetic check. We directly test how much circularity is introduced in this process by shifting the balance profiles by ±1 m and then repeating the geodetic calculations. The propagation of mass-balance errors into the calculated thickness changes is <0.10 m glacier-wide for that ±1 m change in the input mass-balance profiles. Hence, the geodetic check is still largely independent of the mass-balance inputs, and glacier-wide errors in the glaciological balances could only be explained by independent, similar-magnitude, systematic errors in the measured surface elevations.

A geodetic check does not necessarily confirm the glaciologically calculated spatial distribution of mass balance. The mass-balance profile could contain compensating errors and be either steeper or flatter than our estimate. Because such errors would have important consequences for our interpretation of the thinning rates within the upper basin, we experimentally determined how much error would be required in our estimate of the upper basin's mass balance to account for the enhanced thinning we found there. The ~3 m of thinning we observed in the upper basin over a 5-year period could be explained by surface mass balance alone if we overestimated the upper basin balance by 0.6 m w.e. a⁻¹. Errors of that magnitude in the local mass-balance profile are plausible given the sparsity of mass-balance observations. But glacier-wide mass continuity would require an opposite error of 15 m (over the 5-year period) averaged over the entire area below the upper basin. In other words, compensating mass-balance errors could explain the observed thinning only if we also underestimated the balance by 3 m w.e. a⁻¹ over the much smaller area below the upper basin. This scenario would require net accumulation over most of the ablation area, and ablation in the upper basin where seasonal snow remains at the end of the summer ablation season and can thus be ruled out. Consequently errors in the mass-balance profile within the upper basin are likely to be very small or compensate each other within the upper basin, and in either case the flux analysis is still fundamentally correct.

6. RESULTS

First, we provide results from the glaciological mass-balance observations at the point, branch and glacier-wide scales for 2008–2015. Second, geodetic solutions and the distributions of thinning are presented for 1957–2010 and for 2010–2015. Additional data are presented in the supplement.

6.1. Glaciological mass balance

Between 2008 and 2015 we observed a cumulative glacier-wide mass balance of −4.7 m w.e. (Fig. 3a). Measured seasonal point balances ranged from +2.8 ± 0.2 to −5.2 ± 0.2 m w.e. a⁻¹ (data at doi.org/10.5066/F7MP51CB). Inferred annual balance profiles generally had consistent shapes from year to year and similar relationships between the branches (Fig. 2). Glacier-wide annual balances ranged from −1.4 to +0.6 m w.e. a⁻¹, with the most positive balances in 2008 and 2012 and the most negative in 2013 and 2015 (Fig. 3a, Table 1).

Table 1. Branch-wide and glacier-wide mass-balance results (m w.e. a⁻¹)

6.2. Thinning

Glacier thickness declined in most areas over both intervals. Raw (coregistered but not including adjustments for firn density and seasonal aliasing) glacier-wide average surface elevation change between 1957 and 2010 was −0.62 ± 0.55 m a⁻¹. From 2010 to 2015 average surface elevation change was −0.84 ± 0.11 m a⁻¹ (Fig. 4).

Fig. 4. Raw surface elevation change on Eklutna Glacier for (a) 1957–2010 and (b) 2010–2015. The color bars indicate elevation changes in meters per year. White areas within the glacier outline indicate data gaps. The upper basin of the main branch (>1360 m in 2010) is demarcated by the black line. Each panel has the same extent as Figure 1c, and shows the extent of the glacier corresponding to the first year of the interval.

Over the 2010–2015 interval our estimates of the change in firn density, late season ablation, and ice flow were generally small (Fig. 5; underlying data in Figs S3–S5), but they had a net effect of reducing the thinning rate in the upper elevations and consequently provide a more conservative estimate of the distributed thinning. With minor exceptions at the highest elevations, the annual ice-equivalent surface elevation changes from 1957 to 2010 and 2010 to 2015 are negative over most of the glacier's elevation range, including the broad upper basin that dominates the accumulation area of the main branch (Fig. 6). Near the terminus, surface elevation change rates are near −4.0 m a⁻¹, and approach 0 at highest elevations in both branches. In the broad upper basin in the main branch the change in surface elevation is −1.0 ± 0.8 m a⁻¹ from 1957 to 2010, and −0.6 ± 0.1 m a⁻¹ from 2010 to 2015. Within the main branch the glacier is thinning at all elevations below 1550 m, encompassing 92% of the glacier area. In the west branch, by comparison, the changes in surface elevations near 1600 m were −0.2 ± 0.8 m a⁻¹ from 1957 to 2010 and −0.6 ± 0.1 m a⁻¹ from 2010 to 2015, but the total area below 1600 m includes only 65% of the branch area. Analysis of elevation changes from the additional surface elevation datasets in 2007, 2011, 2012 and 2014 reveal similar patterns of thinning and are presented in the supplement.

Fig. 5. Measured and calculated changes in surface elevation, by 50 m elevation bin, for the main and west branches of Eklutna Glacier between 2010 and 2015. Both panels show the 2010 glacier surface hypsometry in gray with values on the upper horizontal axis. (a) Red x's show the mean value of the raw surface elevation difference 2010 to 2015. These are derived from measurements shown in Figure 4. Error bars reflect measurement uncertainty. Pink lines show the difference in ablation that occurred between the surface elevation measurement date and fall mass-balance measurement date in 2010 and 2015. Blue lines show the difference in elevation changes due to ice flow (from flux divergence rates) between the surface elevation and mass-balance measurement dates in 2010 and 2015. Cyan lines show the change in firn pore-space due to changes in firn density profiles during the interval from 2010 to 2015. (b) The same adjustments to the measured surface elevation differences as panel (a), shown with increased scale for detail.

Fig. 6. Mean annual thinning rates, by 50 m elevation bin, for the main and west branches of Eklutna Glacier, for 1957–2010 and 2010–2015. The shaded areas indicate the uncertainty. The bottom of the main branch upper basin is shown by the dashed line.

6.3. Geodetic check

The geodetic mass balance from 1957 to 2010 is −0.52 ± 0.46 m w.e. a⁻¹ (Table 2). From 2010 to 2015 the glacier-wide geodetic mass balance is −0.74 ± 0.10 m w.e. a⁻¹, which is indistinguishable from the glaciological result of −0.73 m w.e. a⁻¹. Evaluating the main and west branches separately we found −0.72 ± 0.08 and −0.77 ± 0.12 m w.e. a⁻¹, which are also not distinguishable from the respective glaciological mass balances of −0.67 and −0.81 m w.e. a⁻¹ (Table 2). Because the surface elevations are both better constrained and more widely distributed than our mass-balance measurements, this agreement between the glaciological and geodetic balances supports the glaciological mass-balance time series, and provides no justification for a geodetic calibration.

Table 2. Comparison of the mean glaciological and geodetic mass balances (m w.e. a^–1) for 2010–2015

6.4. Ice flux

Balance flux, thinning flux and total ice flux for 2010–2015 show that the thinning rates in the upper basin of the main branch of Eklutna Glacier account for half of the total flux out of the basin (Fig. 7). In the main branch the balance flux becomes negative below 1390 m, 6 km upstream of the present terminus and within the upper basin. The thinning flux is greater than the balance flux at the present ELA and increases rapidly through the upper basin, suggesting the lower glacier is presently sustained through thinning of the upper basin. The long-term thinning flux (1957–2010) is similar, although with greater uncertainties. In the west branch, the balance flux becomes similarly negative below 1397 m, although the thinning flux from the higher elevations (>1500 m) makes only a minor contribution to the total ice flux. The main contribution is from thinning along the lower trunk of the glacier, well below the present ELA.

Fig. 7. The mean annual ice, balance and thinning fluxes 2010–2015 (blue), and the mean annual thinning flux for 1957–2010 (red). Uncertainties for the thinning fluxes are shown by the shaded areas.

7. DISCUSSION

We measured a cumulative mass balance of −4.4 ± 0.2 m w.e. at Eklutna Glacier from 2010 through 2015, despite the positive balance year in 2012. We compared this mass-balance record with surface elevations measured in 1957, 2010 and 2015, and found mean annual ice-equivalent surface elevation change of −0.59 ± 0.55 and −0.84 ± 0.11 m a⁻¹ for the 1957–2010 and 2010–2015 periods, respectively. The mass loss manifests as thinning over most of the glacier hypsometry, including much of the accumulation zone within the main branch's upper basin. Changes in firn density, and the late season ablation and ice flow reduced the apparent thinning at higher elevations, but the rates of surface elevation change in the main branch were still −0.6 ± 0.1 m a⁻¹ or faster, indicating thinning over 92% of the branch area. This pattern was similar to that from 1957 to 2010, albeit with higher uncertainties.

7.1. Hypsometry

The shift in the glacier hypsometry to lower elevations has implications for future mass balance and the glacier-derived water supply. In the main branch, the relatively modest thinning rates in the upper basin represent a large disconnect between surface mass balance and ice flow that appears to be both ongoing and long-lived. Another way of saying this is that the mass loss in the main branch is primarily manifest as thinning within the upper basin, which has mitigated terminus retreat. This process is analogous to a surge, albeit over a much longer period of time and at a much slower rate. The ice flux analysis suggests that even without continued thinning or additional climate change, the supply of ice to the lower glacier is unsustainable. Ultimately, the equilibrium state with current climate may be governed by how long the positive feedback between mass balance and thinning persists, but the large discrepancies between the balance and thinning fluxes suggest that the glacier will retreat until it occupies only a small fraction of its present-day extent.

Many mountain glaciers have high mass turnover rates and steep surface slopes, both of which promote short response times and a close linkage to the present climate (Harrison and others, 2001; Oerlemans, 2001). The broad, low-slope upper basin on the main branch is both large and critically located coincident with the present range of ELAs, which make this glacier more susceptible to a significant altitude-mass-balance feedback (Cuffey and Paterson, 2010; Huss and others, 2012).

The thinning is resulting in a downward shift in the hypsometry, decreasing the amount of area above the ELA and increasing the area below it (Fig. 8). The persistent thinning in the upper basin suggests that Eklutna is an intermediate between steeper mountain glaciers and something like the Yakutat Icefield in Southeast Alaska, where the glacier has thinned almost completely out of the accumulation zone (Trüssel and others, 2013). Unlike the Yakutat Glacier, which will likely disappear completely over the next 60–100 years (Trüssel and others, 2015), the highest elevations in the Eklutna basin may maintain an active cirque glacier. However, <10% of the main branch area both has a stable surface elevation and is above the 2008–2015 mean ELA. If we assume that Eklutna Glacier would need a typical accumulation area ratio of 0.58 (Dyurgerov and others, 2009) to return to equilibrium with the current climate, then the reduction in area would exceed 80%.

Fig. 8. 2010–2015 change in hypsometry. The change in area within each elevation bin expressed as a percent of total branch area. Individual lines are shown for three surfaces to show the effect of uncertainty in surface elevation change. Mean ELAs for 2010–2015 are shown in black.

7.2. Stability

The rate and spatial pattern of thinning outweighs the stabilizing influence of terminus retreat. Consider the impact of observed geometric changes in the main branch Eklutna Glacier (we exclude the west branch from this analysis because of its greater uncertainties, and because it does not exhibit the accumulation zone thinning under discussion here) on a static mass-balance profile between 2010 and 2015. The balance profiles from those 2 years provide two snapshots of glacier climate: one from a near equilibrium year (2010 main branch, glacier-wide mass balance −0.29 m w.e. a⁻¹), the other from a negative year (2015 main branch, −1.33 m w.e. a⁻¹). We used first one, then the other balance profile to calculate the glacier-wide mass balance on both the 2010 and 2015 glacier surfaces, as if the glacier was subjected, at each of those times, to exactly the same climate. In both cases, the results are similar: between 2010 and 2015, the glacier's main branch geometry changed in such a way as to change the mass balance by −0.04 m w.e. a⁻¹. While the magnitude of the change over that time period is small, it does show that the positive feedback associated with glacier thinning has in recent years been greater than the negative feedback of terminus retreat. This is an unstable response (Böðvarsson, 1955). A comparable analysis for the past 50 years is less robust, due to larger uncertainties on the 1957 DEM, but is consistent with the possibility that Eklutna Glacier has been undergoing unstable retreat for at least a half-century.

7.3. Water resources

Our results show that Eklutna Glacier, like most mountain glaciers in Alaska (Arendt and others, 2002; Loso and others, 2014; Larsen and others, 2015), has lost considerable mass over the last half-century. Already, the consequences of these changes for downstream water users have been substantial. Most directly, changes in ice volume have impacted total runoff in the basin. If we assume that the Anchorage Municipal Light and Power mean reported inflow into Eklutna Lake from 2000 to 2009 of 0.304 km³ represents a long-term average, then ice loss enhanced cumulative reservoir inflow by 5 ± 4% from 1957 to 2010 and 7 ± 1% from 2010 to 2015. Given the uncertainties we cannot say if 2010–2015 represents an actual increase in the contribution from negative mass balances, but it is clear that negative mass balances have made at least some contribution over the long term. Annual contributions were ~13% in 2013 and 2015. This deglaciation discharge dividend will ultimately diminish as the shrinking glacier eventually returns to a rough equilibrium with the new climate and annual mass balances trend towards zero (e.g. Fountain and Tangborn, 1985; Collins, 2008; Huss, 2011).

Future runoff predictions will need to understand how long the instability will continue. Both the glacier size and the ice-flux analysis are consistent with a substantial deglaciation discharge dividend over the next few decades, but the eventual loss of the deglaciation discharge dividend is unavoidable over longer timescales. Our results demonstrate that predicting the evolving geometry in detail will require understanding the ice dynamics governing upper basin thinning, and that accurate forecasts will require coupled models that capture the changing glacier hypsometry (Schäfer and others, 2015).

8. CONCLUSIONS

The hypsometry of Eklutna Glacier, with 60% of the area in a broad, low-slope basin near and above the recent mean ELA, is driving a response to climate that has implications for both the future of the glacier and for downstream runoff. The cumulative glacier-wide mass balance from 2010 to 2015 was −3.7 ± 0.2 m w.e., with annual balances ranging from −1.4 m w.e. a⁻¹ (in 2013 and 2015) to +0.6 m w.e. a⁻¹ (in 2012). Balance-flux estimates suggest that over 6 km of terminus retreat is required to approach equilibrium with the present climate. With a continued altitude-mass-balance feedback the glacier could lose >80% of its present area, and possibly more with continued climate change.

The mass lost from Eklutna Glacier has important consequences for runoff to a downstream reservoir that is used – in its entirety – to provide municipal water and hydropower for the city of Anchorage. The runoff from net-mass loss alone averaged 7 ± 1% of total inflow to the reservoir from 2010 to 2015, and in 2013 and 2015 reached ~13%. The eventual diminishing deglaciation dividend will necessarily result in decreased water availability. Our results corroborate the growing number of studies warning of the consequences of glacier retreat for downstream water supplies, but suggest the use of caution before making simplistic assumptions about the nature of that retreat.