Lemon Creek Glacier

Lemon Creek Glacier (LCG; 11.6 km²; 58°23′N, 134°24′W) is located in the Juneau Icefield in the Coast Range of southeast Alaska (Fig. 1a), and is a temperate valley glacier (Reference Marcus, Chambers, Miller and LangMarcus and others, 1995). The B _a of the LCG has been monitored since 1953 by the Juneau Icefield Research Program (JIRP) (Reference Pelto and MillerPelto and Miller, 1990). LCG extends from 820 to 1400 m a.s.l. From the head of the glacier to the mean ELA at 1050–1100 m a.s.l. (1998–2010) (Reference Zemp, Nussbaumer, Gärtner-Roer, Hoelzle, Paul and HaeberliWGMS, 2011) (annual variations in ELA are illustrated in Fig. 2a), the glacier flows northward, and in the ablation zone it turns westward, terminating at ∼820 m a.s.l. LCG surface slope changes from ∼4° in the accumulation area to ∼18° at the termini. The glacier terminus retreated on average 10–13 m a⁻¹between 1998 and 2009. For LCG the observed B _a was on average −0.44 m w.e.a⁻¹ from 1953 to 2011 and −0.51 m w.e.a⁻¹ from 1996 to 2011 (Reference Pelto and MillerPelto and Miller, 1990; Reference Zemp, Nussbaumer, Gärtner-Roer, Hoelzle, Paul and HaeberliWGMS, 2011); winter and summer balances are not determined separately. LCG is located in a sub-Arctic region of temperate maritime climate, with annual precipitation of ∼3000–4000 mm and an average annual air temperature at the ELA of −1°C.

Fig. 2. (a) Area elevation band distribution for 2001 derived from Landsat ETM+ Mosaic and observed balance gradients of LCG for 1997/98, 2002/03, 2003/04, 2004/05, 2005/06, 2006/07, 2007/08, 2008/09, 2009/10 and 2010/11, the same years when Landsat TM imagery was obtained to estimate TSL trends (Table 2). (b) Area elevation band distribution for 1999 and 2011 derived from Landsat 7 ETM+ Mosaic and observed balance gradients of MG for 1998/99, 1999/2000, 2000/01, 2001/02, 2005/06, 2007/08 and 2011/12, the same years when Landsat TM/ETM+ imagery was obtained to estimate TSL trends (Table 3).

Mittivakkat Gletscher

Mittivakkat Gletscher (MG; 26.2 km² in 2011; 65°41′N, 37°48′W) is located in the Ammassalik region, southeast Greenland (Fig. 1b), and is a temperate glacier (Reference Knudsen and HasholtKnudsen and Hasholt, 1999). It extends from 160 to 880 m a.s.l. Since the end of the Little Ice Age around AD 1900, MG has undergone almost continuous retreat (Reference Knudsen, Nønberg, Yde, Hasholt and HeinemeierKnudsen and others, 2008; Reference MernildMernild and others, 2011a), in which the glacier area decreased by 18% (1986–2011) (Reference Mernild, Malmros, Yde and KnudsenMernild and others, 2012), volume decreased by 30% (1986–2011) (Reference MernildMernild and others, 2013) and the mean surface slope increased from 0.095 m m⁻¹ (=5.4°) to 0.104 m m⁻¹ (5.9°) (1986–2011).

