Mass change rate – uncertainties

Sources of uncertainty for our measurements are the uncertainty on the elevations incoporated into the regression, uncertainty on the ELA, the effect of varying the maximum deviation allowed from the first elevation, different density scenarios, and uncertainty of the penetration depth of the C-band SRTM DEM. These are covered in turn below.

The uncertainty associated with the for each pixel is calculated from the model covariance matrix (e.g. Reference Aster, Borchers and ThurburAster and others, 2005), which accounts for the uncertainties on the elevations incorporated into the regression. The 95% confidence interval for the volume change rate uncertainty is calculated using the formula (e.g. Reference Howat, Smith, Joughin and ScambosHowat and others, 2008a). U is the total ‘volume’ of uncertainty, calculated by taking the uncertainty at each pixel, multiplying it by the area of the pixel (to determine a ‘volume’ of uncertainty for that pixel) and then adding together ‘volume’ of uncertainty for each pixel where is calculated. N is the number of independent pixels (e.g. Reference Howat, Smith, Joughin and ScambosHowat and others, 2008a), which we determine by dividing the total area by the area over which off-ice are correlated (e.g. Reference Rolstad, Haug and DenbyRolstad and others, 2009). We estimate the scale at which the are independent by finding the area at which the variance of the off-ice begins to ‘flatten’ (for details of the method see Reference Rolstad, Haug and DenbyRolstad and others, 2009; Reference Willis, Melkonian, Pritchard and RamageWillis and others, 2012b), which we estimate to be 810 m ×810 m (see Reference Willis, Melkonian, Pritchard and RiveraWillis and others, 2012a, Reference Willis, Melkonian, Pritchard and Ramageb; Reference Melkonian, Willis, Pritchard, Rivera, Bown and BernsteinMelkonian and others, 2013, for examples). The total contribution from the uncertainty of individual elevations is ~0.05 Gt a^–1.

We use ELAs of 925, 1100 and 1050 m for Taku, Menden-hall and Lemon Creek glaciers, respectively (Reference PeltoPelto, 2000; Reference Motyka, Neel, Connor and EchelmeyerMotyka and others, 2003; Reference Boyce, Motyka and TrufferBoyce and others, 2007; Reference PeltoPelto and others, 2008). A regional ELA of 1000 m (e.g. Reference Larsen, Motyka, Arendt, Echelmeyer and GeisslerLarsen and others, 2007) is used for the rest of the icefield. The choice of ELA determines the accumulation and ablation zone for each glacier basin, which affects our mass change rate in two ways. First, it determines the percentage of the calculated using +4.4/–10 m a⁻¹ versus +4.4/–5 m a⁻¹ allowed deviation from the first elevation. Second, it determines the percentage of the for which a density of 700 kg m ^–3 vs 900 kg m ^–3 is assumed. We assume an ELA uncertainty of ±200 m, and difference the mass change rates at the upper and lower bounds. This difference (0.18 Gt a^–1) is taken as the uncertainty due to ELA for our mass change rate.

The deviation allowed from the first elevation in the time series has a large impact on the mass loss rate (e.g. Reference Willis, Melkonian, Pritchard and RiveraWillis and others, 2012a; Reference Melkonian, Willis, Pritchard, Rivera, Bown and BernsteinMelkonian and others, 2013). The unsymmetric cut-off (+4.4/–5 m a⁻¹ in the accumulation zone, +4.4/–10 m a^–1 in the ablation zone) may bias our results towards thinning, but an unsymmetric cut-off is more physically justified than a symmetric cut-off (Reference Melkonian, Willis, Pritchard, Rivera, Bown and BernsteinMelkonian and others, 2013). One year of accumulated precipitation undergoes densification into firn and eventually ice, so +4.4 m a^–1 is chosen as a reasonable upper limit on the maximum thickening expected over the large areas covered in this study. We note that a point measurement may yield a higher than +4.4 m a ^–1 but we would not expect a single point to be representative of sustained thickening rates over many hundreds of square kilometers. Figure 4 highlights the importance of applying this cut-off, showing that cloud-influenced elevations would seriously degrade the quality of our regression for each pixel. Tests using a higher positive cut-off produced discontinuous and incoherent ‘splotches’ of extreme positive that are unrealistic.

Fig. 4. Elevation values and for randomly selected pixels over Field Glacier. The leftmost elevation in each graph is the SRTM elevation at that pixel. Blue lines indicate the calculated for each pixel. Elevation points bolded red are excluded from calculation. The bottom right panel shows the map (location in Fig. 5), with numbered circles indicating the location corresponding to each graph.

To assess the uncertainty due to the choice of maximum deviation allowed, we find the mode of the distribution of elevation differences between all ASTER DEMs and the SRTM DEM, normalized by dividing by the time interval between each ASTER DEM and the SRTM DEM (Reference Melkonian, Willis, Pritchard, Rivera, Bown and BernsteinMelkonian and others, 2013). The mode is not affected by the choice of allowed deviation, and so is an independent measure to which we can compare our regression-derived rates. We are not suggesting the mode rate represents the best estimate of the average . It is, rather, approximately what the rate would be if the distribution were Gaussian around the peak . The ‘mode rate’ is –0.10 m w.e. a ^–1. We incorporate the difference between the mode rate and our regression rate (0.08 Gt a ^–1) that results from using the same density and penetration assumptions as our regression rate to our uncertainties to fully account for any possible bias.

We assume volume change in the ablation zone is at a density of 900 kg m ^–3 (e.g. Reference McGeeMcGee, 1995; Reference Rignot, Rivera and CasassaRignot and others, 2003; Reference Berthier, Arnaud, Baratoux, Vincent and RémyBerthier and others, 2004, 2009; Reference Truffer, Motyka, Hekkers, Howat and KingTruffer and others, 2009) and volume change in the accumulation zone has a density of 700 kg m ^–3 (e.g. Reference ZempZemp and others, 2010; Reference GardnerGardner and others, 2013). However, the density of the lost material is likely variable. We take the difference between the mass change rate obtained using our density scenario and the rate produced by assuming all mass change is at a density of 900 kg m ^–3 as the uncertainty due to density (0.13 Gt a⁻¹), similar to Reference Ka¨a¨b, Berthier, Nuth, Gardelle and ArnaudKa¨a¨b and others (2012).

