21st-century increase in glacier mass loss in the Wrangell Mountains, Alaska, USA, from airborne laser altimetry and satellite stereo imagery

Abstract

Alaskan glaciers are among the largest regional contributors to sea-level rise in the latter half of the 20th century. Earlier studies have documented extensive and accelerated ice wastage in most regions of Alaska. Here we study five decades of mass loss on high-elevation, land-terminating glaciers of the Wrangell Mountains (~ 4900 km2) in central Alaska based on airborne center-line laser altimetry data from 2000 and 2007, a digital elevation model (DEM) from ASTER and SPOT5, and US Geological Survey topographic maps from 1957. The regional mass-balance estimates derived from center-line laser altimetry profiles using two regional extrapolation techniques agree well with that from DEM differencing. Repeat altimetry measurements reveal accelerated mass loss over the Wrangell Mountains, with the regional mass-balance rate evolving from –0.07 ± 0.19 m w.e. a–1 during 1957–2000 to –0.24 ± 0.16 m w.e. a–1 during 2000–07. Nabesna, the largest glacier in this region (˜1056 km2), lost mass four times faster during 2000–07 than during 1957–2000. Although accelerated, the mass change over this region is slower than in other glacierized regions of Alaska, particularly those with tidewater glaciers. Together, our laser altimetry and satellite DEM analyses demonstrate increased wastage of these glaciers during the last 50 years.

Introduction

Alaskan glaciers cover ˜12% of the world’s glacier area outside of the polar ice sheets and are major regional contributors of sea-level rise (Arendt and others, 2002; Gardner and others, 2013). Among the Alaskan glaciers, the largest contribution comes from the maritime glaciers around the Gulf of Alaska (Motyka and others, 2003; Muskett and others, 2003; Arendt and others, 2006; Larsen and others, 2007; Berthier and others, 2010). Higher sensitivity to increasing temperatures of these maritime glaciers, often located at lower altitudes (e.g. Hock and others, 2009), and effects of tidewater glacier dynamics (O’Neel and others, 2005; McNabb and others, 2012) account for the extensive ice wastage over these regions. The mass-balance rates of the maritime glaciers of the Coast Mountains (–0.65 ± 0.14 m w.e. a^–1) and the Chugach Mountains (–0.64 ± 0.07 m w.e. a^–1) during the period 1962–2006 are twice as negative as those of the more continental glaciers in the Alaska Range (–0.30 ± 0.09 m w.e. a^–1; Berthier and others, 2010). The focus of our study is recent mass change of the Wrangell Mountains glacier complex. These mountains are tall, form an orographic barrier to the storms in the north Pacific and are characterized by high amounts of precipitation (Benson and others, 2007; Kanamori and others, 2008).

The Wrangell Mountains are located between the southeast Alaska maritime glaciers and the continental glaciers of the Alaska Range. Very few studies have attempted to address the mass balance of this region separately (Arendt and others, 2002, 2013). A mass-balance rate of –0.17 ± 0.07 m w.e. a^–1 during 1957–2000 was inferred by comparing laser altimetry results to the US Geological Survey (USGS) elevation contour map, using a single large glacier (Nabesna) to represent the whole region (Arendt and others, 2002). A recent GRACE (Gravity Recovery and Climate Experiment) mascon analysis covering the period 2004–10 indicated that summer, winter and annual mass-balance variability and mass loss was smaller over the Wrangell Mountains than in other regions within the Gulf of Alaska (Arendt and others, 2013). However, the GRACE mascon solution also included a small portion of the St Elias Mountains and may be biased by signal leakage from surrounding ice masses. In the pan-Alaska assessments of multi-decadal mass change by Berthier and others, (2010), the continental glaciers of the Wrangell Mountains were combined with the more coastal glaciers of the St Elias Mountains. Thus, the negative regional mass-balance rate of –0.47 m w.e. a^–1 for the combined ice area of the Wrangell and St Elias Mountains was driven mainly by the rapid mass loss from the large maritime glaciers of the St Elias Mountains (e.g. Bering Glacier) (Berthier, 2010), glaciers flowing into Icy Bay (Muskett and others, 2008) or the Yakutat Icefield glaciers (Trüssel and others, 2013). Furthermore, none of these papers has examined different time periods to assess whether the glacier mass-balance rate of the Wrangell Mountains has evolved with time. A rationale for studying this glacier complex separately is the lack of tidewater glaciers, whose dynamics tend to complicate the understanding of mass-balance response to climate fluctuations.

In this study, we use airborne laser altimetry measurements over three glaciers, 1957 USGS maps and satellite stereo imagery to estimate the mass-balance rate of the Wrangell Mountains glacier complex between 1957, 2000 and 2007. We also test if center-line laser profiling of three glaciers is able to correctly capture their glacier-wide and the regional mass-balance rates. Elevation changes derived from laser altimetry data are extrapolated to the whole region using a technique in which we normalize the elevation range prior to extrapolation. We then compare these results with the more commonly used elevation-dependent extrapolation technique and also with sequential digital elevation model (DEM) analysis from ASTER (Advanced Spaceborne Thermal Emission and Reflection Radiometer) and SPOT5 (Satellite Pour l’Observation de la Terre) stereo imagery. We investigate three periods (1957–2000, 1957–2007 and 2000–07) to assess potential changes in glacier-wide and regional specific mass-balance rate.