For MG the B _a has been observed for 17 years since 1995/96, and the winter and summer balances individually for only 11 years (B _a was measured over the study area:17.3 km² in 1999 and 15.9 km² in 2011). B _a is −1.01 ± 10.74 m w.e. a⁻¹ (1995/96 to 2011/12), with a mean winter balance of 1.16 ± 0.20 m w.e.a⁻¹ and a mean summer balance of −1.99 ± 0.40 m w.e.a⁻¹ (1995/96 to 2001/02, 2004/05 to 2005/06, 2007/08 and 2011/12). In 2010/11, B _a was at a record setting, −2.45 m w.e.: about two standard deviations below the mean, and 0.29 m w.e. more negative than the previous observed record low B _a in 2009/10 (Reference Mernild, Knudsen, Hanna, Richter-Menge, Jeffries and OverlandMernild and others, 2011b). The loss of 1.63 m w.e. in 2011/12 was the fourth highest loss since 1995, and three of the four highest recorded B _a losses have occurred in the last three years included in this study. Since 1995, the mean ELA has risen from around 500 m a.s.l. to 750 m a.s.l. (Reference MernildMernild and others, 2011a updated). Figure 2b shows annual variations in ELA. MG is considered to be located in the Low Arctic (Reference Born and BöcherBorn and Böcher, 2001), and in a relatively wet and snowy part of Greenland (Reference Mernild and ListonMernild and Liston, 2010). An air temperature analysis reveals that the mean annual air temperature for MG was −2.2°C (1993–2011) at 515 m a.s.l. (Reference Mernild, Hansen, Jakobsen and HasholtMernild and others, 2008 updated), and a trend analysis of standard seasonal averages shows the following increases in seasonal air temperature for 1993–2011: 2.9°C in winter, 0.9°C in spring, 2.6°C in summer and 1.0°C in autumn (Reference Hanna, Mernild, Cappelen and SteffenHanna and others, 2012). The mean annual precipitation varied in the range 1400–1800 m m w.e.a⁻¹ (1998–2006) (Reference Mernild, Hansen, Jakobsen and HasholtMernild and others, 2008).

Satellite method at Lemon Creek Glacier

The TSL on LCG is readily identifiable on 34 scenes acquired in 1998 and 2003–11, and visualized using the US Geological Survey (USGS) Globalization Viewer software (Table 1). LCG falls in path/row 58/19 and 57/19; all images are false-color RGB (red, green, blue) composites, bands 3, 4 and 5, with a 2% linear stretch applied. The 7.5 min quadrangle digital elevation model (DEM) from the United States Geological Survey was used (further information: http://eros.usgs.gov/#/Guides/dem). The TSL is manually digitized for each scene. The exception is when the TSL rises to 1200 m a.s.l. or is <900 m a.s.l.: in both cases the surface slopes increase, leading to higher error margins. The satellite images were georeferenced in ArcMap 9.3 using five topographically unique reference points. The five ground control points (GCPs) are part of the benchmark survey network for the Juneau Icefield; their position is determined in the field using rapid static and real-time GPS equipment with an accuracy of 0.01 m horizontally and 0.05 m vertically. The registration errors between the Landsat 5 TM and 7 ETM+ products were 24 m root-mean-square error (RMSE) based on five GCPs. The image spatial resolution of 30 m and the registration of 24 m combined with mean surface gradients of 0.04–0.08 mm⁻¹ yields an error of ±1–4 m in TSL elevation, with a mean of 1.56 m. The data frame containing imagery and base map was transformed to North American Datum (NAD) 1983, Universal Transverse Mercator (UTM) zone 8N to ensure spatial accuracy for measurements.

Satellite method at Mittivakkat Gletscher

For MG the satellite imagery data were available through the USGS ‘EarthExplorer’ online database. The area of MG is covered by two Landsat overpasses path/row 231/14 and 232/14. The TSL on MG was identified using imagery from Landsat 5 TM and 7 ETM+ having a ground resolution of 30 m (Table 1). The TSL was manually digitized from the 26 scenes (Table 3, further below) by composing a false-color image from bands 2, 3 and 5, to maximize the snow-cover contrast in the image. A DEM was extracted from the Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) Global Digital Elevation Model Version 2 (GDEM v2), providing a vertical average precision of ∼12 m over Greenland (Reference TachikawaTachikawa and others, 2011). The vertical error is expected to be closer to the GDEM v2 standard ±8.7 m precision due to the gentle slope of the majority of the glacier from where the measurements were taken (Reference TachikawaTachikawa and others, 2011). The lateral error associated with GDEM v2 is a little more than half a pixel (17 m). ASTER GDEM v2 is a product of the US Ministry of Economy, Trade, and Industry and NASA. The overall registration errors between the Landsat 5 TM and 7 ETM+ products were measured to be 21 m RMSE based on 31 GCPs (first order). The differences between ASTER GDEM v2 and Landsat 5 TM were 22 m RMSE based on 25 GCPs, and between ASTER GDEM v2 and Landsat 7 ETM+ were 20 m RMSE (23 GCPs). The vertical error produced by the registration errors was found to be 1.8–3.1 m, averaging 2.2 m, with a spatial resolution of 30 m and mean surface gradients of 0.06–0.10 m m⁻¹ . The data for MG were projected in World Geodetic System 1984 (WGS84), UTM zone 24N. The accuracy of the Landsat imagery was validated by in situ GPS measurements taken in the field from several years, and all measurements were within half a pixel (15 m) (Reference Mernild, Malmros, Yde and KnudsenMernild and others, 2012).