The penetration depth of the C-band SRTM into snow and ice is an additional source of uncertainty. Reference Rignot, Echelmeyer and KrabillRignot and others (2001) cited uncertainties of ±2 m for C-band penetration on temperate ice in Alaska; we take 0 m penetration (the minimum possible penetration) and 4 m penetration to be end-member cases and examine the effect on our mass change rate. Assuming 0 m penetration depth reduces the mass loss rate by +0.38 Gt a^–1, assuming 4 m of penetration increases the mass loss rate by –0.38 Gt a^–1. We therefore add 0.38 Gt a ^–1 to our uncertainties to account for penetration depth uncertainty. Our consideration of the penetration depth uncertainty is conservative, as we have already tried to correct for the penetration effect.

Taking the square root of the sum of the squared uncertainties gives us an overall uncertainty of 0.46 Gt a^–1 (equivalent to an area-averaged mass balance of 0.12 m w.e. a ^–1). Uncertainties beyond the uncertainty from the regression itself are added as an average rate to the uncertainties for individual glacier basins in Table 1.

ASTER (optical)

Feature- or pixel-tracking has been used to track glacier speeds at numerous locations (e.g. Reference Scambos, Dutkiewicz, Wilson and BindschadlerScambos and others, 1992; Reference Stearns and HamiltonStearns and Hamilton, 2005; Reference Howat, Joughin, Fahnestock, Smith and ScambosHowat and others, 2008b; Reference Willis, Melkonian, Pritchard and RamageWillis and others, 2012b).

To prepare the images we first use NASA’s Automatic Co-Registration and Orthorectification Package (AROP) (Reference Gao, Masek and WolfeGao and others, 2009) to co-register ASTER level 1B images to a Landsat reference image and orthorectify them using the SRTM DEM (Reference Willis, Melkonian, Pritchard and RiveraWillis and others, 2012a). We then produce sub-pixel offsets using the ‘ampcor’ normalized amplitude cross-correlation routine from ROI_PAC (Reference Rosen, Hensley, Peltzer and SimonsRosen and others, 2004; Reference Willis, Melkonian, Pritchard and RamageWillis and others, 2012b). We produce offset measurements every 120 m (8 pixels) using a reference window (matrix used to calculate the cross-correlation coefficient) of 480 m ˟ 480 m (32 ˟32 pixels) (Reference Willis, Melkonian, Pritchard and RiveraWillis and others, 2012a; Reference Melkonian, Willis, Pritchard, Rivera, Bown and BernsteinMelkonian and others, 2013). We experimented with offset samplings of <120 m (which increases computation time), but this did not increase our effective resolution or fill gaps in our spatial coverage.

The results are post-filtered by excluding offsets with a signal-to-noise ratio (SNR, which is the peak cross-correlation coefficient divided by the average cross-correlation coefficient) below a manually selected threshold. The threshold is determined by finding the maximum SNR that does not exclude coherent glacier velocities. An elevation- dependent correction (determined from ‘bedrock’ velocities) is applied to the speeds to correct for the elevation-dependent bias due to imprecise co-registration/orthorectification (Reference Ahn and HowatAhn and Howat, 2011). Applying an elevation-dependent correction significantly improves the contrast between bedrock offsets and glacier offsets (e.g. Reference Melkonian, Willis, Pritchard, Rivera, Bown and BernsteinMelkonian and others, 2013).

Each ASTER pair covers multiple glacier basins. One or both of the images in a pair may be cloudy over some basins. We esimate the uncertainty for each pair from motion on ice-adjacent ‘bedrock’ (Reference Willis, Melkonian, Pritchard and RamageWillis and others, 2012b; Reference Melkonian, Willis, Pritchard, Rivera, Bown and BernsteinMelkonian and others, 2013), which should be zero. Ice-adjacent ‘bedrock’ is defined as non-water, off-ice area within 5 km of a glacier, excluding the closest 600 m (5 pixels) to ensure glacier motion does not contaminate the ‘bedrock’. We calculate ice-adjacent bedrock speeds for each basin within the results from an individual image pair. We clip the ice-adjacent bedrock speeds within each basin at the mean plus or minus two standard deviations, then iteratively recalculate the mean and standard deviation to provide an estimate of uncertainty due to poor co-registration and DEM errors. We take the basin with the minimum mean plus one standard deviation of ice-adjacent bedrock speeds to be representative of well-correlated bedrock for the pair and assign this as the pair’s overall uncertainty. We feel this is justified, as the textured, low surface slopes of the glaciers usually display much more coherent, slowly spatially varying motions than the steeper bedrock.

Velocity uncertainty from pixel-tracking is well described in the literature (e.g. Reference Scherler, Leprince and StreckerScherler and others, 2008; Reference Ahn and HowatAhn and Howat, 2011). Using the techniques suggested by other studies and further developed in our previously published material (outlined above), we estimate velocity uncertainties to be around 0.1–0.2 m d ^–1. The magnitude of our uncertainty is higher than that of some other studies (e.g. Reference Scherler, Leprince and StreckerScherler and others, 2008) that have also applied pixel-tracking to ASTER imagery to estimate glacier velocities. Reference Scherler, Leprince and StreckerScherler and others (2008) took the mean and standard deviation of off-ice, or ‘stable ground’ motion as a measure of uncertainty, with ‘stable ground’ being defined by them as all displacements from –10 m to +10 m. We apply the method of Reference Scherler, Leprince and StreckerScherler and others (2008) to one of our ASTER image pairs (4 October 2004 to 8 May 2005) to compare with the uncertainty obtained using our method. Our method yields an uncertainty of ±0.09 m d^–1, while the method of Reference Scherler, Leprince and StreckerScherler and others (2008) provides an uncertainty of ±0.03 m d⁻¹. Consistently lower uncertainties would be expected using the method of Reference Scherler, Leprince and StreckerScherler and others (2008), because cutting off displacements at 10 m is generally a more severe restriction than our 2 clipping of adjacent off-ice displacements. However, they noted that removing glacier displacements using glacier outlines, which is part of our method, would be preferable to assuming all displacements from –10 m to +10 m are stable ground. We therefore consider our approach to be a reasonable and conservative estimate of the uncertainty for coherent motion, i.e. motion that has a high SNR and is spatially consistent with neighboring motions.