Study Area

The Wrangell Mountains, central Alaska (˜61.0–62.45º N, 141.0–145.0º W; Fig. 1), are covered by an ice area of ˜4900 km² spanning 450–4998 m a.s.l. The vast majority of the glaciers in this range are land-terminating, but proglacial lakes have developed at the front of some of the glaciers.

Fig. 1. Wrangell Mountains glaciers with outlines from 1957 and laser altimetry flight lines from 2007 and 2000, color-coded to show location of the profiles. The surveyed glaciers are labeled Nabesna, Kennicott and Nizina. Colored shading marks unmeasured subregions used in extrapolation method B. The purple profile over Nabesna Glacier is located closest to the Wrangell volcano. Contours from the 1957 USGS DEM are displayed every 500 m.

We have laser altimetry data from three glaciers in this region: Nabesna, Kennicott and Nizina–Regal–Rohn (here-after referred to as Nizina for simplicity) glaciers. Nabesna (˜1056 km²), the largest glacier in this region, is also among the largest land-terminating glaciers in the world. It is˜120 km long and fed by more than 40 tributaries. Li and others, (2008) have observed short-term velocity changes of Nabesna Glacier using synthetic aperture radar interferometric techniques. Kennicott Glacier (˜393 km²) is a 43 km long glacier located on the south-facing slope of the Wrangell Mountains. Bartholomaus and others (2008, 2011) have studied the effects of an active subglacial water system on the short-term surface elevation and velocity changes of the glacier in summer. Nizina Glacier (˜429 km²) is also located on the south-facing slope of the Wrangell Mountains. Together, the three profiled glaciers cover˜1855 km², corresponding to ˜39% of the total ice area.

The north-facing glaciers in this mountain range tend to terminate at elevations greater than 800 m a.s.l., while glaciers on the south side terminate at lower elevations (e.g. Kennicott Glacier at ˜425 m a.s.l. and Nizina Glacier at ˜650 m a.s.l.) (Fig. 1) due to considerably higher snow accumulation rates on the ocean-facing side (Benson and others, 2007). The annual precipitation on the ocean-facing side is 3–5 m w.e. a^–1, several times larger than on the north side (Benson and others, 2007). Precipitation on the Wrangell Mountains is episodic in nature, correlates with the seasonality and timing of the individual storms in the coastal regions and is due to regional weather systems rather than the local convective systems common to mountainous regions (Kanamori and others, 2008).

The mountain range is volcanic in origin, with an active shield volcano (Mount Wrangell) and several old and eroded volcanoes. On the summit caldera of the Wrangell volcano, a basal ice melting rate of –0.64 m a^–1 and a volcanic heat flux of 7 W m^–2 were inferred from ice-penetrating radar measurements (Clarke and others, 1989). Glacier calorimetric measurements at this location also suggested increased volcanic heat flux along the north crater (20–140 W m^–2) within the past decade (Benson and others, 2007). The mean annual temperature measured 10 m below the glacier surface of the summit caldera is ˜20ºC, so the ice melt observed in the summit caldera is exclusively due to basal melt (Benson and others, 2007). Therefore, some of the glaciers may have unusual basal rheology due to their proximity to an active volcano (Sturm and others, 1991).

Datasets

Airborne laser altimetry

Laser altimetry profiles were acquired on Nabesna and Kennicott Glaciers on 2 June 2000 and 20 June 2007, while Nizina Glacier was profiled only on the latter date. Altimetry data were collected along glacier center lines for most of the major tributaries of these glaciers. The flight lines for both years and the outline of the profiled glaciers are shown in Figure 1.

The airborne laser altimeter profiling system consisted of a nadir-pointing infrared range-finder, a gyro and a compass mounted on a small aircraft. The airplane was flown 20– 100 m above the glacier surface (Sapiano and others, 1998). The range-finder measures the distance from the aircraft to the glacier surface. The gyro and the compass measure the orientation of the aircraft every second. The absolute positions of the aircraft are determined by kinematic GPS. The laser altimetry shots are spaced at 1.2 m intervals on the surface at a typical aircraft speed of 30 m s^–1. The beam diameter is 0.18 m at a distance of 100 m (Echelmeyer and others, 1996). The altimetry profiles were usually acquired on clear days for better GPS solutions and accurate elevation measurements. The accuracy of GPS solutions depends on the number of satellites available, and, for the present analysis, measurements with error of >0.2 m are not considered. The accuracy of the profiling system is typically 0.3 m depending on the surface slope of the glacier. It decreases as the laser encounters steep slopes or icefalls (Echelmeyer and others, 1996; Arendt and others, 2008). All data points are referenced in International Terrestrial Reference Frame 2000 (ITRF00), and coordinates are projected to a Universal Transverse Mercator (UTM) projection (ellipsoid WGS84), zone 7. Elevation data are recorded as height above ellipsoid. The vertical coordinates of the altimetry profiles were transformed using the US National Geodetic Survey model GEOID 99 (Alaska).

USGS DEM