Calculations

For LCG and MG the snow ablation rates were calculated based on both the TSL–mass-balance-gradient method and the snow-pit–satellite method (Fig. 3a and b). For the TSL–mass-balance-gradient method snow ablation rates were calculated from the rise in TSL (Tables 2 and 3), where the TSL migration rates were multiplied with the field-determined balance gradients near the TSL. For the snow-pit– satellite method, the snow ablation rates were calculated based on observed snow loss in snow pits (Tables 2 and 3), where the snow depths were divided by the time interval for the TSL to transect the snow pits. For example, if the snowpack depth on 1 July in a snow pit was 1.4 m w.e., and on 12 August the TSL reached the snow pit, then it took 42 days to ablate 1.4 m w.e. of snow, yielding the snow ablation rate. At LCG, snow-pit excavations were conducted for the years 1998 and 2003–12 (Table 2) and at MG for the years 1999, 2000–02, 2006, 2008 and 2012 (Table 3).

Fig. 3. Estimated snow ablation rates for LCG (a) and MG (b) based on the TSL–mass-balance-gradient method and the snow-pit– satellite method (Tables 2 and 3).

Table 2. Dates and elevation of TSL observations on LCG. The TSL elevation rate is the mean daily rise in TSL elevation since the previous observation date at least 15 days earlier. Snow-pit depth (SWE) is from the original mid-July excavation on the date indicated; the accumulated snow has completely ablated. Snow ablation rate is the ablation needed to remove the snow-pit accumulation by the observation date

Table 3. Dates, mean elevation and standard deviation of satellite-derived TSL on MG. Also shown are observed snow-pit and estimated snow ablation values. Where the TSL elevation rate is negative it indicates that the TSL has moved down-glacier. For the years 2000, 2001, 2002 and 2008, no snow ablation rates are estimated, either because of an insufficient number of snow pits or no available Landsat imagery for estimating TSL

The annual ELAs for LCG and MG were estimated based on second-order polynomial regression between the TSL elevation and the TSL date (remembering that for any specific date of TSL observation the transient mass balance at the TSL is zero), where ELA faces the highest calculated TSL at the end of the ablation season (Fig. 4a and b). In many years the time difference between available and usable imagery where the TSL is visible at the end of the ablation season can be several weeks, which is why the approach of utilizing a second-order polynomial regression is an advantage. The annual AAR was calculated from the estimated ELA. Both glaciers were partitioned into elevation bands (LCG 50 m elevation bands and MG 100 m elevation bands), and the AAR for a given observed year was determined based on the glacier area above and below the ELA: for LCG a fixed area of 11.6 km² was used for the entire period, and for MG an assumed linearly decreasing area from 17.3 km² (1999) to 15.9 km² (2011) was used.

Fig. 4. (a) Satellite-derived LCG TSL elevations throughout five summer periods: 2004, 2006, 2007, 2010 and 2011. Only summer periods are shown from where TSL was estimated by satellite at least three times, including at the beginning of the accumulation season. The ELA is well estimated by TSL observations, except in 2011 when no observations occur within 15 days of the end of the melt season. (b) Satellite-derived MG TSL elevations throughout four summer periods: 1999, 2006, 2008 and 2012. Only summer periods are shown from where TSL was estimated by satellite at least three times, including the beginning of the accumulation season.

To reconstruct the LCG and MG B _a, linear regressions between the estimated AAR and observed B _a were used (for the entire MG the B _a is considered to be accurate within ∼15% (Knudsen and Hasholt, 2004; Reference MernildMernild and others, 2011a)). A linear regression between the AAR and B _a gives the relation

(1)