ALOS (radar)

We use L-band SAR pixel-tracking (e.g. Reference RignotRignot, 2008; Reference Strozzi, Kouraev, Wiesmann, Sharov and WernerStrozzi and others, 2008; Reference Rignot, Mouginot and ScheuchlRignot and others, 2011; Reference Burgess, Forster and LarsenBurgess and others, 2013) from five 46 day, ascending ALOS pairs to produce ice velocities for the Juneau Icefield. The ALOS SAR images have an initial pixel resolution of 3.3m (azimuth) by 8.3 m (range). The offsets have an effective resolution of approximately 150 m (azimuth) by 200 m (range), based on the step size of 50 pixels (azimuth) by 25 pixels (range). SAR images cover a broader area than the optical images and are not limited by cloud cover, providing velocities at many glaciers with no ASTER observations. SAR pixel-tracking also performs well in the snow-covered, high-altitude accumulation zone where optical images lack trackable features.

Raw ALOS SAR data are processed using ROI_PAC, and offsets are produced by ‘ampcor’. The results are SNR-filtered and run through the elevation-dependent correction routine. The SNR threshold we use for ALOS is less than ASTER; it depends on various aspects of the processing (chip size, sensor type) and is adjusted so that coherent areas of offset over glaciers are not excluded. ALOS interferometry does not yield much velocity information due to the relatively large motions, long separation between scenes, and changing surface characteristics.

Parallel flow (ERS)

We use ROI_PAC to process 1 day SAR image pairs taken by ERS to provide additional velocity measurements over the southern part of the JIF, including Taku Glacier. An ascending track interferogram from 28 to 29 October 1995 and a descending track interferogram from 29 to 30 October 1995 are unwrapped at a resolution of 41 m × 46 m, then georeferenced and down-sampled to a resolution of 300 m × 300 m. The interferograms are unwrapped using the branch-cut method (Reference Goldstein and WernerGoldstein and others, 1988) and power spectrum filter (Reference Goldstein and WernerGoldstein and Werner, 1998). Combining the results of these two unwrapped interferograms and applying the surface parallel-flow assumption – that all motion is parallel to the surface gradient of the glacier, obtained from the SRTM DEM – yields east–west and north–south velocities. Reference Joughin, Kwok and FahnestockJoughin and others (1998) provided the equations used to convert the unwrapped interferograms to velocities.

GPS