For multi-decadal elevation changes, altimetry profiles and DEMs derived from satellite stereo imagery are compared to USGS topographic maps of 1: 63 000 scale established from aerial photographs acquired in 1957. These maps have a contour interval of 100 ft (30.48 m). The 15 min USGS DEMs used for elevation change calculations in this paper are derived from these maps and are directly provided by the USGS with a 2^” grid spacing, corresponding to about 60 m in latitude and 30 m in longitude. The horizontal and vertical datums of the USGS DEMs are respectively in NAD 27 and National Geodetic Vertical Datum (NGVD 29) elevations relative to mean sea level. The USGS DEMs are converted to WGS84 to be consistent with the profiles. To make sure that all data share the same reference for elevation, we have converted the NGVD29 vertical datum to EGM96 by subtracting 2.5 m from the USGS DEMs. This empirical value for the vertical datum conversion is based on the average elevation difference between the USGS DEMs and EGM96 along Ice, Cloud and land Elevation Satellite (ICESat) altimetry profiles (Zwally and others, 2002) over ice-free surfaces in all southeast Alaska (Berthier and others, 2010). The vertical offset applied here, 2.5 m, is the mean for all southeast Alaska. This vertical offset varies by less than ±1.5 m when computed separately for the different glacierized subregions considered by Berthier and others, (2010). EGM96 and GEOID99 over Alaska are similar (personal communication from D. Roman, 2013), so the exact same USGS DEM can be used in both laser altimetry and DEM differencing analyses.

DEMs from ASTER and SPOT5 stereo-imagery

For sequential DEM analysis of ice volume change, we use four DEMs derived from SPOT5-HRS (Korona and others, 2009) and ASTER (Fujisada and others, 2005) pairs of stereo images acquired between 2003 and 2006 (Table 1). When available, images as close as possible to the end of the ablation period (mid-September in Alaska) are selected to minimize errors owing to seasonal elevation changes. The SPOT5 DEM has a grid spacing of 40 m, the ASTER DEMs of 30 m. Three of these DEMs are identical to those used in a previous pan-Alaskan volume change study (Berthier and others, 2010). For the present study, we include an additional DEM derived from one ASTER stereo pair acquired on 16 September 2003 to improve coverage in the eastern part of the Wrangell Mountains, around Nizina Glacier. Both ASTER and SPOT5 DEMs were automatically derived from stereo imagery without ground control points and thus contain artifacts and some planimetric/altimetric biases. For all satellite DEMs, unreliable elevations due to clouds in the imagery, shadows and lack of features in the accumulation areas affected ˜11% of the ice-covered area sampled by the DEMs. They are masked using the score channel that indicates the quality of the correlation during DEM generation. Planimetric shifts between the SPOT5/ASTER DEMs and the USGS DEMs are corrected by minimizing the standard deviation of their elevation differences off glaciers (e.g. Berthier and others, 2007). Vertical biases are estimated on and off glaciers using the ICESat data acquired closest in time to the acquisition date of each satellite DEM (Table 1).

Table 1. Characteristics of the Spot5-HRS and ASTER stereo-images used to derive DEMs and statistics of the DEM calibration against ICESat data. ‘ICESat’ indicates the ICESat laser periods used to evaluate the satellite DEMs. ‘Mean’ and ‘Std dev.’ refer to the mean and standard deviation of the elevation differences between the DEMs (before calibration) and the ICESat data, respectively. N is the number of points where the differences have been computed

Area–altitude distribution

The outlines of the glacier complex are from 1957 and 2010 (Randolph Glacier Inventory version 3.0, Arendt and others, 2012). The area in each elevation band is calculated from the 1957 USGS DEM and glacier outlines to obtain the initial area–altitude distribution. Landsat images from 2010 are used to update the glacier outlines and calculate area changes and front retreat rates between 1957 and 2010. The 2010 glacier outline is also used to create a new area– altitude distribution from the older 1957 DEM after correcting for the area changes. For volume-change estimates, an average of the original (1957 USGS DEM) and the updated area–altitude distribution (Landsat 2010) is used. The total area used in the volume-change calculations consists of averages of the original, (1957) and the updated Landsat, (2010) area.

Methodology

Elevation changes along altimetry flight lines

Elevation changes over multi-decadal timescales (1957– 2000, 1957–2007) are calculated by differencing the altimeter-measured surface elevations from those of the 15 min USGS DEMs. Earlier studies (e.g. Sapiano and others, 1998) only used the profiling data where the ground track of the profile intersected a map contour. Here we use all available laser shots along the flight line. The USGS DEM is interpolated at the location of the laser shots for elevation values from 1957. The values of the interpolated DEM are then differenced from the corresponding elevation values obtained from the laser altimetry to estimate elevation change along the flight lines. For profile-to-profile elevation change estimates, we identify all data points of the later flight line, (2007) that fall within a 20 m wide circular window centered over the earlier data point. These data points are used to calculate the elevation change by subtracting it from the older profile. The elevation changes within a 30.48 m contour interval of the original USGS contour map are then averaged to create an elevation-dependent profile of elevation changes for each glacier.