where s is the slope and AAR₀ is the AAR when B _a = 0. Zero values of AAR were excluded from the regression (only MG experienced years where AAR = 0), since AAR and B _a are not linearly related when net ablation occurs all over the glacier surface (Reference MernildMernild and others, 2011a). Additional information about the LCG B _a program and methods is provided by Reference Marcus, Chambers, Miller and LangMarcus and others (1995), Reference Sapiano, Harrison and EchelmeyerSapiano and others (1998) and Reference Miller and PeltoMiller and Pelto (1999), and about the MG B _a program by Reference Knudsen and HasholtKnudsen and Hasholt (2008) and Reference MernildMernild and others (2011a).

Snow ablation rates and ELA reconstruction

The TSL for LCG was observed for 39 dates during the period that defines 25 time periods during which satellite observations are at least 15 days apart (Table 2). For MG the numbers were 25 dates within 11 time periods (Table 3). For LCG the observed positive TSL migration rates varied from 2.9 ± 0.9 m d⁻¹ (for the 2004 ablation season) to 3.9 ± 0.0 m d⁻¹ (2005) (having a migration rate of up to 5.2 m d⁻¹ between subsequent satellite observations), with a mean for all ablation periods of 3.8 ± 0.6 m d⁻¹ (positive rates indicate when TSL is moving towards higher elevations, and, here and below, ± is stated as plus or minus one standard deviation). At the beginning of the accumulation season, negative LCG TSL migration rates occurred within the range −1.6 to −15.4 m d⁻¹ (Table 2), indicating a lowering in the TSL elevation between September and October. The mean TSL migration rate on LCG of 3.8 m d⁻¹ compares well with the mean migration rate of 3.7 m d⁻¹ on nearby Taku Glacier (Reference PeltoPelto, 2011), a temperate maritime valley glacier located in the Juneau Icefield (671 km²; 58.4° N, 134.1° W), ∼20 km to the northeast of LCG. The larger area of Taku Glacier allows the use of high temporal resolution Moderate Resolution Imaging Spectroradiometer (MODIS) imagery for accurate TSL identification. This provides additional dates closer to the end of the ablation season, and allows application of the TSL migration rate for well-constrained estimates of the snow ablation rates and the location (elevation) and date of the annual ELA.

For MG the observed positive TSL migration rates varied from 5.6 ± 0.1 m d⁻¹ (for the 2000 ablation season) to 14.9 ± 15.1 m d⁻¹ (2012) (having a migration rate up of to 37.3 m d⁻¹ between subsequent satellite observations), with a mean for all ablation periods of 9.4 ± 9.1 m d⁻¹ (Table 3). At the beginning of the MG accumulation season, from the end of August to September/October, TSL migration rates ranged from −0.5 to −18.4 m d⁻¹, illustrating the lowering rate of the TSL (Table 3). The TSL migration rate was used to determine snow ablation rates using both methods: the TSL–mass-balance-gradient method and the snow-pit–satellite method (see Section 3). For LCG, based on the TSL–mass-balance method, the snow ablation rates varied from 0.023 to 0.039 m w.e. d⁻¹, averaging 0.028 ± 0.004 m w.e. d⁻¹, whereas snow ablation rates based on the snow-pit–satellite method varied from 0.025 to 0.038 m w.e. d⁻¹, averaging 0.031 ± 0.004 m w.e. d⁻¹ (Fig. 3a; Table 2). The JIRP ablation measurements for LCG during the 2004–10 ablation seasons, over a total period of 162 days, yield a mean snow ablation rate of 0.031 m w.e.d⁻¹, which is in accordance with calculations: the estimated snow ablation rates for LCG were significantly identical (97.5% quartile; based on the null hypothesis). The similarity of the TSL and field snow ablation rates supports the concept that remote-sensing TSL observations (which can be extended over longer time periods and are not simple point measurements), together with field snow-pit observations, offer a useful approach for estimating annual ablation rates, which are important in assessing changes in glacier mass balance in the Juneau Icefield region.