Transverse and longitudinal profiles of velocities and elevations on the JIF have been collected using rapid-static and real-time kinematic GPS surveying from 1993 through 2007 (e.g. Reference McGeeMcGee, 1995; Reference LangLang, 1997; Reference McGeeMcGee and others, 2007) and made available on the web (http://www.crevassezone.org). Several profiles are surveyed at least once per year to monitor year-to-year elevation change, and ~1000 stakes have been surveyed twice in a given year (within 7–14 days) to measure velocities (Reference McGeeMcGee and others, 2007; Reference McGeeMcGee, 2009). These surveys are always performed during July and August (Reference McGeeMcGee and others, 2007; Reference McGeeMcGee, 2009). The GPS displacements are estimated to have horizontal and vertical accuracy of 1 and 5 cm, respectively (Reference McGeeMcGee, 1995; Reference McGeeMcGee and others, 2007).

Taku Glacier

Taku Glacier has an area-averaged mass balance of +0.44 ±0.15 m w.e. a^–1 for 2000–10, and its continued advance corroborates the positive we find at its front. It contributes +0.34 ±0.12 Gt a ^–1 to the total mass change rate for the JIF.

Mendenhall Glacier

We find an area-averaged mass balance of –0.48 ±0.19 m we. a⁻¹ for Mendenhall Glacier from 2000 to 2013. Peak thinning rates of 9 m a⁻¹ are observed between 2000 and 2013 near the calving front of Mendenhall Glacier. These rates are comparable to the 1995–2000 thinning rate of 9ma^–1 in the terminus region produced by differencing center-line airborne laser profiles in 1995, 1999 and 2000 (Reference Motyka, Neel, Connor and EchelmeyerMotyka and others, 2003).

Juneau Icefield versus other icefields

The JIF is closer to being in a state of mass balance (–0.13 ±0.12 m w.e. a ^–1) than any other Alaskan icefield. GRACE mass loss estimates for glacierized Alaska given by Reference GardnerGardner and others (2013) (for 2003–09, area 87 100 km²) and Reference ArendtArendt and others (2013) (for 2003–10, area 82 505 km²) are equivalent to area-averaged mass balances of –0.61 ±0.21 and –0.79 ±0.13 m w.e. a ^–1, respectively.The longer-term mass balances for Alaska in general, and southeast Alaska in particular, are also more negative than for the JIF. Reference Berthier, Schiefer, Clarke, Menounos and RémyBerthier and others (2010) found an area-averaged mass balance of –0.48 0.10 m w.e. a –¹ for Alaskan glaciers from 1962 to 2006 (area 87 862 km²), and Reference Larsen, Motyka, Arendt, Echelmeyer and GeisslerLarsen and others (2007) found an area-averaged mass balance of –1.03 ±0.27 m w.e. a ^–1 from 1948/1961/1982/1987 to 2000 for glaciers in southeast Alaska (area 14 580 km²).

Comparison with DEM-differencing results for the Juneau Icefield

Comparison of our results with those of Reference Berthier, Schiefer, Clarke, Menounos and RémyBerthier and others (2010) and Reference Larsen, Motyka, Arendt, Echelmeyer and GeisslerLarsen and others (2007) suggests the possibility that the thinning rate has slowed at the JIF over the last 10 years compared to earlier years. Reference Berthier, Schiefer, Clarke, Menounos and RémyBerthier and others (2010) estimated an area-averaged mass balance between 1964 and 2006 of –0.48 0.±14 m w.e. a –¹ for the JIF (personal communication from E. Berthier, 2013). This is more negative than our average (–0.13 0.12 m w.e. a ^–1). Part of the difference is due to the lack of coverage over Taku Glacier (a large positive contributor) in Reference Berthier, Schiefer, Clarke, Menounos and RémyBerthier and others (2010). Removing Taku Glacier from our results produces an area-averaged mass balance of –0.31 0.17 m w.e. a ^–1.While this is still more positive than the area-averaged mass balance found by Reference Berthier, Schiefer, Clarke, Menounos and RémyBerthier and others (2010), the difference is of the same magnitud as the uncertainties.

Reference Larsen, Motyka, Arendt, Echelmeyer and GeisslerLarsen and others (2007) estimated an area-averaged mass balance of –0.62 m w.e. a^–1 for the JIF from 1948/1982/1987 to 2000. Their estimate includes Taku Glacier, meaning its mass balance is significantly more negative than our estimate and likely more negative than that of Reference Berthier, Schiefer, Clarke, Menounos and RémyBerthier and others (2010), who again lack the positive contribution of Taku Glacier.

Spatially, the difference between our and those of Reference Larsen, Motyka, Arendt, Echelmeyer and GeisslerLarsen and others (2007) and Reference Berthier, Schiefer, Clarke, Menounos and RémyBerthier and others (2010) appears to be concentrated at higher elevations (Reference Larsen, Motyka, Arendt, Echelmeyer and GeisslerLarsen and others, 2007, Fig. 7; Reference Berthier, Schiefer, Clarke, Menounos and RémyBerthier and others, 2010, supplemental figure S1g), which would be expected if the difference is driven by accumulation. We discuss this in more depth below, where we examine the connections between decadal-scale climate variability and thinning at the JIF.

Taku Glacier

We find an area-averaged mass balance of 0.44 ±0.15 m w.e. a^–1 for Taku Glacier from 2000 to 2010. We compare this rate with those from other studies below, highlighting the decadal-scale variability of mass balance at Taku Glacier, and note the possible effect of differences in the spatial extent of measurements between studies.

Reference Pelto, Kavanaugh and McNeilPelto and others (2013) used surface mass-balance measurements collected by JIRP from 1946 to 2011 to estimate the average mass balance for Taku Glacier, with surface mass balance defined as the ‘difference between the net snow accumulation and net ablation over one hydro-logic year’. They found a 1946–88 area-averaged mass balance of ~0.4 m w.e. a⁻¹. Reference Larsen, Motyka, Arendt, Echelmeyer and GeisslerLarsen and others (2007) differenced DEMs over a somewhat longer time period, 1948–2000, and found an area-averaged mass balance of ~0.4 m w.e. a –¹. Both of these are close to our more recent, 2000–10 area-averaged mass balance.

The 1988–2006 results of Reference PeltoPelto and others (2008) show a shift to mass loss at Taku Glacier, with an area-averaged mass balance of –0.14 m a^–1. The anomalously large snowfalls at the JIF in late 2006 and early 2007 (e.g. Reference Luthcke, Arendt, Rowlands, McCarthy and LarsenLuthcke and others, 2008) noted above, which are not included in the measurements of Pelto and others (2008), account for a small part of the difference between our area-averaged mass balance and their 1988–2006 rate. These years are included in the surface mass-balance results of Reference Pelto, Kavanaugh and McNeilPelto and others (2013), who found a 2000–10 mass balance of –0.04 m a ^–1 for Taku Glacier, more positive than the 1988–2006 rate but still more negative than our rate of 0.44 0.15 m w.e. a ^–1.

A significant part of the difference is therefore due to factors other than the different time periods covered. This is also evident from calculating the average mass balance using the data in Reference Pelto, Kavanaugh and McNeilPelto and others (2013) for the study period of Reference Larsen, Motyka, Arendt, Echelmeyer and GeisslerLarsen and others (2007). Tabulating the mass-balance estimates of Reference Pelto, Kavanaugh and McNeilPelto and others (2013) for 1948–2000 yields an area-averaged mass balance of 0.26 m w.e. a⁻¹, somewhat less than Reference Larsen, Motyka, Arendt, Echelmeyer and GeisslerLarsen and others (2007).

Reference Pelto, Kavanaugh and McNeilPelto and others (2013) did not have any surface massbalance measurements in the ablation zone, where we find overall thickening, which may account for part of the difference between their rate and the rates found by our study and Reference Larsen, Motyka, Arendt, Echelmeyer and GeisslerLarsen and others (2007). Reference Pelto, Kavanaugh and McNeilPelto and others (2013) acknowledge this as ‘the principal (source of) error’ in their mass-balance record. The lowspatial density of surface massbalance measurements is another acknowledged source of error (Reference PeltoPelto and others, 2008, Reference Pelto, Kavanaugh and McNeil2013). However, we do not peculate in what direction this would bias their results.

Mendenhall Glacier

Reference Motyka, Neel, Connor and EchelmeyerMotyka and others (2003) found mass loss rates from 1948 to 1995 and 1995 to 1999 at Mendenhall Glacier that are two to three times our 2000–13 area-averaged mass balance of –0.48±0.19m w.e. a⁻¹. They estimated an area-averaged mass balance of −0.86±0.10 m w.e. a⁻¹ from 1948 to 1995 (by differencing a 1948 topographic map with a 1995 airborne laser profile) and –1.29±0.03 m w.e. a⁻¹ from 1995 to 1999 (repeat airborne laser profiles). In 1999–2000 they found a dramatic, positive shift to an areaaveraged mass balance of +0.85 ±0.11 m w.e. a ^–1. Reference Motyka, Neel, Connor and EchelmeyerMotyka and others (2003) attributed this to the anomalously cool and wet summer of 2000. The positive mass balance in 2000, along with high snowfall in late 2006 and early 2007, likely contributes to the more positive average mass balance we observe for Mendenhall from 2000 to 2013.

The surface mass-balance data available during our study period are sparse, but yield a more positive average mass balance than the preceding time periods given in Reference Motyka, Neel, Connor and EchelmeyerMotyka and others (2003). Reference Boyce, Motyka and TrufferBoyce and others (2007) extended the results of Reference Motyka, Neel, Connor and EchelmeyerMotyka and others (2003), estimating surface mass balances for 2003, 2004 and 2005 of –1.8, –1.2 and –0.9 m w.e. a ^–1, respectively. The average mass balance obtained by adding the values in Reference Boyce, Motyka and TrufferBoyce and others (2007) for 2000, 2003, 2004 and 2005 is –0.63 ±0.15 m w.e. a ^–1 assuming independent, annual uncertainties on mass balance of 0.3 m w.e. a^–1). Nevertheless, the effect of a positive balance year on the average surface mass balance is significant, and the high snowfall mentioned above for 2006–07 would have a similar effect to 2000 on the average for Mendenhall Glacier.