Elevation changes along each laser altimetry profile of the three surveyed glaciers for the different time periods indicate large small-scale spatial variability in the elevation changes (Fig. 2a–f). This small-scale variability may be due to spurious returns from small crevasses, debris, steeper surface slopes, low-lying thin cloud or fog and errors in the USGS DEMs. There are also significant differences in the elevation changes between different profiles that are located within the same elevation band but on different tributaries of a surveyed glacier. For most elevation ranges, we smooth the data for each glacier by averaging all elevation change values of all available profiles within each 30.48 m contour interval of the original USGS maps to obtain the elevation change ∆z_i for each elevation band, i, of the same glacier. A mean ∆z curve is compiled for each of the three profiled glaciers from the ∆z_i values from all tributaries of the glacier and is used to obtain the glacier-wide mass balances.

Fig. 2. Elevation-change rates for all available laser shots over the three profiled glaciers for three time periods. The colored lines refer to the laser altimetry profiles along various tributaries (same colors as the laser profiles in Fig. 1). The black curve is the final elevation-change curve used for computing glacier-wide balances. The green curves in (c), (f) and (g) are cumulative area derived from 30.48 m elevation bins (% of total area) as a function of elevation.

In some cases, manual adjustments are necessary to determine the mean ∆z curve and these are obtained after carefully scrutinizing elevation changes from all available profiles. In the ablation zone (lower elevation), sometimes an elevation band does not have enough points due to loss of laser returns from a steep or highly crevassed region (e.g. Fig. 2f at 675–750 m a.s.l.). In such lower-elevation cases we linearly interpolate the elevation change values from the neighboring bands.

At upper reaches of the glaciers, multi-decadal elevation change values with very high scatter and those that indicate unrealistic thinning or thickening are discarded. In these cases and for upper elevation bands lacking reliable data, we assume no elevation change (e.g. Fig. 2a and b). Where the elevation change pattern looks unrealistic, we cross-check the pattern with repeat altimetry change rates of 2000–07 to identify unreliable measurements. Over such areas, only multi-decadal elevation change patterns that are consistent with those observed in the repeat altimetry measurements are retained in the final ∆z curve. For example, we discard the blue profile (Fig. 2a) that shows anomalous thickening and thinning rates during 1957–2007 between 2000 and 3000 m a.s.l. on Nabesna Glacier. This profile is located on the western tributary, originating from the summit caldera of the Wrangell volcano (Fig. 1). This profile was surveyed for the first time in 2007, so we do not have profile-to-profile data for comparison. As this profile is located on a tributary originating from the summit caldera of the Wrangell volcano, the convective heat flow due to the volcano may affect the basal melt rates and ice dynamics and produce anomalous surface elevation changes along this profile. However, as this hypothesis lacks evidence, we discard the elevation changes along this profile. The unrealistic elevation changes may be due to possible errors in the 1957 USGS maps resulting from lack of topographic controls and surface texture in these regions (e.g. Aðalgeirsdóttir and others, 1998).

We also assume that the surface-elevation changes are equal to glacier thickness changes, hence we neglect any elevation changes that may occur due to isostatic rebound or subglacial erosion/deposition or changes in glacier surface elevation due to subglacial water movement.

Glacier-wide mass-balance rates of profiled glaciers

To extrapolate the measured elevation changes to the entire glacier, we assume that the elevation change along the center line also represents the entire basin of the glacier (Arendt and others, 2002, 2006). We assess the validity of this assumption by comparing our elevation change estimates with those obtained using the satellite DEM differ-encing method (see Discussion).

The volume change for each glacier is calculated by summing the product of the elevation change for each elevation band, ∆z _i , and the area in the corresponding elevation band, A_i . The elevation of each band is taken from the original 30.48 m contours from the 1957 USGS DEM. The volume change is translated into a mean glacier-wide mass-balance rate Ḃ (m w.e. a^–1) by

(1)

where ʹp _i and p _w are the densities of glacier ice and water respectively, A _total is the glacier area, I is the number of elevation bands and n is the number of years between surveys. For regional mass-balance estimates, we have used 4847.25 km² as the total glacierized area of the Wrangell Mountains. This area is the mean area from the USGS maps, (1957) and the satellite images, (2010).

The density factor used to convert volume change to mass change estimates can vary widely depending on the climate conditions during the time interval of the observations, and is difficult to estimate. Assuming a constant density of 900 kg m^–3, as frequently done in the literature, may lead to a systematic overestimation of mass loss when firn area or thickness has decreased (Sapiano and others, 1998; Huss, 2013). Following these two studies, we use a conversion factor of p _i = 850 ± 60 kg m^–3. Huss, (2013) recommended this value for periods longer than 5 years, with stable mass-balance gradients, presence of firn area and volume changes significantly different from zero.

Regional extrapolation techniques

We apply two techniques of regional extrapolation to derive the volume change and mass balance of the unmeasured glaciers in the Wrangell Mountains from the measured glaciers.

Method A: elevation-dependent extrapolation

This technique has been used extensively in previous studies for mass-balance estimates from airborne laser altimetry (e.g. Sapiano and others, 1998; Arendt and others, 2002, 2006). The measured elevation change within an elevation band is applied directly to the same elevation band of the unmeasured glaciers. The specific mass balance is computed according to Eqn (1) using ∆z from the measured glaciers and A_i and I from the unmeasured areas.

Method B: extrapolation using normalized elevation range