For MG the snow ablation rates showed more variability than for LCG, with rates in the range 0.037–0.072 m w.e. d⁻¹, averaging 0.051 ± 0.018 m w.e. d⁻¹ (based on the TSL–mass-balance-gradient method), and 0.028–0.073 m w.e. d⁻¹, averaging 0.047 ± 0.019 m w.e. d⁻¹ (based on the snow-pit– satellite method) (Fig. 3b; Table 3). However, the estimated snow ablation rates for MG were significantly identical (97.5% quartile; based on the null hypothesis). At MG no direct field snow ablation measurements have been conducted to validate the estimated snow ablation values, but in future mass-balance model simulations the calculated snow ablation rates have the potential to be compared against simulated ablation rates. Reference Mernild, Liston, Hasholt and KnudsenMernild and others (2006) presented simulated daily snow and ice melt rates using the modeling software package SnowModel (Reference Liston and ElderListon and Elder, 2006; Reference Mernild, Liston, Steffen and ChylekMernild and others, 2010; Reference Liston and MernildListon and Mernild, 2012) for the period 1999–2004, and these simulated rates (0.03–0.04 m w.e. d⁻¹) were on average slightly lower than the estimated snow ablation rates presented in this study (Fig. 3b). A reason for this could be that the TSL–mass-balance-gradient and snow-pit–satellite methods concern the ablation rate of the snowpack (melt, evaporation and sublimation), whereas SnowModel simulations only include surface melt rates from the snowpack and the bare glacier ice (SnowModel simulations forced by mean daily climate data).

The satellite-derived TSL and dates for LCG and MG provide a dataset for estimating the annual ELA and the date of the end of the ablation season. In Figure 4a and b, seasonal variations of TSL are shown for LCG (2004, 2006, 2007, 2010 and 2011) and MG (1999, 2006, 2008 and 2012). A second-order polynomial regression between the day of year (DOY) and TSL on an annual scale indicates that ELA for LCG varies between 1020 m a.s.l. in 2007 and 1170 ma.s.l. in 2011, and the ablation seasons ended between DOY 234 (22 August) and DOY 260 (17 September) for the years where satellite-derived TSL observations through and beyond the ablation season were available. These estimated ELA conditions for LCG are in accordance with annual fieldwork observations (Fig. 2a), although ELA is overestimated on average by ∼50 m a.s.l. compared to direct observations; much of this overestimation occurred from 2011, where no TSL observations were available within 15 days of the end of the melt season. For MG, the estimated annual ELA was located between 460 m a.s.l. (2008) and 780 m a.s.l. (2012), and the ablation season ended between DOY 227 (15 August) and DOY 230 (18 August). The estimated ELA was significantly identical to MG annual field observations (97.5% quartile; based on the null hypothesis), although the estimated ELA was on average underestimated by ∼90 m a.s.l. (Fig. 2b), giving reason to believe that the method presented here is useful for ELA estimations at both MG and LCG. For the years 1999, 2006, 2008 and 2012 the MG ablation season ended within 3–4 days in mid-August. At MG the B _a observations are conducted in early/midAugust, which seems to be a good time for capturing the majority of the ablation season (at least for the four years 1999, 2006, 2008 and 2012 as illustrated in Fig. 4b); however, surface melt occurred considerably later in particular years (e.g. until late October in 2010 and late September in 2012).

AAR and B _a reconstruction

AAR varies greatly from one year to another (Table 4); however, for a period long enough to filter out extremes but shorter than the timescale of adjustment to glacier equilibrium, it gives a measure of the health of the glacier (Reference CogleyCogley and others, 2011). For LCG, observed AAR varied between 0.07 (1998) and 0.82 (2000) for the period 1997/98 to 2011/12, averaging 0.57 ± 0.24, while AAR for MG varied between 0.75 (2003) and 0.00 (e.g. 2012) for the period 1998/99 to 2011/12, averaging 0.15 ± 0.22 (Table 4). MG experienced AAR = 0 six times within the last 14 years, including the three most recent years in this study (2010, 2011 and 2012). According to Reference Dyurgerov, Meier and BahrDyurgerov and others (2009), glaciers and ice caps in equilibrium with the local climate typically have an AAR of 0.5–0.6, with a global average of 0.579 ± 0.009. Reference PeltoPelto (2010) identified that glaciers having a frequent AAR = 0 lack a persistent accumulation zone and cannot survive.

Table 4. Observed and TSL satellite-derived AAR for LCG (1998–2011) and MG (1999–2012)