We spatially isolate the difference between our average and the rates of Reference Motyka, Neel, Connor and EchelmeyerMotyka and others (2003) by separating out our ‘lower glacier’ using the SRTM DEM (<800 m a.s.l.; Reference Motyka, Neel, Connor and EchelmeyerMotyka and others, 2003). Reference Motyka, Neel, Connor and EchelmeyerMotyka and others (2003) measured a lower glacier thinning rate of 1 m a^–1 between 1948 and 1982 and twice that between 1982 and 2000 (measured by differencing the 2000 center-line airborne laser profile with a center-line profile taken from 1982 topographic data). Our lower glacier thinning rate for 2000– 13 is 2.2 ± 0.3 m a ^–1. This is approximately the same as the 1982–2000 rate found by Reference Motyka, Neel, Connor and EchelmeyerMotyka and others (2003), therefore the difference between our 2000–13 and their 1982–2000 rate is concentrated at higher elevations.

Comparison of climate changes to thinning

We examine snowfall and summer temperature data from the Juneau International Airport station (provided on the US National Oceanic and Atmospheric Administration website, http://www.arh.noaa.gov/cliMap/akClimate.php) to determine if there is an association between changes in these climate variables and the various estimates of mass balance/ for the JIF given above. Figure 8 shows the tabulation of these data into average yearly snowfall (October–April) and average summer (June–August) temperature for study periods spanning 1948–2013 relative to 1970–2000 averages. The periods 1948–1982/2000 have the highest snowfall and lowest summer temperatures, 1982/1987–2000 and 1995– 99 have the lowest snowfall and highest summer temperatures, and 2000–2009/2010/2013 are intermediate, close to the 1970–2000 averages.

These trends are reflected in the average ELA of Taku Glacier, which we tabulate from data provided in Reference Pelto, Kavanaugh and McNeilPelto and others (2013). Taku’s average ELA is 917 m for 1948–82, 931 m for 1948–2000, 958 m for 1982–2000, 999 m for 1987–2000, 1089 m for 1995–99, and 949–959 m for 2000– 2009/2011. The higher thinning rates observed at Menden-hall from 1995 to 1999 by Reference Motyka, Neel, Connor and EchelmeyerMotyka and others (2003) are also consistent with temperature data that show this study period has the highest average summer temperature.

Higher average summer temperatures and lower average snowfall from 1982/1987–2000 versus 1948–2000 and 2000–2009/2010/2013 are likely also a factor in the higher thinning rates found by Reference Larsen, Motyka, Arendt, Echelmeyer and GeisslerLarsen and others (2007). Comparing Figure 2 and Figure 7a of Reference Larsen, Motyka, Arendt, Echelmeyer and GeisslerLarsen and others (2007), there is a clear jump in thinning rates for the JIF at the transition from the 1948 DEM to the 1982/1987 DEMs. Taku Glacier is largely covered by the 1948 DEM in their study, and they find an average similar to ours (~0.4 m w.e. a⁻¹).Llewellyn Glacier, covered by the 1987 and 1982 DEMs, is a clear example of the higher thinning rates they obtained using those more recent DEMs. We find an area-averaged mass balance for Llewellyn Glacier that is close to 0 (see Table 1), whereas Reference Larsen, Motyka, Arendt, Echelmeyer and GeisslerLarsen and others (2007) estimated a of roughly –0.7 km³ a⁻¹, equivalent to an area-averaged mass balance of –1.3 to –1.4 m w.e. a^–1. The difference in appears to be concentrated almost exclusively at higher elevations, where they found higher thinning than at the front of Llewellyn Glacier.

We examine the accumulation-area ratio (AAR) at Llewellyn Glacier and its sensitivity to the ELA in order to see if it is consistent with the lack of thinning that we observe at higher elevations. Figure 9 illustrates the effect of changing the ELA on the AAR for the six largest glaciers on the JIF (including Llewellyn), as well as Mendenhall Glacier. Only 1.5% of Llewellyn’s area is below 800 m, compared with 10–20% for the other glaciers. Assuming the ELA we use in our study (1000 m, also the average ELA used by Larsen and others 2007), the AAR of Llewellyn Glacier is 0.96, the highest of any of the six largest glaciers on the JIF. None of the others have an AAR above 0.9. Raising the ELA to 1400 m only decreases the AAR of Llewellyn Glacier to 0.82, the highest AAR assuming a 1400 m ELA of any of the glaciers shown in Figure 9. This suggests that rising temperatures would have a smaller impact at Llewellyn Glacier than at the other large glaciers on the JIF. Lower average snowfall from 1982/1987 to 2000 may therefore be a larger factor in high-elevation thinning at Llewellyn over that period relative to other glaciers.

Fig. 9. AAR vs elevation for the six largest glaciers on the JIF and Mendenhall Glacier. The green dots show AAR at the ELA we use, the blue dots show the AAR at 800 m, and the red dots show the AAR at 1400 m. Elevations are from the SRTM DEM.

A high AAR, relative insensitivity of AAR to ELA changes, and higher average snowfall and lower average summer temperatures from 2000 to 2009/2010/2013 versus 1982/1987–2000 are all consistent with the lower thinning rates we find at the higher elevations of Llewellyn Glacier compared with Reference Larsen, Motyka, Arendt, Echelmeyer and GeisslerLarsen and others (2007).

ASTER vs GPS speeds

Several observations of velocity using surveyed stakes and GPS spanning 1950–2006 were made at Taku Glacier along transverse Profile 4, near the equilibrium line (Reference McGeeMcGee and others, 2007; Reference PeltoPelto and others, 2008). No significant interannual variations in flow were observed despite changes to the mass balance of the glacier (Reference McGeeMcGee and others, 2007; Reference PeltoPelto and others, 2008). Pelto and others (2008) speculated that the velocity stability is due to the great thickness of Taku Glacier, which is up to 1500 m deep (Reference Nolan, Motyka, Echelmeyer and TrabantNolan and others, 1995). This thickness produces a high basal shear stress that limits basal sliding, which is more likely to be temporally variable (Reference Nolan, Motyka, Echelmeyer and TrabantNolan and others, 1995)). Reference Nolan, Motyka, Echelmeyer and TrabantNolan and others, 1995) and Reference PeltoPelto and others (2008) therefore concluded that motion of Taku Glacier is likely primarily due to internal deformation. However, seasonal speed variations, which are generally the result of changes in sliding (Reference PeltoPelto and others, 2008), are observed at Taku Glacier (Reference Nolan, Motyka, Echelmeyer and TrabantNolan and others, 1995); Reference Truffer, Motyka, Hekkers, Howat and KingTruffer and others, 2009). This suggests that sliding is a non-negligible component of motion at Taku Glacier, at least in the ablation zone.