It is commonly observed that glaciers, irrespective of their altitudes, have their highest thinning rates close to the terminus. This poses a problem for elevation-dependent extrapolation where the measured glaciers have different terminus elevations such as in the Wrangell Mountains. Therefore, we normalize the elevation range (Schwitter and Raymond, 1993; Arendt and others, 2006; Johnson and others, 2013) according to

(2)

Here Z is the normalized elevation, z is the elevation of the original USGS maps and z _h and z _t are the elevations at the head and terminus of the glacier respectively. Using Eqn (2), each measured glacier is divided into 100 normalized elevation bins ranging from 0 to 1 with a spacing of 0.01. Thus all the measured glaciers are of equal normalized length. The mean ∆z curve of each glacier is obtained by taking the mean of the elevation changes within the normalized elevation bands, and volume change is calculated for each measured glacier using Eqn (1) (Fig. 3).

Fig. 3. Elevation-change rates for each 30.48 m elevation bin versus elevation (a–c), and normalized elevation range (e–g) for the three surveyed glaciers and time periods, used for calculating the glacier-wide balance of the surveyed glaciers and for extrapolation to unmeasured areas of the Wrangell Mountains. Area–altitude distributions for various domains are given in (d) and (h).

For extrapolation to unmeasured glaciers, we group the unmeasured glaciers based on their terminus elevations. We identify five subregions with glaciers of similar terminus elevations (Fig. 1). The elevations of each of these subregions are normalized using Eqn (2) such that each region again has 100 elevation bands ranging from 0 to 1. A new area– altitude distribution is obtained based on the 0.01 spacing of the normalized elevation bands. The mean ∆z curve compiled from the elevation changes of the three measured glaciers is then applied to the new area–altitude distributions calculated for each of the five subregions. The volume change and the specific mass-balance rate are calculated using Eqn (1).

Mass-balance estimates from sequential DEM analysis

The elevations from the 1957 USGS DEM are subtracted from those of the recent satellite DEMs to estimate the glacier volume change. There are data gaps over ˜11% of the ice-covered area due to poor correlation during the generation of the satellite DEM or inconsistent elevation contours in the USGS maps, mostly in the textureless accumulation area. Over these data gaps, the elevation change of the unmeasured areas is assumed to equal the regional mean elevation change at the same altitude. Volume changes are converted to mass change assuming a mean density of 850 kg m^–3.

Error estimation from the laser altimetry analysis

Errors in the USGS topographic maps are one of the largest sources of inaccuracies in volume change estimation and include both random and systematic errors. The random errors are calculated as in previous studies (Arendt and others, 2002, 2006). The vertical errors involved in the USGS maps (E _USGS) are ±15 m in the ablation zone and ±45 m in the accumulation zone (Aðalgeirsdóttir and others, 1998; Arendt and others, 2006). The uncertainty is large if the ice area is relatively flat and devoid of rock outcrops and consequently introduces error due to lack of contrast leading to ‘floating contours’. This is the largest source of error in the USGS DEM (Arendt and others, 2002). The altimetry system error (E _LA) is ±0.3 m. Ambiguities involved in the dates of the aerial photographs used in making the USGS topographic maps (E _MD) are ±2.5 m (Arendt and others, 2002). The volume-to-mass conversion error (E _D) is ±60 kg m^–3 (Huss, 2013). These errors are allowed to propagate through the elevation bins and then summed in quadrature to estimate the random errors involved in the glacier volume changes.

Extrapolation of measured elevation changes to unmeasured areas is another major source of error in the laser altimetry analysis. Those extrapolation errors are due to (1) the assumption that one or a few center-line profiles represent the whole glacier, i.e. the profile-to-glacier error (E _PG), and (2) the assumption that the three surveyed glaciers are representative of the unmeasured glaciers, i.e. the regional extrapolation error (E _EXT). To quantify the profile-to-glacier error, we compute the difference between the mass-balance rates of the three surveyed glaciers from DEM differencing and the mass-balance rates calculated by sampling the map of elevation differences at the location of the laser altimetry measurements only (last two columns of Table 1). The standard deviation of this difference (0.05 m w.e. a^–1) is used as the profile-to-glacier error.

The regional extrapolation error (E _EXT) is obtained from the satellite and USGS DEM differencing. E _EXT is computed from the difference between the area-weighted mass-balance rates of the three profiled glaciers and the mass-balance rate of the unsurveyed area (˜2997 km²). The area-weighted mass-balance rate of the three profiled glaciers is more negative than the mass-balance rate of the unsurveyed area, and the difference, 0.15 m w.e. a^–1, is used as the regional extrapolation error.

The profile-to-glacier error and the extrapolation error are then added to the random errors in Eqn (3) for total error (E _tot).

(3)

Systematic errors in this study include the USGS map errors and the extrapolation errors.

Error estimates from the DEM differencing method