In Table 4, annual TSL satellite-derived AAR values are listed for LCG (2004, 2006, 2007, 2010 and 2011) and MG (1999, 2006, 2008 and 2012) and compared against field observations. For LCG the TSL satellite-derived method on average underestimated AAR by 0.16 (16%) compare to observations, and for MG AAR overestimated on average by 0.14 (14%), but since the respective error bars overlap, there is no significant difference. In Figure 5a and b, TSL satellite-derived AAR is plotted against B _a for LCG and MG. The TSL satellite-derived AAR B _a trend lines (red lines) follow observed values and trend lines (black lines). If additional observations are added to the trend lines, B _a can be substituted by the satellite observations, once sufficient data exist to better constrain late-season TSL behavior and hence annual AAR determination. Based on the LCG satellite-derived AAR and AAR₀ conditions (Fig. 5a), expected changes in the LCG area and volume can be derived from

Fig. 5. (a) LCG observed AAR and B _a trend line (dashed black line) from 1997/98 to 2010/11, and TSL satellite-derived AAR and observed B _a trend line (dashed red line) based on data from 2003/04, 2005/06, 2006/07, 2009/10 and 2010/11. Also illustrated are standard errors for each dataset. (b) MG observed AAR and B _a trend line (dashed black line) from 1995/96 to 2011/12 (zero values of AAR are excluded from the regression, as AAR and B _a are only linearly related when ELA is located within the elevation range of the glacier), and TSL satellite-derived AAR and observed B _a trend line (dashed red line) based on data from 1998/99, 2005/06, 2007/08 and 2011/12. Also illustrated are standard errors for each dataset.

(2)

the ratio of the current AAR to its equilibrium value. The fractional changes in area (p _s) and in volume (p _v) are calculated from

(3)

(4)

where γ = 1.36 for valley glaciers, derived empirically and from theory (Reference Bahr, Dyurgerov and MeierBahr and others, 2009). Based on the LCG trend between the TSL satellite-derived AAR and B _a, LCG has an estimated AAR₀ of 0.57 that is comparable to the observed AAR₀ value of 0.67 (Fig. 5a). Reference Dyurgerov, Meier and BahrDyurgerov and others (2009) computed AAR₀ for 86 glaciers and ice caps by using linear regression between AAR and B _a, showing an average value of 0.58 ± 0.01. The resulting LCG AAR (0.42; Table 4) and AAR₀ values (0.57; Fig. 5a) (based on the TSL satellite-derived AAR relationship to B _a) indicate that LCG will lose 26 ± 3% of its present area and 34 ± 3% of its volume typically over several decades or longer if current climate conditions in the region of LCG persist. Based on observed AAR (0.57; Table 4) and AAR₀ (0.67; Fig. 5a), LCG will respectively lose 15 ± 1% and 20 ± 2%. Similar area and volume fraction calculations were conducted for MG, indicating that MG based on the TSL satellite-derived AAR (0.33; Table 4) and AAR₀ (0.66; Fig. 5b) will lose about 50 ± 6% of its present area and 61 ± 5% of its volume if current climate conditions in southeast Greenland persist. MG is significantly out of balance with climate, and far below the global AAR mean, and will likely lose a significant amount of its current area and volume even in the absence of further climate changes. Based on α _r calculations from observations in Reference MernildMernild and others (2011a) (AAR = 0.16 (Table 4) and AAR₀ = 0.61 (Fig. 5b)), MG will lose 74 ± 8% of its current area and 84 ± 7% of its volume over several decades or longer if current climate conditions persist. For both glaciers the satellite-estimated fractional areal and volume losses seem to point out the extent to which the glaciers are out of balance with present-day climate observations.

An expansion of the study by adding satellite-derived annual glacier conditions, i.e. ELA, AAR and B _a, is desirable to better quantify the presented relationships, to increase accuracy and further validate the findings at LCG and MG. Also, so-called ‘transient’ area-averaged mass balances can be computed and related to concurrent transient ELA and AAR values; this method assumes that the relationship between transient values of mass balance and ELA and AAR in the course of one season is identical to the relationship between B _a and ELA and AAR at the end of the mass-balance year over many years (Reference Hock, Kootstraa and ReijmerHock and others, 2007).