1. Introduction
The Greenland ice sheet (GrIS) has rapidly lost mass during the past decade (e.g. Reference Jacob, Wahr, Pfeffer and SwensonJacob and others, 2012), with approximately one-third to one-half of this deficit attributed to accelerated discharge from marine-terminating outlet glaciers (e.g. Reference Van den BroekeVan den Broeke and others, 2009; Reference ShepherdShepherd and others, 2012; Reference Enderlin, Howat, Jeong, Noh, Van Angelen and Van den BroekeEnderlin and others, 2014). Although some studies suggest that the relative importance of dynamic losses from the GrIS may slow at centennial timescales (Reference GoelzerGoelzer and others, 2013), recent evidence indicates that they are likely to contribute substantially to 21st-century sea-level rise (Reference NickNick and others, 2013; Reference StockerStocker and others, 2013). Consequently, identifying the dominant controls on outlet glacier dynamics is critical for accurate prediction of 21st-century sea-level rise, yet substantial uncertainty remains (Reference StockerStocker and others, 2013). Greenland outlet glaciers can retreat rapidly in response to external forcing (e.g. Reference Howat, Joughin, Fahnestock, Smith and ScambosHowat and others, 2008; Reference JoughinJoughin and others, 2008a), but this can be strongly modulated by local topographic controls (e.g. basal topography and along-flow fjord width variations) (Reference JamiesonJamieson and others, 2012; Reference Enderlin, Howat and VieliEnderlin and others, 2013; Reference Carr, Stokes and VieliCarr and others, 2014). These factors can produce substantial variations in local retreat rates and highlight the pitfalls of extrapolating retreat patterns into the future and between glaciers (Reference McFadden, Howat, Joughin, Smith and AhnMcFadden and others, 2011; Reference Moon, Joughin, Smith and HowatMoon and others, 2012; Reference Carr, Stokes and VieliCarr and others, 2014). Although the potential impact of basal topography on marine-terminating glacier dynamics has long been recognized (e.g. Reference WeertmanWeertman, 1974; Reference Meier and PostMeier and Post, 1987), its influence on the contemporary and near-future behaviour of major Greenland outlets has not been systematically assessed due to the limited availability of high-resolution subglacial topographic data.
Here we investigate the multi-decadal behaviour of Humboldt Glacier (HG; Fig. 1), northern Greenland, in relation to basal topography and atmospheric and surface oceanic forcing (the term ‘surface oceanic’ includes sea ice and sea surface temperatures). HG is the widest marine-terminating outlet in Greenland (∼90 km; Reference Weidick, Williams and FerrignoWeidick, 1995) and drains ∼3% of the ice sheet (Reference Rignot and KanagaratnamRignot and Kanagaratnam, 2006). Previous studies reported negative mass balance, grounding-line retreat and thinning during the 1990s (Reference Joughin, Fahnestock, Kwok, Gogineni and AllenJoughin and others, 1999; Reference AbdalatiAbdalati and others, 2001; Reference Rignot, Gogineni, Joughin and KrabillRignot and others, 2001), but HG’s recent dynamics and the factors driving its behaviour have yet to be investigated. Moreover, ice thickness profiles indicate a deep bed beneath the northern section of the terminus (Reference Joughin, Fahnestock, Kwok, Gogineni and AllenJoughin and others, 1999), and more recent work suggests that the bed is below sea level within ∼140 km of the calving front (Reference Morlighem, Rignot, Mouginot, Seroussi and LarourMorlighem and others, 2014). Here we use new high-resolution (500 m) bed data coupled with numerical modelling to assess the impact of basal topography on HG’s response to external forcing.
Retreat of Humboldt Glacier (HG) for selected years between 1975 and 2012 for (a) the entire terminus and (b) the northern section. Frontal positions are colour-coded by year and are as close as possible to 31 July. Image source: Landsat, USGS GLOVIS (28 June 2000). Inset 1: Location of HG and Qaanaaq/Qaanaaq Mittarfik meteorological stations on the Greenland ice sheet. Inset 2: Location of HG catchment within surrounding topography and glacial features. Image source: MODIS Aqua, USGS GLOVIS (5 August 2012).
2. Methods
2.1. Calving front positions
Frontal positions were digitized manually from synthetic aperture radar (SAR) Image Mode Precision imagery, provided by the European Space Agency (1992–2012), and from Landsat data from the US Geological Survey’s (USGS) Global Visualization Viewer (GLOVIS) (1975–2012). Following previous studies (e.g. Reference Moon and JoughinMoon and Joughin, 2008), frontal positions were digitized from sequential images, within a reference box of fixed width and upstream extent (Fig. 2), and mean retreat was calculated by dividing the change in reference box area by its width. Positions were calculated relative to 6 August 1975 (the earliest image date). Errors were evaluated by repeatedly digitizing ten sections of rock coastline proximal to HG’s terminus from a subset of 15 SAR and Landsat images, and resulted in a mean error of 42.4 m.
Basal topography and surface elevation profiles for HG. (a) Basal topography of HG, with frontal positions (colour-coded as in Fig. 1), transects (thick black line), catchment areas for model transects (black-and-white dashed line), frontal position reference boxes (thin dark grey line) and potential basal pinning point (circled in black). (b) Mean ice velocity field for HG (colour scale) from MEaSUREs data (Reference Joughin, Smith, Howat, Scambos and MoonJoughin and others, 2010b), overlain on basal topography (grey contours), and potential basal pinning point (circled in black). (c) Change in ice velocity for HG between winter 2000/01 and 2008/09 (red-yellow: acceleration; blue: deceleration) from MEaSUREs data (Reference Joughin, Smith, Howat, Scambos and MoonJoughin and others, 2010b), overlain on basal topography (grey contours), and potential basal pinning point (circled in black). (d, e) Surface elevation, basal topography and flotation elevation for transects along (d) transect 1 and (e) transect 2.
2.2. Atmospheric and oceanic data
2.2.1. Air temperature data
Surface air temperature (SAT) data were obtained from Qaanaaq (77°28′N, 69°13′W) and Qaanaaq Mittarfik (77°29′N, 69°23′W) meteorological stations (Fig. 1, inset 1) and were provided by the Danish Meteorological Institute at a 3 hourly temporal resolution (Reference Carstensen and JørgensenCarstensen and Jørgensen, 2011). Data were filtered to account for missing values and were only used to calculate monthly/annual averages if the criteria were met (Reference CappelenCappelen, 2011), namely, (1) no more than two consecutive records were missing in a day; (2) no more than three records in total were missing in a day; (3) daily averages were available for 22 or more days per month; and (4) monthly averages were available for all months of the year. In order to extend the temporal coverage of the data, records were used from Qaanaaq between January 1996 and August 2001 and from Qaanaaq Mittarfik between September 2001 and December 2010. The two stations are located 1.9 km apart and at the same elevation (16 m.a.s.l.). The variation between the two stations was assessed by comparing mean monthly values for the period of data overlap (August 2001–October 2004), and the average difference was 0.28°C. These data were used to calculate mean summer (June–August (JJA)) and mean annual air temperatures and the number of positive degree-days per year.
As meteorological data were not available for the entire study period and Qaanaaq/Qaanaaq Mittarfik are located ∼190 km from HG (Fig. 1, inset 1), reanalysis data were also used to assess air temperature changes. The data were obtained from US National Centers for Environmental Prediction (NCEP)/US National Center for Atmospheric Research (NCAR) Reanalysis 1, which has a spatial resolution of 2.5° (∼250 km × 280 km at 79°N) (Reference KalnayKalnay and others, 1996), and the European Centre for Medium-Range Weather Forecasts’ ERA-Interim reanalysis, which has a spatial resolution of 0.75° (∼75 km × 80 km at 79° N) (Reference DeeDee and others, 2011). NCEP/NCAR data are used for the period 1975–2010, to coincide with the frontal position record, and ERA-Interim data from the start of the dataset in 1979 until 2010. In both cases, air temperature values from the 700 mbar geopotential height were used, as their correlation with meteorological station data is better than values from 2 m height (personal communication from A. Gardner, 2013) and the influence of sea surface temperature (SST) on air temperature values is reduced (Reference Moholdt, Wouters and GardnerMoholdt and others, 2012).
2.2.2. Sea-ice data
Sea-ice data were extracted from charts provided by the the US National Ice Center (NIC) (http://www.natice.noaa.gov/), which are compiled from a range of remotely sensed and directly measured data sources and have a spatial resolution of up to 50 m. Sea-ice values were sampled from a polygon extending the full width of the glacier terminus, as defined in the land mask of the sea-ice product, and 50 m perpendicular to it. The estimated accuracy of sea-ice concentrations is ±10% (Reference Partington, Flynn, Lamb, Bertoia and DedrickPartington and others, 2002), and the use of multiple sources and manual interpretation is thought to improve upon the accuracy provided by a single data source.
2.2.3. Sea surface temperature data
SST data were obtained from two sources: (1) monthly SST climatology products acquired by the Moderate Resolution Imaging Spectroradiometer (MODIS) at spatial resolution of 5 km, which were provided by the NASA Ocean Color Project (http://oceancolor.gsfc.nasa.gov/), and (2) the Reynolds SST analysis dataset Version 2 at a spatial resolution of 0.25° (∼25 km × 28 km at 79°N), provided by the US National Oceanic and Atmospheric Administration (NOAA) (Reference Reynolds, Smith, Liu, Chelton, Casey and SchlaxReynolds and others, 2007). SST values were sampled from all grid squares within 25 km of the terminus of HG from the MODIS data. Due to the comparatively coarse resolution of the Reynolds data, SSTs were sampled from the grid squares closest to HG’s terminus.
2.2.4. Statistical analysis
For all atmospheric and oceanic forcing data, the statistical significance of trends was assessed using the t-statistic. The t-statistic is used to evaluate whether the coefficient associated with a given independent variable is significantly different from zero, given the variation in the data, and is calculated by dividing the coefficient by the standard error. In this case, the t-statistic is used to evaluate whether the trend in a given forcing factor (e.g. mean annual air temperatures) over time is significantly different from zero, taking into account the variations in individual values of that forcing factor. The t-statistic has an associated p-value, which tests the probability of obtaining a value of the t-statistic that is at least as extreme as the one observed, if the null hypothesis is true. Here we use a significance interval of 0.05 (i.e. a confidence interval of 95%), so that p-values ±0.05 show that the coefficient is significantly different from zero.
Trends in mean annual and mean summer (JJA) air temperatures were calculated from both reanalysis datasets for the following time periods: the entire study period (1975–2010 for NCEP/NCAR data and 1979–2010 for ERA data); prior to rapid retreat (1975/79–98); and after the onset of rapid retreat (1999–2010). In addition, a paired t-test was used to test for statistical difference between mean annual air temperature for the period 1999–2008 (the decade following onset of rapid retreat) and the preceding decades between 1949 and 1998 (all preceding decades for which data were available) (Table 1). NCEP reanalysis data were used, as they provide the longest record. A one-tailed test and a significance interval of 0.05 were used, so that test results with a p-value of 0.05 or less show that a given decade was significantly cooler than the period 1999–2008.
Decadal mean air temperatures and paired t-test results, comparing air temperatures between 1999 and 2008 with preceding decades. Test results with a p-value of ±0.05 show that the decade was significantly cooler than 1999–2008
2.3. Basal topography, surface elevation and ice velocity data
Airborne radar measurements of ice thickness were obtained and processed as part of the Greenland Outlet Glacier Geophysics (GrOGG) project (Reference PalmerPalmer and others, 2013). Data were acquired in May 2012 using the High-Capability Radar Sounder (HiCaRS). Survey lines were flown parallel to the ice-flow direction, at a constant height of 800m above the ice surface, and perpendicular to the ice-flow direction, at a constant altitude for each perpendicular survey line. Ice thickness measurements from previous airborne surveys by the Center for Remote Sensing of Ice Sheets (CReSIS), University of Kansas, USA (Reference GogineniGogineni and others, 2001), and NASA’s Operation IceBridge (Reference Koenig, Martin, Studinger and SonntagKoenig and others, 2010) were incorporated to improve data coverage. Bed elevations were calculated by interpolating all ice thickness data onto a 500 m grid, removing areas >5 km from any survey line and then subtracting them from the combined laser-altimeter/ Greenland Ice Mapping Project (GIMP) ice surface digital elevation model (DEM) (Howat and others, 2012), also resampled to 500m (Reference PalmerPalmer and others, 2013). The basal topographic data are the highest spatial resolution data available for the study region, with a grid spacing of 500 m, compared to the most recent Greenland bed elevation dataset which has a resolution of 1–2.5 km (Reference BamberBamber and others, 2013). These data are previously unpublished and are yet to be incorporated into ice-sheet-wide compilations of basal data (e.g. Reference BamberBamber and others, 2013).
Flight-line crossover points were used to evaluate the vertical errors in ice thickness measurements and basal elevations within the study area. The root-mean-square (rms) for the ice thickness measurements was 8.4m (from 27 crossovers) for the HiCaRS radar data and 5.8 m (from 13 crossovers) for the Operation IceBridge/CReSIS data. As detailed above, the bedrock topography was derived by subtracting ice thickness measurements from the GIMP DEM (Howat and others, 2012), so errors from both datasets need to be accounted for when calculating the errors associated with basal elevations. The mean ICESat (Ice, Cloud and land Elevation Satellite) validation error for the GIMP DEM tiles covering HG is 17.7m (Howat and others, 2012). This value is combined with ice thickness errors in quadrature, as the two error sources are independent, to give a total rms for the basal elevations of 21.9 m for the HiCaRS data and 20.6 m for the Operation IceBridge/CReSIS data. This compares favourably with the error values of Reference BamberBamber and others (2013), which range between ±10m and ±300 m.
Surface elevation data were obtained from the GIMP DEM, which is constructed from ASTER (Advanced Space-borne Thermal Emission and Reflection Radiometer) and SPOT-5 (Satellite Pour l’Observation de la Terre) DEMs for the study area (Howat and others, 2012). The data are registered to ICESat elevations from 2003 to 2009, so the DEM has a nominal date of 2007. The rms error (RMSE) across the ice sheet, relative to ICESat, is ±10 m. This ranges from ∼±1 m over most ice surfaces to ±30m in areas of high relief (Howat and others, 2012). As detailed above, the mean RMSE in the study region is 17.7m (Howat and others, 2012). Given the nominal date of the DEM, surface profiles were taken from the closest frontal position to midsummer in 2007 (24 July 2007), in order to avoid including sea ice in the elevation profile. Surface elevations were sampled at 500m intervals from the 2007 terminus, along transects following the approximate ice flow direction.
Surface and basal elevations were sampled at 500m intervals from the terminus, along two flow-parallel transects (Fig. 2). The transects represent the two end-members of the basal topography at HG: transect 1 follows the deepest route through the basal trough in the northern section, while transect 2 overlies the shallowest basal topography, beneath the southern section of the terminus (Fig. 2). Basal and surface elevations were used to calculate flotation elevations along each transect, assuming ice density is 910 kg m−3 and ocean density is 1028 kg m−3. Ice surface velocities were obtained from the winter ice-sheet-wide velocity maps for the GrIS from the MEaSUREs (Making Earth Science Data Records for Use in Research Environments) programme (Reference Joughin, Smith, Howat, Scambos and MoonJoughin and others, 2010a). A mean velocity field for HG was calculated using all available velocity maps (winter 2000/01, 2005/06, 2006/07, 2007/08 and 2008/09). Velocity changes were calculated between winter 2000/01 and 2008/09. Only velocity changes greater than ±15 m a−1 are included, based on error estimates from Reference Joughin, Smith, Howat, Scambos and MoonJoughin and others (2010a) and from visual inspection of areas of rock coastline proximal to HG’s terminus, which should demonstrate no discernible acceleration/deceleration yet show changes of ±15 m a−1 in the dataset.
3. Numerical Modelling
3.1. Model description
For the numerical modelling experiments we use a one-dimensional flowline model which is described in detail in Reference Nick, Van der Veen, Vieli and BennNick and others (2010) and has previously been used to investigate the behaviour of large Greenland outlet glaciers and their response to external forcing (Reference Vieli and NickVieli and Nick, 2011; Reference NickNick and others, 2012, Reference Nick2013). It calculates the evolution of frontal positions, ice surface, ice flow and stress field for the two along-flow transects detailed above (Fig. 2). In the model the driving stress (τ
d) is balanced by resistive stresses from the base (τ
b) and lateral margins (τ
lat) and by the along-flow stress gradients, 
, in the ice flow direction (x):
We assume a nonlinear sliding relation (Reference WeertmanWeertman, 1957) that is a function of the effective pressure at the bed (N). The stress balance (Eqn (1)) gives the following expression for depth- and width-averaged ice flow (u), for an ice thickness of (H) and half-width (W):
where ρ i is ice density (910 kg m−3), g is acceleration due to gravity, S is the ice surface, β is the basal sliding coefficient (Reference WeertmanWeertman, 1957), n and m are the exponents for ice flow and sliding relations, respectively, and are both set to 3, and f lat is the buttressing factor. A is the flow rate factor (Reference GlenGlen, 1955) that relates to ice rheology, and is set to a value of 1.16 × 10−18 Pa−3 a−1, applicable to an ice temperature of −30°C. Equation (2) is solved for the width-averaged ice flow (u) by iterating for effective viscosity (ν):
Variations in surface accumulation (a) and width (W) with distance along flow are explicitly accounted for in calculations of surface elevation change along the flowline:
The model allows the glacier terminus to move freely using a physically based calving model (Reference Nick, Van der Veen, Vieli and BennNick and others, 2010), and the approach to simulating grounding line motion is consistent with boundary layer theory (Reference SchoofSchoof, 2007), overcoming previous issues relating to model numerics (Reference Vieli and PayneVieli and Payne, 2005). The horizontal grid resolution adjusts with each time step, so that the grounding line position can be accurately tracked over time (Reference Vieli and PayneVieli and Payne, 2005). The model includes a dynamic calving model, which is the best currently available and is based on the depth of both surface and basal crevasses: calving occurs when surface and basal crevasses penetrate the full ice thickness (Reference Nick, Van der Veen, Vieli and BennNick and others, 2010). Using this calving criterion, the terminus is not necessarily at flotation, but can instead be above flotation. Where the glacier is at flotation, which is usually the case very close to the calving front, the velocity boundary condition at the calving front is calculated from the longitudinal stress that balances the difference in hydrostatic pressure between the ice front and the ocean (Reference Vieli and PayneVieli and Payne, 2005), and is given by
where horizontal velocity gradient with along-flow distance x is evaluated at the terminus. The parameter ρ w is the density of ocean water (1028 kg m−3)and H f is ice thickness at the terminus. The parameter f s is used to scale the rate factor (A) in the perturbation experiments described in Section 2.4.2. It is used to simulate sea-ice-induced changes in longitudinal strain rates and is set to 1 in the reference state, where no perturbation is applied (Reference NickNick and others, 2013).
3.2. Perturbation experiments
The modelling does not aim to reproduce or predict the evolution of HG, but instead explores differences in the sensitivity of the two transects to external forcing (Fig. 2). We built the initial reference states from remotely sensed data (Table 2), and values were calculated separately for each transect (e.g. the width was calculated across the area draining into the northern and southern sections respectively, using the drainage basin extent for HG from Reference Joughin, Fahnestock, Kwok, Gogineni and AllenJoughin and others (1999) to delineate the outer margins and a combination of the surface DEM and Landsat imagery to identify the divide between the two basins). The model was run with these input parameters for 2000 years and had an initial grid resolution of 500 m. The basal sliding coefficient was then adjusted so that modelled retreat rates approximate those observed at each transect between 1976 and 1999 (i.e. prior to rapid retreat).
Input parameters and data sources for numerical modelling experiments
To these initial states we applied perturbations associated with changes in atmospheric and sea-ice forcing, which our observations suggest are primary controls (Section 4; Fig. 3). Air temperature increase is approximated by increasing the water-level parameter within crevasses and hence increasing the crevasse depth (Reference Nick, Van der Veen, Vieli and BennNick and others, 2010, Reference Nick2013). Following Reference NickNick and others (2013), changes in the potential buttressing effect of sea ice, which could result from weakening and/or thinning of the sea ice, are applied by altering the rate factor (A, here 1.16 × 10−18 Pa−3 a−1) using the parameter f s (Eqn (5)). Initially, f s is set to a value of 1, when no perturbations are applied. Higher values of f s increase the longitudinal strain rate at the terminus and thus simulate a reduction in sea-ice buttressing (Reference NickNick and others, 2013). SSTs showed no significant change during the study period and are therefore unlikely to be a primary control on retreat. They were therefore not included in the numerical modelling experiments.
Relative frontal position and forcing factors at HG. (a) Frontal position of HG, relative to 6 August 1975 for the full terminus (black), northern (red) and southern section (blue). (b) Mean annual and mean summer (JJA) air temperatures from reanalysis data. (c) Mean annual and mean summer (JJA) air temperatures from reanalysis data. (d) Mean annual and mean summer (JJA) air temperatures and positive degree-days (PDD) at Qaanaaq/Qaanaaq Mittarfik. (e) Mean seasonal sea-ice concentrations for December–February (DJF), March–May (MAM), June–August (JJA) and September–November (SON). (f) Mean summer (JJA) sea surface temperatures from MODIS (blue) and Reynolds (green line).
Step changes in crevasse water depth and sea-ice buttressing were applied to each transect after an initial phase of 100 years with no perturbation: crevasse water depth was increased by 5, 10, 15, 20 and 25 m, and reduced sea-ice buttressing was simulated by increasing the factor f s to 1.05, 1.1, 1.2, 1.5 and 2.0, so that f s = 2.0 increases the terminus longitudinal strain rate by a factor of two (Reference NickNick and others, 2013). The model was then run for a further 100 years. Crevasse water depth is difficult to estimate, so we use values up to 10m as potentially realistic estimates (cf. Reference Cook, Zwinger, Rutt, O’Neel and MurrayCook and others, 2012) and then explore more extreme values. Note that we assess the differing sensitivity of the two sections to the same perturbation, rather than predicting the response of a single profile to forcing of a particular magnitude, so the absolute perturbation values selected should have minimal impact on results.
4. Results
Terminus-wide retreat rates averaged 81 m a−1 between 1975 and 2012 (Figs 1 and 3). Retreat showed two distinct phases: rates were relatively low between 1975 and 1999 (mean 37 m a−1), but thereafter increased substantially (mean 162 m a−1) for the period 1999–2012. The behaviour of the northern and southern sections of the terminus differed markedly (Fig. 1). Between 1975 and 2012, mean retreat rates were an order of magnitude greater in the northern section of the terminus (147 m a−1) than in the southern section (15 m a−1). This difference persisted between 1999 and 2012, when retreat rates averaged 289 m a−1 in the northern section, compared to 35 m a−1 in the southern section.
During the period between winter 2000/01 and winter 2008/09, the northern section flowed significantly faster (150–570 m a−1) than the southern section (<150 m a−1) (Fig. 2b), as also found in previous work (Reference Rignot, Gogineni, Joughin and KrabillRignot and others, 2001; Reference Joughin, Smith, Howat, Scambos and MoonJoughin and others, 2010b). The highest flow rates occurred at the terminus, ∼22 km south of the northern margin, and were located ∼7 km south of the area of deepest basal topography and most rapid retreat (Fig. 2b). Velocity changes between winter 2000/01 and 2008/09 varied spatially: acceleration occurred at the southern end of the large embayment formed in the northern section of the terminus and reached a maximum of 37 m a−1 (Fig. 2c). The maximum deceleration was observed close to the northern margin of HG’s terminus, near the tip of the northern arm of the embayment, and reached values of up to 31 m a−1 (Fig. 2c).
The northern section of the terminus is situated in a basal trough extending up to 72 km inland (Fig. 2a). Its depth is generally >300 m, reaching a maximum of 475 m below sea level ∼39 km inland from the terminus. In contrast, the bedrock beneath the southern section is comparatively shallow, reaching a maximum depth of 220 m (Fig. 2a). On the northern transect, the terminus is close to flotation up to 6.5 km inland (Fig. 2d), whereas the southern section only approaches flotation within 500 m of the terminus (Fig. 2e).
Frontal retreat coincided with a statistically significant linear warming trend in mean annual air temperature of 0.05°C a−1 (p < 0.01) in both the NCEP/NCAR and ERA-Interim reanalysis data between 1975 (NCEP/NCAR)/1979 (ERA) and 2010 (Fig. 3b; Table 3). A significant warming trend of between 0.06 (NCEP/NCAR) and 0.07°C a−1 (ERA) (p < 0.01) was also observed in both datasets for summer (JJA) air temperatures (Fig. 3c; Table 4). No statistically significant trends in mean annual or mean summer air temperatures were identified in either dataset prior to the onset of rapid retreat (1975/79–98) (Fig. 3b and c; Tables 3 and 4). In contrast, significant warming was observed in both datasets between 1999 and 2010, coincident with rapid retreat, with mean annual air temperatures warming by between 0.12 (ERA) and 0.17 (NCEP/NCAR)°C a−1 and mean summer temperatures increasing by 0.17 (ERA) to 0.19 (NCEP/NCAR)°C a−1 (Fig. 3b and c; Tables 3 and 4). In addition, mean annual air temperatures were anomalously high during the year before the onset of retreat (1998), at 1.3°C above the 1975–2010 mean (Fig. 3b).
Linear trends in mean annual air temperatures for NCEP/NCAR and ERA-Interim reanalysis data for the periods 1975/79–2010 (entire study period), 1975/79–2010 (prior to rapid retreat) and 1999–2010 (rapid retreat). For each trend, the t-statistic, p-value and R 2 value are given. t is calculated by (coefficient/standard error) and tests whether the coefficient is significantly different from zero, given the variability in the data. p gives two-tailed p-values, which test the probability of obtaining a value as great as, or greater than, the observed value if the null hypothesis is true (i.e. the coefficient value is zero). A p-value of 0.05 was used to identify statistically significant results, which are shown in bold. R 2 provides a measure of how well the linear trend fits the data points
Same as Table 3, but for mean summer (JJA) air temperatures
Paired t-tests showed that mean air temperature for the 10 years following the onset of retreat (1999–2008) was significantly higher than any of the preceding 10 year intervals from 1949 to 1999, at a confidence interval of ±0.05 (Table 1). Although the meteorological record is incomplete (spanning 1996–2010), it showed a similar trend of 0.1°C a−1 (p = 0.02; Fig. 3d). Atmospheric warming was particularly marked in summer (JJA), increasing by 0.2°C a−1 between 1996 and 2010 (p < 0.01) (Fig. 3d). The number of positive degree-days (PDDs) almost trebled, from 208 in 1996 to 597 in 2010 (Fig. 3d). Interannual data demonstrate that the onset of seasonal retreat closely followed air temperature rising above freezing (Fig. 4a and c), with retreat generally beginning in June, and air temperature rising above freezing in late May.
Glacier frontal position relative to 1975 (black crosses) and seasonal climatic/oceanic controls for the northern (a, b) and southern sections (c, d). (a, c) Mean monthly air temperatures for Qaanaaq/Qaanaaq Mittarfik: red: >0°C; blue: <0°C. (b, d) Mean monthly sea-ice concentrations: blue: fast ice (100%); red: all other values.
Frontal retreat was concurrent with reduced summer (JJA) and autumn (September–November (SON)) sea-ice concentrations (Fig. 3a and e). Between 1995 and 2010, summer sea-ice concentrations reduced by 1.5% a−1 (p = 0.00), and autumn values declined by 0.9% a−1 (p = 0.05). The largest seasonal retreats occurred during 2007 and 2009, when summer sea-ice concentrations were low (Fig. 4b and d). SSTs from both the MODIS and Reynolds datasets showed no statistically significant trends during the study period (Fig. 3f).
Modelled retreat rates for step reductions in sea-ice buttressing were more than an order of magnitude greater on the northern than on the southern transect (Fig. 5). For example, retreat was 13 times greater on the northern transect when f s was increased to 1.1. Increasing crevasse water depth by 5 and 10 m produced retreat rates 15 and 16 times greater on the northern section, respectively. For larger increases in crevasse water depth, however, the difference in retreat rates between the two sections was less pronounced (Fig. 5).
Numerical modelling experiments showing changes in relative frontal position over time for the northern (transect 1) (a, c) and southern (transect 2) transects (b, d). Step increases in crevasse water depth (a, b) and reductions in sea-ice buttressing (simulated using front-stress perturbations) (c, d) were applied after 100 years.
5. Discussion
5.1. Climatic and oceanic controls on glacier frontal position
The association between frontal position change and atmospheric warming (Figs 3 and 4; Tables 3 and 4) suggests that seasonal and interannual retreat at HG were linked to higher air temperatures. This is consistent with the coincidence of atmospheric warming from the mid-1990s with glacier retreat and dynamic changes elsewhere in Greenland (e.g. Reference Howat, Joughin, Fahnestock, Smith and ScambosHowat and others, 2008; Reference Moon and JoughinMoon and Joughin, 2008; Reference Howat and EddyHowat and Eddy, 2011; Reference Bevan, Luckman and MurrayBevan and others, 2012). Elevated summer air temperatures at HG would enhance meltwater availability and may thus increase the frequency of hydrofracture of terminus crevasses, promoting calving and retreat. Hydrofracture of crevasses has previously been linked to retreat at Jakobshavn Isbræ, via its impact on terminus crevassing (Reference Sohn, Jezek and Van der VeenSohn and others, 1998) and/or weakening of lateral shear margins (Reference Vieli and NickVieli and Nick, 2011). Evidence from satellite imagery supports this mechanism, as numerous water-filled lakes and crevasses are present within ∼25 km of HG’s northern terminus and are particularly prevalent over the lower ∼ 7 km (Fig. 6b), which is near flotation (Fig. 2c). Notably, however, no water is observed in rift-like crevasses that form within a few hundred metres of the ice front, and from which the large, tabular icebergs subsequently calve (Fig. 6c and d). We therefore suggest that hydrofracture may not necessarily cause the icebergs to detach from the terminus, but instead may open and deepen fractures near the front, thus facilitating calving once the ice reaches the terminus. This is consistent with our numerical modelling results, which indicate that the northern section of HG is highly sensitive to raising the water-level parameter within crevasses and hence increasing the crevasse depth. In addition to hydrofracture, meltwater may facilitate retreat by increasing the outflow of subglacial melt at the calving front: buoyant subglacial meltwater plumes are thought to increase submarine melt rates by promoting a compensatory inflow of warm water at depth (e.g. Reference Rignot, Koppes and VelicognaRignot and others, 2010; Reference JenkinsJenkins, 2011; Reference StraneoStraneo and others, 2011) and are visible at the terminus of HG (Fig. 6a). However, the lack of temperature and salinity data from HG’s fjord precludes more detailed assessment of the effect of these plumes.
Subglacial meltwater plume and the pattern of iceberg calving/rift formation in northern section of HG’s terminus. Background Landsat scene acquired 22 July 1999. (a) Turbid meltwater plume exiting the calving front. Background Landsat scene acquired 28 June 2003. (b) Water-filled crevasses and lakes close to the ice margin. Background Landsat scene acquired 22 July 2000. (c, d) Large tabular icebergs calving back to rifts over a 2 day period. Background Landsat scenes acquired 22 July 1999 (c) and 24 July 1999 (d). All Landsat scenes were provided by Global Land Cover Facility.
Retreat at HG also coincided with a statistically significant decline in summer (JJA) sea-ice concentrations (Fig. 3). This may contribute to retreat by extending the duration of seasonally high calving rates (Reference JoughinJoughin and others, 2008a; Reference Amundson, Fahnestock, Truffer, Brown, Lüthi and MotykaAmundson and others, 2010; Reference Howat, Box, Ahn, Herrington and McFaddenHowat and others, 2010). Numerical modelling supports the assertion that sea-ice concentrations are an important control on terminus behaviour at interannual timescales (Figs 3 and 5). The influence of sea ice is partially supported by intra-annual data, because the onset of seasonal retreat sometimes coincides with springtime sea-ice loss (e.g. 2005 and 2009; Fig. 4). However, this relationship is not perennial (e.g. 2007), and seasonal retreat persisted after winter sea-ice formation in 2003 and 2007 (Fig. 4). This concurs with satellite observations of calving while sea ice is present (Fig. 6b and c). Moreover, large seasonal retreats occurred during years of high summer sea-ice concentrations (e.g. 2003), and a return to pre-retreat concentrations in 2002 and 2004 had little impact on retreat rates (Fig. 4). This suggests that sea ice influences frontal position at HG, but that this relationship is complex and may be insufficient to strongly modulate the frontal position during ongoing retreat.
5.2. Topographic controls on terminus response to forcing
Despite being subject to the same external forcing, the northern and southern sections of HG exhibited very different retreat rates (Fig. 1). We attribute this difference to the underlying basal topography (Fig. 2). The bed beneath the northern section is deeper (max. depth 475 m) than the southern section (max. depth 220 m) and is reverse-sloping in places, with the ice approaching flotation along the lowest 6.5 km (Fig. 2). This agrees with observations from the northern section: (1) calving is predominantly via large, tabular icebergs which are linked to floating termini (Reference Amundson, Fahnestock, Truffer, Brown, Lüthi and MotykaAmundson and others, 2010); (2) the terminus frequently calved back to large surface rifts (Fig. 6), which are associated with near-floating ice (Reference JoughinJoughin and others, 2008b); and (3) the surface profile is nearly flat (∼0.4°) near the terminus (Fig. 2d).
The near-flotation and deeper basal topography of the northern section have important implications for HG’s response to forcing. Basal shear stresses are low for areas close to flotation, so the relative contribution of longitudinal stresses to the force balance is expected to be higher than in well-grounded areas (Reference Echelmeyer, Harrison, Larsen and MitchellEchelmeyer and others, 1994). As the northern section is wide (∼45 km), it is unlikely that these stresses are transferred to the lateral margins. Increased longitudinal stresses produce higher longitudinal stretching rates, which together with dynamic thinning promote crevasse formation (Reference Van der VeenVan der Veen, 1998a). Thus HG’s terminus is expected to be more vulnerable to calving via hydrofracture and full-thickness crevassing, as both basal and surface crevasses could form (Reference Van der VeenVan der Veen, 1998a,Reference Van der Veenb). Moreover, variations in longitudinal stresses associated with changes in sea-ice buttressing would have a greater influence on retreat rates and calving.
Where the terminus is close to flotation, longitudinal stresses at the calving front increase linearly with ice thickness, and stretching rates increase with thickness to the power 3 (Reference SchoofSchoof, 2007). Due to its deeper bed and thicker ice (Fig. 2), the northern section therefore undergoes higher stretching rates and longitudinal stresses at the terminus, which will promote crevasse formation, calving and retreat. Furthermore, flotation close to the terminus will, at least temporarily, increase the area exposed to submarine melting. Although less significant, the grounding line of the northern section is located in deeper water and would therefore experience subsurface temperatures that are ∼0.11°C higher, due to the pressure dependence of the melting point of ice (Reference Rignot and SteffenRignot and Steffen, 2008; Reference Rignot, Koppes and VelicognaRignot and others, 2010).
5.3. Velocity changes
Previous retreats of Greenland outlet glaciers with substantial floating sections, such as Jakobshavn Isbræ (e.g. Reference Joughin, Abdalati and FahnestockJoughin and others, 2004, Reference Joughin2008b) and Alison Glacier (Reference Carr, Vieli and StokesCarr and others, 2013), have been accompanied by marked acceleration of inland ice. In contrast, HG exhibits a dichotomous pattern of acceleration within the main embayment and deceleration near to its northern lateral margin (Fig. 2c). We suggest that the observed deceleration close to the northern lateral margin may result from the terminus retreating onto a small area of shallower bedrock, which may act as a basal pinning point and thus cause deceleration by increasing resistive stress (Fig. 2a and c). This idea is supported by the lower retreat rates in this area, in comparison to the northern section as a whole (Fig. 2). The additional resistance to flow here may also account for the location of the maxima in ice velocity and acceleration to the south of the area of highest retreat and deepest basal topography (Fig. 2a–c).
The acceleration in the main embayment of the northern section (37 m a−2) between winter 2000/01 and 2008/09 is smaller than changes recently observed on Hagen Bræ and Academy Glacier, northern Greenland: the two glaciers accelerated by ∼80 m a−2 between winter 2000/01 and 2005/06, following retreats of 5.4 km and 0.9 km, respectively (Reference Joughin, Smith, Howat, Scambos and MoonJoughin and others, 2010b). The relatively muted response of HG’s inland ice to retreat may be due to the potential pinning point at the northern margin (Fig. 2a and c) and/or the large width of the near-floating northern section, which would limit the transfer of resistive stresses to the lateral margins and would thus reduce the impact of retreat on ice velocities in this area. However, the data indicate that changes in frontal position at HG do influence inland ice velocities to some extent (Fig. 2c), in contrast to nearby Petermann Glacier (Reference NickNick and others, 2012) and C.H. Ostenfeld Glacier (Reference Joughin, Smith, Howat, Scambos and MoonJoughin and others, 2010b), where calving of large sections of ice tongue appears to have limited effect on ice velocities. This suggests that continued retreat at HG would cause acceleration of inland ice, but that its impact is less than on certain other glaciers in northern Greenland, potentially due to its pinning point. Notably, however, the basal topography deepens behind the potential pinning point at HG (Fig. 2a), highlighting the possibility of rapid retreat and ice acceleration in this area, once the terminus retreats beyond that point.
5.4. Numerical modelling of glacier sensitivity to forcing
Numerical modelling results confirm the increased sensitivity of the northern section to forcing (Fig. 5). For all magnitudes of sea-ice-induced increases in longitudinal strain rate, modelled retreat rates were at least an order of magnitude greater on the northern section (Fig. 5). This likely reflects the greater influence of changes in longitudinal stress on the northern section. Increasing crevasse depth by raising crevasse water levels (as a proxy for enhanced surface melt) by 5–10m also resulted in an order-of-magnitude difference in retreat rates between the two sections (Fig. 5). This is consistent with observations: satellite imagery shows that numerous water-filled crevasses and lakes are present on the northern and southern sections, within ∼25 km of the calving front, but retreat rates are an order of magnitude higher on the northern section (Figs 1 and 3). In the numerical model, the response of the two sections of the terminus was more similar when crevasse water depth was increased by 15 m or more (Fig. 5). This increased sensitivity of the southern section likely reflects the local bed and ice surface topography, because the ice is thinner 0.5 km inland than at the terminus, so that a greater step change in crevasse water depth may cause more rapid retreat.
Our results demonstrate that the modelled glacier is strongly sensitive to crevasse water level, via its influence on crevasse depth. This is consistent with previous studies that applied a similar model at Columbia Glacier, Alaska, USA (Reference Cook, Zwinger, Rutt, O’Neel and MurrayCook and others, 2012), and Jakobshavn Isbræ (Reference Vieli and NickVieli and Nick, 2011). It highlights the need to further assess the relationship between crevasse depth, crevasse water level and calving: model outputs are highly dependent on the input water level, due to its control on crevasse depth, but this parameter is poorly constrained (Reference Cook, Zwinger, Rutt, O’Neel and MurrayCook and others, 2012). In particular, future work could investigate whether hydrofracture influences glacier behaviour via its effect on calving events at the terminus and/or by deepening crevasses up to a few kilometres inland, particularly on floating sections, thus creating weaknesses that promote calving once the ice reaches the terminus.
The influence of submarine melt rates on HG’s frontal position is not modelled because there are very limited subsurface ocean temperature data from the fjord. However, due to its location at the extreme northwest of the GrIS, the contribution of submarine melting to mass loss at HG is likely to be much smaller than for other, more southerly glaciers where retreat has been linked to incursion of warm Atlantic Water (AW) (Reference Holland, Thomas, De Young, Ribergaard and LyberthHolland and others, 2008; Reference ChristoffersenChristoffersen and others, 2011). This is supported by measurements from the neighbouring Petermann Glacier fjord (Fig. 1), which recorded Atlantic-derived water from the Arctic with a maximum temperature of 0.2°C in 2009 (Reference Johnson, Münchow, Falkner and MellingJohnson and others, 2011; Reference StraneoStraneo and others, 2012). This is considerably lower than values reported for other major Greenland outlet glaciers, including Jakobshavn Isbræ (1.7–3.3°C) (Reference Holland, Thomas, De Young, Ribergaard and LyberthHolland and others, 2008; Reference Motyka, Truffer, Fahnestock, Mortensen, Rysgaard and HowatMotyka and others, 2011; Reference StraneoStraneo and others, 2012), Kangerdlugssuaq Glacier (1–1.8°C) (Reference ChristoffersenChristoffersen and others, 2011; Reference StraneoStraneo and others, 2012), Helheim Glacier (3.5–4.5°C) (Reference StraneoStraneo and others, 2012), 79 Glacier (1.0°C) (Reference StraneoStraneo and others, 2012) and central west Greenland (max 1.2°C) (Reference Rignot, Velicogna, Van den Broeke, Monaghan and LenaertsRignot and others, 2011).
Mass conservation calculations at Petermann Glacier suggest that basal melting is a major component of mass loss, despite the comparatively low ocean temperatures (Reference Rignot and SteffenRignot and Steffen, 2008). However, calving losses from Petermann Glacier were estimated at 0.59 km3 a−1, (Reference HigginsHiggins, 1991; Reference Rignot and SteffenRignot and Steffen, 2008), and applying the same calculation to the northern section of HG gives a calving rate five times greater than at Petermann Glacier (3.05 km3 a−1).This suggests that the relative contribution of submarine melting to ice loss at HG is much lower than at Petermann Glacier, due to higher calving rates, and that the role of submarine melting may be limited due to HG’s oceanographic setting. Furthermore, evidence from satellite imagery suggests that subglacial meltwater plumes are present during summer (Fig. 6a), particularly in the northern section of HG. This is not yet incorporated into numerical models, so perturbation experiments where submarine melt rates are simply increased are unlikely to be realistic. We therefore highlight the acquisition of subsurface ocean temperature data and the inclusion of plume flow in numerical models as a key area for future work.
6. Future Outlook
The basal trough beneath HG has important implications for its response to 21st-century climate change. The trough extends up to 72 km inland and has sections where the bed slope deepens up-glacier (Fig. 2). Simple linear regression of time versus frontal position on the northern section (1999–2012) gives a retreat rate of 427 m a−1 (R 2 = 0.95). If this trend continues, the terminus will remain in the trough for ∼170 years, highlighting the potential for sustained ice loss from HG during the 21st century and beyond, notwithstanding any nonlinear behaviour that could produce even higher retreat rates. One potential area where nonlinear behaviour could occur is close to the northern margin, where the location of the front on a pinning point may be currently causing deceleration (Fig. 2a and c). However, the bed deepens inland, suggesting that terminus recession from this pinning point may initiate rapid retreat and ice acceleration in this area and may produce substantial changes in the spatial pattern of velocity on HG.
It is likely that many other outlet glaciers are situated in similar overdeepenings to HG (Reference Cook and SwiftCook and Swift, 2012; Reference Morlighem, Rignot, Mouginot, Seroussi and LarourMorlighem and others, 2014), which could promote enhanced dynamic loss. Note, however, that the northern section of HG is substantially wider (∼45 km) than the majority of Arctic outlet glacier fjords (∼5 km), where the relative influence of lateral stresses on dynamics is likely to be greater (Reference RaymondRaymond, 1996). The influence of basal topography should therefore be considered in combination with fjord width, when interpreting and/or forecasting outlet glacier change in other locations (Reference Gudmundsson, Krug, Durand, Favier and GagliardiniGudmundsson and others, 2012; Reference JamiesonJamieson and others, 2012; Reference Enderlin, Howat and VieliEnderlin and others, 2013; Reference Carr, Stokes and VieliCarr and others, 2014).
7. Conclusions
Following almost three decades of gradual recession (37 m a−1),retreat at Humboldt Glacier accelerated markedly after 1999 (162 m a−1), coincident with significant increases in mean annual and mean summer air temperatures (Fig. 3a–d). Sea-ice decline may also influence retreat, but its role is not straightforward, because calving can occur when the sea ice is present and the relationship between sea ice and seasonal retreat rates is variable (Fig. 4). Subglacial topography strongly modulated HG’s behaviour and generated retreat rates that were an order of magnitude greater in its northern section, which is underlain by a major basal trough that extends up to 72 km inland (Fig. 2). Ice velocities were consistently higher in the northern section, but showed a dichotomous pattern of temporal change, potentially due to the presence of a basal pinning point near the northern lateral margin. Numerical modelling sensitivity experiments demonstrate that this differing response persists for moderate increases in crevasse water depth (as a proxy for atmospheric warming) and reduced sea-ice buttressing, simulated by front-stress perturbations (Fig. 5). We conclude that basal topography is an important influence on contemporary outlet glacier behaviour at seasonal to decadal timescales. A uniform basal topography cannot be assumed for individual glaciers, and these variations need to be considered when assessing glacier dynamics and response to forcing. At Humboldt Glacier, the basal trough is likely to continue to facilitate high retreat rates and mass loss during the 21st century.
Acknowledgements
Funding for this work was provided by a Durham Doctoral Studentship to J.R.C. Radio-echo sounding data were acquired and processed through UK Natural Environment Research Council (NERC) grant NE/H020667 to J.A.D. and P.C. and a G. Unger Vetlesen grant to the University of Texas Institute for Geophysics (UTIG). GrOGG laser altimetry was supported by NNXAD33G to D.D.B. This paper is UTIG contribution No. 2733. S.S.R.J. was supported by UK NERC fellowship NE/J018333/1.