The main sources of uncertainties for the DEM differencing method are the vertical random errors in the USGS (15–45 m; Aðalgeirsdóttir and others, 1998), ASTER (±15 m; Fujisada and others, 2005) and SPOT5 (±10 m; Korona and others, 2009) DEMs and errors in USGS map dates (±3.5 years) (Table 1; see also Berthier and others, 2010). These different elevation errors were summed quadratically and divided by the square root of the number of map elevation contours to obtain the total error in the elevation changes. The 1σ vertical errors for the DEMs from the literature are confirmed by their comparison with ICESat data. After correction of the vertical biases with ICESat, the 1σ error is ˜8 m for the SPOT5 DEM and 11–17 m for the ASTER DEMs (Table 1). The quality of the adjustment of the 2006 SPOT5 DEM using ICESat data and the consistency of this vertically adjusted SPOT5 DEM with the 2007 laser altimetry measurements was further checked in the upper reaches of Nabesna Glacier. In less than 1 year (between September 2006 and June 2007), only small elevation differences are expected. The mean difference is only –0.4 m, with a standard deviation of 5.9 m (N = 38 228). For the USGS DEM, on the ice-free terrain, the 1 σ error (also estimated using ICESat data) is 25 m, an intermediate value between the published error in the ablation (15 m) and accumulation areas (±45 m) of glaciers. By adjusting the old and new elevation dataset (USGS and satellite DEMs) to a common altimetric reference (ICESat laser profiles) on the surrounding ice-free terrain (Nuth and Ka¨a¨b, 2011), we minimized systematic elevation errors due to poor geodetic control (Berthier and others, 2010, supplementary table S2). For unsurveyed areas, only 11% in the Wrangell Mountains, we assumed that errors doubled those calculated on surveyed areas. A ±10% error was also included for the total ice-covered area. Further details of our error analysis are provided by Berthier and others (2010, supplementary text).

Results

Elevation changes and specific mass balances of three profiled glaciers

The specific mass-balance rate is negative for all investigated glaciers and time periods (Table 2), ranging from –0.07 ± 0.09 m w.e. a^–1 (Nabesna 1957–2000) to –0.46± 0.05 m w.e. a^–1 (Kennicott 2000–07). Both glaciers with repeat altimetry measurements in 2000 and 2007 (Nabesna and Kennicott) show strong acceleration in mass loss.

Table 2. Glacier area and mass-balance rate for the three profiled glaciers as well as the entire glacierized area of the Wrangell Mountains for various time periods. Mass-balance rates are derived from three methods including center-line laser altimetry, sequential DEM analysis and elevation differences obtained from the sequential DEMs sampled along the altimetry profiles (‘Profile only’ column). Laser altimetry results are derived from two extrapolation methods (A and B; see text). ‘Profile only’ indicates that the map of elevation differences has been sampled at the location of the laser altimetry measurements and the elevation differences are then processed using the same extrapolation techniques as the laser altimetry data

Elevation changes of Nabesna Glacier show increased thinning rates for surface elevations below 2200 m a.s.l. for the more recent period 2000–07 (Fig. 2c–e). The specific mass-balance rate of Nabesna is about four times more negative in 2000–07 (–0.29 ± 0.06 m w.e. a^–1) than in 1957–2000 (–0.07 ± 0.09 m w.e. a^–1) (Table 2). Using Landsat images, we find a ˜3% area shrinkage between 1957 and 2010 (Table 2) while the glacier front has retreated by ˜0.5 km.

The mass balance of Kennicott Glacier has also become more negative, with rates dropping from –0.17 ± 0.06 m w.e. a^–1 during 1957–2000 to –0.46± 0.05 m w.e. a^–1 during 2000–07 (Fig. 2d–f). The area reduction was ˜2.5% and the terminus receded by an average of 0.15 km during 1957–2010.

Nizina is the only glacier that has two nearly parallel altimetry profiles along its trunk below˜ 1200 m a.s.l. The measurements sampled almost all elevation bands up to the highest elevations (Fig. 2g). The average thinning rates of the two profiles below 1200 m are very similar at –1.45 and –1.46 m a^–1, and the standard deviation of elevation change for this lower part of the two profiles below 1200 m is 0.14 m a^–1. The final ∆z curve is obtained by taking an average of these two profiles. Of the three surveyed glaciers, Nizina Glacier exhibits the most negative mass-balance rate (–0.44 ± 0.06 m w.e. a^–1) for the 1957–2007 period (Table 2). As this glacier was profiled for the first time in 2007, the rate of mass loss between 2000 and 2007 could not be assessed. This glacier had the highest shrinkage in area (˜4%) of the three profiled glaciers. The terminus position retreated by an average of 2.8 km from 1957 to 2010. Sequential DEM analysis also confirmed that, of the three glaciers surveyed by laser altimetry, Nizina had the largest mass loss during 1957–2005 (Fig. 4).

Fig. 4. Thickness change measurements from DEM differencing during the period 1957–2005. The colored boxes indicate the dates of the different satellite DEMs. Date format is month-day-year. The laser altimetry profiles from 2007 are indicated by black lines. Data gaps are indicated in white.

Mass-balance rate of the entire Wrangell Mountains glacier complex