Most of the JIRP GPS measurements were made in areas that are covered with snow for long periods during the year (Reference McGeeMcGee, 2009). The low-contrast snow cover makes acquiring reliable ASTER motions difficult. An exception is the JIRP ‘Longitudinal A’ profile, which runs down the center line of Taku Glacier and extends far enough into the ablation zone to overlap with reliable ASTER-derived speeds. Figure 7 compares GPS speeds from 20 to 23 July 2001 with ASTER speeds from image pairs covering October 2004 to May 2005, March 2009 to July 2009, and August 2009 to June 2010 and the 1995 ERS speeds (discussed below). Figure 6 shows the location of the points in Figure 7 as gray triangles.

The GPS speeds are consistently ~0.1 m d ¹ higher than the ASTER speeds. This is likely a seasonal effect, such as that observed by Reference Nolan, Motyka, Echelmeyer and TrabantNolan and others (1995) at their ‘Goat Ridge’ transect and Reference Truffer, Motyka, Hekkers, Howat and KingTruffer and others (2009) further down-glacier at the terminus. The GPS speeds shown in Figure 7 are from a short period (3 days) in late July, whereas the ASTER speeds are months-long averages and are from different years. The velocity stability noted above at Profile 4 is significantly up-glacier from these points, which lie along the trunk of Taku Glacier, well into the ablation zone.

Although the GPS has lower uncertainties, the ASTER results are more extensive, helping to fill in areas missed by GPS, such as the 14 km gap in GPS speeds for Taku Glacier between the lower end of the Longitudinal A profile and the terminus (Reference McGeeMcGee, 2009). ASTER and GPS speeds are complementary, as most ASTER measurements are in the ablation zone and most GPS measurements are in the accumulation zone. GPS speeds steadily increase down the Longitudinal A profile to a maximum of 1.14 m d^–1 at the point furthest down-glacier (Reference McGeeMcGee, 2009), ~10 km from the terminus. ASTER-derived speeds exhibit a maximum along the trunk of Taku Glacier in the bend north of Hole-in-the-Wall Glacier, below the maximum GPS speed of the Longitudinal A profile. Both the March 2009 to July 2009 and May 2010 to June 2010 pairs show this area of maximum speed, and decorrelation over the same area in the October 2004 to May 2005 pair could be due to relatively faster motion there.

Parallel-flow velocities, ALOS pixel-tracking vs GPS

The parallel-flow assumption applied to ERS interferometry extends our temporal coverage of speeds at Taku Glacier back to late October 1995, providing further information on the long-term variability in speeds at Taku Glacier. We check the velocities this produces against GPS velocities measured by the JIRP. The mean difference is 0.08 m d^–1, with a standard deviation of 0.08 m d^–1. Figure 10 gives a visual comparison of the GPS and parallel-flow velocities in map view. Figure 11 shows a one-to-one plot of GPS vs ERS velocities as gray triangles, and GPS vs ALOS velocities as red circles.

Fig. 10. Map of parallel-flow speeds calculated from one ascending and one descending ERS interferogram (from 28 October 1995 to 29 October 1995 and 29 October 1995 to 30 October 1995, respectively) and GPS velocities (red arrows) from JIRP. White arrows indicate ERS velocity at GPS points.