The regional mass-balance rates derived using the normalization technique (method B) are in good agreement with those derived using the elevation-dependent technique (method A) (Table 2). Both methods show accelerated mass loss over the Wrangell Mountains during 2000–07, compared to 1957–2000 (Fig. 3), although the acceleration is less significant than for Nabesna and Kennicott Glaciers taken separately because of the large regional extrapolation errors. The area at lower elevations is very well represented by the three surveyed glaciers (Fig. 3d). Averaging the regional mass-balance rates from both methods, we find –0.07 ± 0.19 m w.e. a^–1 for 1957–2000 and –0.24±0.16 m w.e. a^–1 for 2000–07, and hence a 3.5-fold increase in the rate of mass loss (Table 2). The total area change of the Wrangell Mountains glaciers is ˜ 4% during 1957–2010.

Compared with 1957–2000, the mass-balance rate is more negative in 1957–2007. This is due to the inclusion of data for 2000–07, a period of rapid thinning, and also due to the addition of the measured mass balance for Nizina Glacier which has a more negative multi-decadal mass-balance rate than Nabesna and Kennicott Glaciers. The first effect dominates as the inclusion/exclusion has a very small influence on the regional mass-balance rate. The 1957– 2007 regional mass-balance rate evolves from –0.17±0.18 m w.e. a^–1 to –0.15 ± 0.18 m w.e. a^–1 when the measured profile of Nizina Glacier is excluded, a difference that is well within the error bounds.

Our mass-balance rates for all three time periods are calculated using the hypsometry derived from the 1957 USGS DEM. In order to test the validity of our approach, we recalculated the mass-balance rates using the hypsometry derived from the global DEM version 2 (GDEM v2) derived from ASTER stereo images acquired between 2000 and 2011. The mass-balance rates of Nabesna, Kennicott and Nizina glaciers for 1957–2007 were within 0.01 m w.e. a^–1 of the rates presented in Table 2.

Comparison with sequential DEM analysis

We compare results derived from the laser altimetry data with those from sequential DEM analysis (Fig. 5; Table 2). It should be noted that both laser altimetry and the sequential DEM analysis use the same USGS maps and hence are not independent estimates. The regional specific mass-balance rates during 1957–2007, when all three profiled glaciers were included in the analysis, are in good agreement with those derived from satellite DEM differencing during 1957–2005 (Table 2).

Fig. 5. Elevation changes per 30.48 m elevation bin from laser altimetry (1957–2007) used for computing glacier-wide balances (Fig. 3) compared to those from DEM differencing (1957–2005) sampled along the same profiles and from DEM differencing (averaging the elevation changes over all gridcells of each elevation bin). The area–altitude distribution is shown for the 30.48 m bins.

We use the map of elevation difference obtained by DEM differencing to detect any systematic bias due to center-line sampling of the glaciers by laser altimetry. Following Berthier and others, (2010), the map of elevation differences from DEMs is sampled at the location of the 2007 laser altimetry profiles to identify whether the laser altimetry sampling causes a systematic bias in the glacier-wide mass-balance rate (Fig. 5). The resulting mass-balance rate is almost unchanged for Kennicott, more negative for Nizina and less negative for Nabesna (Table 2, ‘Profile only’ column), but the differences are well within error limits. Hence, for these three glaciers, there is no evidence of a systematic error in glacier-wide balances caused by the assumption that elevation changes per elevation band along the flight lines are representative of the entire elevation band.

We compare the mean elevation change profiles obtained from laser altimetry during 1957–2007 to those derived from DEM differencing (Fig. 5). For elevations below 3000 m a.s.l. the elevation changes derived from the two techniques compare well and the differences are within error bounds for most cases. This comparison is especially important in the upper elevations where the laser altimetry data analysis involves some degree of subjectivity due to the large scatter of individual laser shots and data gaps. In the upper reaches of the glaciers (typically above 3000 m a.s.l.), elevation changes from laser altimetry (extrapolated or assumed to be zero) often deviate from the elevation changes from the DEM differencing method. This is especially the case on Kennicott Glacier (Fig. 5c), where DEM differencing suggests an unexpectedly strong thinning above 3000 m a.s.l. There are no laser altimetry data to determine the accuracy of the satellite DEM in this area and to assess whether this thinning is real. The precision of the USGS and ASTER DEMs is notoriously low in these flat and textureless regions. Fortunately, only 10% of the ice-covered area of the Wrangell Mountains is found above 3000 m a.s.l., and <5% above 3500 m a.s.l., so the glacier-wide and regional mass-balance rates show little sensitivity to the inclusion/exclusion of these uncertain areas.

Discussion

Our detailed analysis of elevation change using both laser altimetry and DEM differencing shows the validity of assuming that center-line elevation changes represent the entire width of the glacier for our study area. Mass-balance rates of the Wrangell Mountains derived from center-line sampling do not exhibit a systematic bias. This is probably because the main tributaries of the three glaciers have been surveyed (Fig. 1). Secondly, an updated glacier outline is available for 2010 to take into account the (relatively modest) change in glacierized area. These two reasons, among others, were cited by Berthier and others, (2010) to explain the overestimate of mass loss by Arendt and others, (2002). A single center-line profile will not adequately represent a glacier if the other major tributaries are not surveyed and the area change is not accounted for. Hence, large differences in potential bias and even its sign can occur between glaciers as illustrated in Berthier and others (2010, supplementary table S4). Similar conclusions were drawn recently for Glacier Bay, southeast Alaska (Johnson and others, 2013). Our study also indicates that, for the Wrangell Mountains, the largest source of uncertainty in the laser altimetry estimate is the regional extrapolation errors, i.e. the need to assume that the three profiled glaciers are representative of the whole region. This is also in agreement with Johnson and others, (2013) who showed that this error can be reduced by sampling a larger number of glaciers.