Fig. 11. One-to-one plot of GPS speeds (x –axis) vs ALOS (red circles) and ERS (gray triangles) speeds (y –axis). Dates are year/month.

The agreement between GPS and ALOS/ERS results illustrates the advantage of using radar to obtain reliable velocities in the accumulation zone, as this is where most GPS velocity measurements are made. ERS parallel-flow velocities are a better match with GPS in the accumulation zone than ALOS pixel-tracking. The time separation of the ALOS pairs leads to decorrelation, reducing the quality of velocities obtained in the accumulation zone. Interferometric synthetic aperture radar (InSAR)-derived velocities are an order of magnitude better than pixel-tracking for a given time separation (e.g. Reference Fialko, Simons and AgnewFialko and others, 2001), so the 1 day ERS InSAR can provide reliable results even over slower-moving areas.

Figure 7 shows the ERS speeds as well as ASTER and GPS for the lowest eight points of Longitudinal Profile A on Taku Glacier. The ERS speeds match the ASTER speeds, and are also consistently ~0.1 m d⁻¹ lower than the GPS speeds, again demonstrating the seasonality of speeds in the ablation zone of Taku Glacier, in contrast to the relative stability at Profile 4 near the ELA. The quality of the match between the ERS and GPS velocities and the match between ERS and ASTER speeds shown in Figure 7 corroborates the ASTER results, as ERS and ASTER are from two different data sources using two different techniques.

Interannual speed changes

We find no interannual acceleration at any of the outlet glaciers on the JIF. Figure 12 shows the lack of significant speed changes at Taku, Gilkey, Llewellyn, Meade and Field glaciers. At Taku Glacier, ASTER speeds from October 2004 to May 2005 closely match ALOS speeds from 2009 to 2011 along the first 6 km of the profile. The profiles at Meade Glacier show no interannual speed change from 2004 to 2009, although summer pairs from 2004 and 2009 indicate a possible seasonal effect. At Llewellyn Glacier there is a possible ‘shift’ in the speed profile from 1995 to 2004 by ~1 km, although with the large relative uncertainties this is not clear. The 1995, 2005 and 2009/2010 profiles all show a consistent top speed of 0.7–0.8 m d⁻¹. Speeds at Field Glacier are consistent, with the exception of small stretches of some high-uncertainty pairs, and no interannual or seasonal variations are discernible. For Gilkey Glacier, all pairs have low speeds at the front; the ALOS pairs from 2007– 08 and 2009–10 match the ASTER pair from June 2004 to August 2004. These are higher than the ASTER speeds from August 2004 to April 2005, and August 2009 to May 2010, starting at ~3.5 km along the profile, indicating possible seasonal variations in speeds at Gilkey. There is, however, no evidence of interannual acceleration, as ASTER and ALOS speeds separated by 5 years closely match. The absence of acceleration outside our uncertainties (~±0.1 m d⁻¹) indicates that the thinning observed at the front of Gilkey, Meade, Field and Llewellyn glaciers is not dynamic, and is instead primarily due to increasing ablation caused by rising temperatures in the region (Reference Stafford, Wendler and CurtisStafford and others, 2000).

Fig. 12. Longitudinal profile of speeds from velocity pairs covering Taku (a), Meade (b), Llewellyn (c), Field (d) and Gilkey (e) glaciers, starting at the front. Dates are year/month.

Taku Glacier – seasonal speed changes

Reference Truffer, Motyka, Hekkers, Howat and KingTruffer and others (2009) used GPS and repeat aerial photography to examine daily, seasonal and interannual speed changes at the terminus of Taku Glacier. Seasonal speed changes there are related to changing basal conditions (Reference Truffer, Motyka, Hekkers, Howat and KingTruffer and others, 2009). Figure 13 shows the SNR-filtered speeds from three ASTER image pairs, a ‘summer’ pair (March to July 2009), ‘non-summer’ pair (October 2004 to May 2005) and a ‘year-round’ average (August 2009 to June 2010). Tracks 1–3 are approximately in the same place as the tracks used by Reference Truffer, Motyka, Hekkers, Howat and KingTruffer and others (2009). Figure 14 shows speeds from three ASTER pairs along these tracks.

Fig. 13. Speeds from ASTER image pairs October 2004 to May 2005 (a), March 2009 to July 2009 (b) and August 2009 to June 2010 (c). The bottom left image shows the ablation zone of Taku Glacier and the location of images (a), (b) and (c). Speeds are higher in (b), which is the ‘summer’ pair.

Fig. 14. Speeds along tracks 1, 2 and 3 (see Fig. 13) for ASTER image pairs October 2004 to May 2005 (red), March 2009 to July 2009 (blue) and August 2009 to June 2010 (green). Uncertainties (μ ± σ) represented as gray bars. All speeds increase with distance from the front, ranging from about 0.05 to 0.7 m d⁻¹. ‘Summer’ speeds (March 2009 to July 2009) are higher than speeds for the other pairs, especially for track 3 (the southwestern track).

The ASTER speeds are approximately the same as Reference Truffer, Motyka, Hekkers, Howat and KingTruffer and others (2009), increasing from 0.1 m d⁻¹ near the front to 0.7 m d⁻¹ about 2000–2500 m upstream. Speeds from the ‘summer’ (blue) pair are 0.2–0.3 m d⁻¹ higher along track 3 than the non-summer and year-round speeds at distances >2000 m from the front. This result matches the findings of Reference Truffer, Motyka, Hekkers, Howat and KingTruffer and others (2009).