Our approach of using only the USGS hypsometry for calculating mass-balance estimates over all time periods is validated using GDEM v2 hypsometry. The differences in the mass-balance rate for the three profiled glaciers using the new hypsometry are within 0.01 m w.e. a^–1 of those calculated using the USGS DEM (Table 2). Thus our approach of using the hypsometry from the old USGS maps works well for a region with relatively modest thinning and retreat, like the Wrangell Mountains.

Our estimate of mass-balance rate for Nabesna, –0.07±0.09 m w.e. a^–1 during 1957–2000, is less negative than the –0.17 ± 0.07 m w.e. a^–1 for the same period estimated by Arendt and others, (2002). The difference between the two estimates is mainly due to the different methods used to derive the mean ∆z curve in the two studies. Arendt and others use only those laser altimetry points that cross a USGS map contour. This method may result in unrealistic elevation change values or data gaps when the laser returns are from a crevasse or steep ice surface. At upper elevations, Arendt and others (2002, supplementary fig. S1) determined a linear relation between elevation change and elevation for measured regions. They then applied that same function to elevations above those where measurements were available (between 2100 and 3100 m a.s.l.) for a smooth transition to zero elevation change line. In our study, the last available laser altimetry elevation change point, at 2621 m a.s.l., sharply joins the zero-change line. The difference between these two estimates, based on the same laser altimetry measurements and the same USGS DEM, shows the range of errors due to the need to (1) fill data gaps (e.g. unsampled elevation bands) and (2) subjectively discard some data in the region of anomalous changes in both studies.

The mass-balance rate of the Wrangell Mountains is less negative than that found in many other regions of Alaska over comparable periods. Arendt and others, (2006) found an average mass-balance rate of –0.74 ± 0.1 m w.e. a^–1 for the Western Chugach Mountains during 1950/57 to 2001/04 that was largely dominated by Columbia Glacier, a large tide-water glacier. The maritime glaciers in the southeastern part of the St Elias Mountains (14 580 km²) including the rapidly thinning Yukatat Icefield (Trüssel and others, 2013) and the Coast Mountains (Glacier Bay, Juneau and Stikine Icefields) had a mass-balance rate of –1.10 ± 0.29 m w.e. a^–1 during 1948/87 and 2000 (Larsen and others, 2007). These glacierized regions of southeast Alaska and northern British Columbia include tidewater, land- and lake-terminating glaciers, and the mass-balance rates are several times more negative than our regional mass-balance rates over the Wrangell Mountains for 1957–2000 or 1957–2007 (Table 1). More negative mass-balance rates were also found for 107 glaciers (42 km²) in the Brooks Range, in the continental climate of the north of Alaska (–0.54 ± 0.05 m w.e. a^–1), for the shorter 1970–2001 period (Geck and others, 2013). During the period 1962–2006, the regional mass-balance rate of all Alaskan glaciers was –0.48 ± 0.10 m w.e. a^–1 (Berthier and others, 2010) compared to –0.17 ± 0.18 m w.e. a^–1 calculated for the Wrangell Mountains in this study.

Conclusions

In this study, we combined airborne laser altimetry, topographic maps from the 1950s and recent DEMs from satellites to study five decades of mass change on the Wrangell Mountains. The regional glacier mass balance of the Wrangell Mountains, where the vast majority of glaciers are land-terminating, is several times less negative than that of other glacierized mountain ranges that include tidewater, lake- and land-terminating glaciers, such as the St Elias, Coast or western Chugach Mountains near the Gulf of Alaska.

One of our extrapolation techniques normalizes the elevation range of the glaciers so that the typically observed pattern of largest thinning rates close to the terminus is preserved. Although this approach is more physical, we find that for the Wrangell Mountains the results are similar to those derived from the widely used approach that assumes an elevation-change/elevation relation for the unmeasured glaciers despite large differences in glacier termini.

We observe no center-line bias for the laser altimetry analysis of the three profiled glaciers in the Wrangell Mountains after comparing with satellite DEM differencing. This is mainly attributed to the facts that (1) an updated glacier outline is available close to the year of the altimetry survey and (2) the different tributaries of the three surveyed glaciers are well sampled by the laser altimetry. The largest source of uncertainty in the regional mass-balance rate is the regional extrapolation error. This error stems from the need to assume that the three surveyed glaciers, which cover˜39% of the total ice area, are representative of all other glaciers in the mountain range.

The glaciers in the Wrangell Mountains have lost mass at an accelerated rate over the past decade as revealed by laser altimetry. The central value of the uncertainty range of the 2000–07 mass loss estimated for the Wrangell Mountains glaciers (–0.24 ± 0.16 m w.e. a^–1) exceeds the central value of the uncertainty range of the 1957–2000 period (–0.07 ± 0.19 m w.e. a^–1) by more than a factor of 3. The acceleration of mass loss is consistent with general patterns of temperature increase in Alaska (Arendt and others, 2009); however, the exact drivers for the observed changes in the Wrangell Mountains require further investigation.