Summer temperatures are higher, leading to increased surface melting. Melt, in turn, may percolate to the glacier bed, increasing the pressure of the basal water system and reducing basal friction (Reference Iken and BindschadlerIken and Bindschadler, 1986) or altering the configuration of the subglacial drainage (Reference Howat, Tulaczyk, Waddington and BjörnssonHowat and others, 2008c). Reduced basal friction can lead to an increase in speed, accounting for higher summer speeds at the front of Taku Glacier. This has been suggested as a mechanism for velocity variations on Glaciar Soler, Northern Patagonia Icefield, by Reference Yamaguchi, Naruse, Matsumoto and OhnoYamaguchi and others (2003), who observed that higher speeds were associated with a higher liquid water input (either from increased melt or rainwater) and postulated that this reached the bed, reducing friction and increasing basal sliding (Reference Iken and BindschadlerIken and Bindschadler, 1986). Reduced basal friction from increased meltwater input to the base also leads to higher summer speeds at Mendenhall Glacier (Reference Motyka, Neel, Connor and EchelmeyerMotyka and others, 2003) and causes seasonal velocity variation at many other temperate glaciers (e.g. Reference Yamaguchi, Naruse, Matsumoto and OhnoYamaguchi and others, 2003; Reference Truffer, Motyka, Hekkers, Howat and KingTruffer and others, 2005; Reference Howat, Tulaczyk, Waddington and BjörnssonHowat and others, 2008c; Reference PeltoPelto and others, 2008; Reference SugiyamaSugiyama and others, 2011). This effect is less significant further up-glacier around the ELA (e.g. Profile 4), where greater thickness might contribute to higher basal friction and more stable velocities (Reference PeltoPelto and others, 2008).

Taku Glacier – flux

Table 2 gives the flux results at transverse Profile 4 on Taku Glacier for the ALOS pairs spanning 2007–11 and the ERS velocities from 1995. There is no significant variation in the flux between the ALOS pairs, which varies from (0.32–0.40) ±0.16 to (0.36–0.45) ±0.17 km³ a⁻¹ from 2007 to 2011, although the ERS flux from 1995 is somewhat greater at (0.51–0.64) 0.16 km³ a⁻¹ due to the higher velocities measured. Reference PeltoPelto and others (2008) use GPS velocities to measure a mean flux of 0.53 km³ a⁻¹ from 1996 to 2004, varying by 0.02 km³ a⁻¹. This value lies within the uncertainty for all of our flux estimates, indicating a stability in the flux at Profile 4 from 1995 to 2011. The precision of the remote-sensing-based techniques is not as great as GPS, but the consistency of the ALOS measurements suggests that the flux is indeed relatively stable.

Mendenhall Glacier – flux

Table 3 shows the fluxes at gates e and g (with g shifted 300 m north for our study) from Reference Motyka, Neel, Connor and EchelmeyerMotyka and others (2003) for the four ALOS pairs from 2007 to 2011. Our estimates do not show a significant change in flux as the glacier continues to thin and retreat. Reference Motyka, Neel, Connor and EchelmeyerMotyka and others (2003) do not give a value for flux at these gates; however, they do list a value of 0.08 km³ a⁻¹ for the ‘flux gate’ (Reference Motyka, Neel, Connor and EchelmeyerMotyka and others, 2003; Fig. 5d).

We take the flux at gate g to approximate calving flux, as the transect is very close to the glacier front. We shift this gate 300 m to the north to avoid profiling a section of Mendenhall Glacier that has already calved off by the time the ALOS imagery was acquired. With this in mind, we find that calving constitutes approximately 2.5 ±2% to 5 ±4% of total volume loss at Mendenhall Glacier, roughly the same as the 6% estimated by Reference Motyka, Neel, Connor and EchelmeyerMotyka and others (2003) and about the same as Reference Boyce, Motyka and TrufferBoyce and others (2007), who found that calving flux was ‘about 2.6% of ice lost by surface melting between 2002 and 2005’. This is significantly less than the percentage of mass loss due to calving at tidewater-cycle glaciers such as Marinelli Glacier on the Cordillera Darwin Icefield, Chile, where ~70% of the mass loss is from calving (Reference Melkonian, Willis, Pritchard, Rivera, Bown and BernsteinMelkonian and others, 2013), and is evidence that mass loss at Mendenhall Glacier continues to be primarily due to warming-induced thinning (e.g. Reference Motyka, Neel, Connor and EchelmeyerMotyka and others, 2003).

Flux at other glaciers on the JIF

All major outlet glaciers show a similar pattern and magnitude of velocities to those of Mendenhall, with higher speeds up-glacier and low front speeds (~0.1 m d⁻¹). Meade, Gilkey and many other glaciers on the JIF are lake-terminating, and tend to have much lower calving rates than tidewater-calving glaciers for equivalent depths (Reference Van der VeenVan der Veen, 2002). Low terminus speeds also suggest low calving rates (Reference Van der VeenVan der Veen, 2002), although without thickness measurements we cannot measure the calving flux to determine its magnitude relative to the overall mass change rate. There are certainly rapidly thinning lacustrine glaciers, but many of the most rapidly thinning (e.g. Glaciar Upsala, Southern Patagonia Icefield) calve into deep lakes and show a different pattern of glacier speeds, with speed rapidly increasing to a maximum at the front (e.g. Reference Skvarca, Raup and AngelisSkvarca and others, 2003; Reference Sakakibara, Sugiyama, Sawagaki, Marinsek and SkvarcaSakakibara and others, 2013). Back-of-the-envelope calculations show that the calving front of Meade and Gilkey Glaciers would have to be 300–500 m thick at the current front width and speed for the calving flux to account for even 10% of the current mass loss rate. Tabular icebergs indicate retreat driven by buoyancy-induced calving as the glacier recedes into deeper waters. Measuring the area loss at the front of Gilkey and Field Glaciers from ASTER imagery and assuming an average thickness of 100 m for ice lost yields a volume loss rate from retreat that is approximately 10% and 25% of the volume loss rate for those glaciers, respectively. Even an average depth of 100 m for the ice lost at the front seems too high for these glaciers; again, without measurements of the ice depth we cannot be certain. Overall, the lack of acceleration and low front speeds show that the pattern (e.g. the areas of maximum thinning behind the front of these glaciers) is primarily due to melting.