Asynchronous behavior of outlet glaciers feeding Godthåbsfjord (Nuup Kangerlua) and the triggering of Narsap Sermia's retreat in SW Greenland

Abstract

We assess ice loss and velocity changes between 1985 and 2014 of three tidewater and five-land terminating glaciers in Godthåbsfjord (Nuup Kangerlua), Greenland. Glacier thinning accounted for 43.8 ± 0.2 km³ of ice loss, equivalent to 0.10 mm eustatic sea-level rise. An additional 3.5 ± 0.3 km³ was lost to the calving retreats of Kangiata Nunaata Sermia (KNS) and Narsap Sermia (NS), two tidewater glaciers that exhibited asynchronous behavior over the study period. KNS has retreated 22 km from its Little Ice Age (LIA) maximum (1761 AD), of which 0.8 km since 1985. KNS has stabilized in shallow water, but seasonally advects a 2 km long floating tongue. In contrast, NS began retreating from its LIA moraine in 2004–06 (0.6 km), re-stabilized, then retreated 3.3 km during 2010–14 into an over-deepened basin. Velocities at KNS ranged 5–6 km a⁻¹, while at NS they increased from 1.5 to 5.5 km a⁻¹ between 2004 and 2014. We present comprehensive analyses of glacier thinning, runoff, surface mass balance, ocean conditions, submarine melting, bed topography, ice mélange and conclude that the 2010–14 NS retreat was triggered by a combination of factors but primarily by an increase in submarine melting.

1. INTRODUCTION

Glaciers along the margin of the Greenland ice sheet (GrIS) have been thinning and retreating over the past several decades (Warren, 1991; Warren and Glasser, 1992; Moon and Joughin, 2008; Leclercq and others, 2012) with tidewater outlet glaciers accounting for much of the ice loss: up to 58% before 2005 (e.g. Rignot and Kanagaratnam, 2006; Van den Broeke and others, 2009; Rignot and others, 2011; Shepherd and others, 2012) and 32% between 2009 and 2012 (Enderlin and others, 2014). Glaciers in southwest Greenland are among those that are experiencing significant ice loss, including those located in Godthåbsfjord (Nuup Kangerdlua), a 200 km long fjord near Nuuk (Fig. 1). However, the timing and magnitude of these changes has varied by glacier type and location. Three tidal outlet glaciers: Kangiata Nunaata Sermia (KNS), Akullerssuup Sermia (AS) and Narsap Sermia (NS), and three land-terminating glaciers: Saqqap Sermia (SS), Kangilinnguata Sermia (KS) and Qamanaarsuup Sermia (QS) drain into inner Godthåbsfjord, also known as Kangersuneq; two additional land-terminating glaciers are located near KNS and drain into a side branch of Godthåbsfjord: Isvand (IL) and Kangaasaruup Sermia (KSS).

Fig. 1. Drainage basins and 2009 surface velocities in the Nuuk region of Greenland overlain on July and August 2014 Landsat imagery. Drainages adopted from Van As and others (2014). Black lines indicate ATM flight paths over Kangiata Nunaata Sermia (KNS), Narsap Sermia (NS), and Saqqap Sermia (SS).

All three tidewater glaciers have experienced thinning and acceleration during the past two decades, similar to tidewater outlet glaciers elsewhere along the GrIS margin (Rignot and Kanagaratnam, 2006; Joughin and others, 2010a). The KNS and AS branches, and land-terminating QS, were joined during the LIA but now constitute separate glaciers as a result of post-LIA retreat (Weidick and others, 2012; Lea and others, 2014a). KNS is the major outlet glacier of the two and herein we will refer primarily to it when making comparisons with NS. KNS and NS have displayed a marked asynchronicity in behavior. KNS was at its LIA maximum in 1761, retreated ~5 km by 1808 and ~22 km by 2012 (Weidick and others, 2012; Lea and others, 2014a). In contrast, NS has just begun retreating from its LIA maximum. Roughly 40 km apart, the two outlet glaciers would seemingly be affected by similar atmospheric and oceanic conditions. This asynchronous behavior in part motivated the present study. Similar asynchronous behavior is well documented in coastal Alaska (e.g., Post and others, 2011; McNabb and Hock, 2014; Truffer and Motyka, 2016) as well as for other GrIS outlet glaciers (e.g., Carr and others, 2013; Bartholomaus and others, 2016).

A number of studies have addressed changing conditions that could impact Godthåbsfjord glaciers. Several have examined ocean conditions (Mortensen and others, 2011, 2013, 2014; Bendtsen and others, 2015a) while van As and others (2014) and Langen and others (2015) used regional climate models to estimate freshwater runoff into the fjord and surface mass balance (SMB). Weidick and others (2012) and Lea and others (2014a, b) discussed the KNS LIA boundary and its post LIA retreat while Lea and others (2014b) modeled the influence of atmospheric forcing on KNS terminus retreat.

The general understanding of tidewater glacier dynamics, once retreat has been initiated, has substantially improved due to an increase in detailed observations (e.g., Walter and others, 2010) and modeling (e.g., Pfeffer, 2007; Vieli and Nick, 2011); however, the mechanisms that initiate retreat remain elusive. Several processes have been suggested including atmospheric warming, oceanic warming and changes in ice mélange conditions (Straneo and others, 2013; Moon and others, 2015; Truffer and Motyka, 2016). Our goal is to document changes occurring at all of the glaciers in Godthåbsfjord and then scrutinize drivers that may have forced these changes. We examine and synthesize a diverse set of data including: (1) DEMs and laser altimetry to determine ice loss, thinning rates and geodetic balance; (2) terminus positions from satellite images and aerial photos; and (3) velocities from satellite radar and optical imagery. We evaluate frontal ablation by calculating terminus ice fluxes and rates of terminus volume change, and examine effects of atmospheric forcing using the HIRHAM5 climate model (Langen and others, 2015). We also investigate ocean conditions in the fjord using previously published data, as well as our own fjord hydrographic and bathymetric data and fjord-surface temperatures (FST) derived from MODIS images as a proxy for ice mélange conditions.

2. DATA AND METHODS

2.1. Digital elevation models

2.1.1. 1985 DEM

The Agency for Data Supply and Efficiency (the successor agency of the Geodetic Institute) recorded aerial stereo-photography covering Godthåbsfjord on 19 July 1985 at a scale of 1:150 000. We produced a 25 m DEM for 1985 from these photos that covers ~1800 km² of the GrIS margin in our study area. Lack of ground control and air photo coverage limited the extent of upglacier coverage. The DEM was produced following methods in Korsgaard and others (2016). The DEM was regridded to 40 m and referenced to the WGS 84 ellipsoid for comparison to SPOT DEMs and laser altimetry.

2.1.2. 2008 SPOT-5 DEM

We obtained SPOT5 imagery and associated DEMs of the region for 22 June 2008 and 02 August 2008 under the SPIRIT Polar-Dali Program (Korona and others, 2009). Ground resolution is 5 m along track and 10 m across track. The DEM has a grid spacing of 40 m and referenced to the WGS 84 ellipsoid. We used the 02 August 2008 imagery, which covers almost the entire study area, including KNS and most of NS, for our analysis and merged a portion of the earlier DEM to fill in areas in the north part of our region. We also masked minor areas of cloud cover and excluded these areas from our analysis.

2.1.3. 2014 Worldview DEM

We obtained Worldview stereoscopic submeter-resolution imagery of our study area from DigitalGlobe Inc. and distributed by the Polar Geospatial Center at the University of Minnesota. The image pairs were acquired on 04 August 2014, 05 August 2014 and 27 July 2014. The August 2014 imagery covers most of the study area while the July pair expands coverage to the north to include SS. We used the Ohio State University DEM extraction software Surface Extraction through TIN-Based Searchspace Minimization (SETSM, Noh and Howat, 2015) to construct the DEM and generate orthoimages. SETSM is fully automated and requires no user input other than the imagery and its accompanying metadata. The DEMs were created at 2 m posting referenced to the WGS 84 ellipsoid and then down-sampled to 40 m through bilinear interpolation to facilitate comparison with our other two DEMs.

2.1.4. DEM accuracy

We assessed DEM accuracy using a variety of ground-truth datasets over land areas, including our own kinematic and static GPS surveys near KNS, Atmospheric Topographic Mapper (ATM) laser altimetry data from August 2008 and April 2014, and May 2011 Land, Vegetation and Ice Sensor (LVIS) data. The latter two datasets provided excellent spatial coverage for our region. We excluded land terrain where we knew the DEMs failed to model the land surface and filtered LiDAR data for slope roughness.

For the 1985 DEM, the results of comparing elevations to the ground-truth datasets at ~21 500 land points showed a Gaussian distribution with a standard deviation (σ _sd ) = 4.5 m and a bias of −1.7 m. Given these results and the close proximity of the points to the ice margins (<5 km), we adjusted the DEM for the bias and assign an error of ± 4.5 m per pixel for the DEM over ice. The error and bias are similar to the values reported by Kjær and others (2012), Kjeldsen and others (2015) and Korsgaard and others (2016), based on DEM co-registration to ICESat data for this region.

Korona and others (2009) reported the accuracy of GrIS SPOT DEMs to be within ±6 m of ICESat for 90% of the data. We used our land based datasets (GPS, ATM and LVIS) to further assess the accuracy of the SPOT DEM for our region of interest. This comparison (19 300 points) showed a Gaussian distribution with a bias of 2.2 m, and σ _sd  = 5.3 m. Serendipitously, the 02 August 2008 ATM flight over KNS coincides with the SPOT 5 acquisition date. These data were used to independently assess the accuracy of the 2008 DEM over KNS ice. The bias of 2.1 m is similar to that found for the land points but the standard deviation (2.3 m) is much lower. We adopt the latter standard deviation for point elevation uncertainty after adjusting for bias.

Using similar methods for the 2014 DEM, we obtained a σ _sd  = 8.0 m with a bias of 0.86 m from the analysis of ~21 300 land points. We also used additional 2014 ATM data over ice for glaciers SS, QS and KSS to check the 2014 DEM. Because of dynamic thinning during the nearly 4-month time difference between the ATM flight and the Worldview imagery, we excluded ATM data over the tidewater glaciers for DEM validation. A correction for ablation was obtained by examining ablation data at three weather stations maintained by GEUS under the PROMICE program (www.promice.org): two are located on QS and the third on KS. We neglect summer emergence velocity. The results from comparing ~6200 ATM ice points gave a σ _sd  = 2.9 m and a bias of −2.0 m, indicating a substantially better accuracy of the 2014 DEM over ice than over land. We adjusted for the bias and adopted the latter standard deviation as a measure of elevation accuracy of ice grid points. We note that Noh and Howat (2015) reported a much smaller error, comparable in accuracy with the LiDAR data, for the original 2 m DEM; however, scaling the DEM to 40 m posting decreased accuracy.

2.1.5. DEM differencing

After correction for errors and biases we differenced the three DEMs to produce elevation change (dZ) maps with grid spacing of 40 m for 1985–2008 and 2008–14. To facilitate comparison between the two time periods, we used the area of overlap for the three DEMs to estimate volume changes (ΔV) and the area averaged water equivalents (ΔZ _ave ) for each glacier and over each period. We neglect any isostatic uplift and assume no changes in bed elevation, such as may be caused by erosion and sediment deposition; these effects are considered minor, likely much less than a few cm per year (Hallet and others, 1996).

To estimate the uncertainty (or standard error) of such calculations when comparing DEMs to determine volume change, we adopt the methodology developed by Rolstad and others (2009) and outlined in Motyka and others (2010). This method uses variograms of the differenced DEMs over adjacent land areas to determine an area of correlation, Ac, which is then taken as a measure of error correlation between the two DEMs over the ice. When comparing the 85–08 and 08–14 DEMs, we found Ac = 0. 8 km² for both, which is considerably smaller than the area, A, for both the entire region of coverage (~2100 km²), and for the individual glacier drainages.

One further adjustment is commonly examined in geodetic calculations: changes in ice and firn density profiles between DEM dates. In our case, the snow line is well above the highest 1985 DEM elevation so that elevation changes are entirely those of an ice surface with a single density, which we assume to be 910 kg m⁻³.

2.2. Airborne laser altimetry

2.2.1. Airborne topographic mapper

Table 1 lists NASA ATM data for KNS and NS (Krabill, 2014). Track lines are shown in Figure 1. The KNS centerline was first profiled by NASA in 1993, then again in 1998 and 2001. Starting in 2008, ATM surveys were conducted on an approximately annual basis (except for 2013) over the centerlines of both KNS and NS. Coverage was extended to SS in 2012. Data span the last two decades for KNS and frequency of acquisition in later years allows investigation of short-term variations in KNS surface elevation. We used the ATM ICESSN product, which fits the swath of laser shots with tiles both across swath and along track – typically a tile represents 500–1000 laser shots fit to a plane that is 40–70 m on a side. The early ATM flights had a narrower track and were commonly flown over or near the centerline of a glacier. However, for KNS, the flight path did not necessarily follow the main flow lines. KNS ATM tracks for 2001 and earlier consisted of one or two lines of tiles fit to the swath of laser returns. In 2008, the swath width was increased yielding three sets of tiles of a similar size across the swath. NASA switched to swath LiDAR in 2010 with returns covering a width of ~300 m. The reported accuracy of ATM data is ±0.3 m.

Table 1. Inventory of NASA ATM flights and DEMs

Our ATM elevation change comparisons use the 1985 DEM as base reference. At KNS, we use the 2008 DEM rather than the 2008 ATM data to calculate ΔZ (surface elevation change) because the 2008 ATM flight line diverged from the path used in other years.

2.2.2. 2011 LVIS

Land, Vegetation, Ice Sensor (LVIS) data of the region was flown in May 2011. Swath width averaged 400 m. Because the flight pattern of LVIS did not follow glacier flowlines or the paths of ATM flights, we elected not to include the LVIS data in our analysis. However, these data did prove useful for ground-truthing the three DEMs.

2.3. Terminus positions

Time series of terminus positions were generated using satellite images and aerial photos for the time span 1985–2015. The bulk of the record was derived using the 15 m panchromatic scenes (Band 8) from Landsat 7 and 8, and in early years using 30 m Thematic Mapper scenes (Band 3) onboard Landsat 5. The record is supplemented with a handful of 60 m Multispectral Scanner images (Bands 2 and 3) from Landsat 4 and 5, and a few 15 m ASTER images (Bands 2 and 3). Landsat scenes were downloaded from the USGS Earth Resources Observation and Science (EROS) Center's Global Visualization Viewer (GLOVIS), ASTER images were obtained from NASA's EOSDIS (Land Processes Distributed Active Archive Center (LP DAAC), 2001) and the 1985 position was derived using an aerial photo. All Landsat and ASTER images and the 1985 photos were orthorectified and used the same ellipsoid (WGS84). Glacier termini were manually digitized and mean front positions calculated following the techniques of Cassotto and others (2015). Center flowlines were determined using ice velocity fields. Points within 1.5 km of KNS and NS flowlines (total width = 3 km; colored points in Supplementary Material Figs S1a and S1b) and within 0.8 km of SS flowline (1.6 km total width) were used to calculate mean front positions; points along slow moving margins (gray points in Figs S1a and S1b) were excluded. These positions were then plotted as a time series to evaluate seasonal and decadal variations in ice-front behavior.

2.4. Ice velocities and frontal ablation

2.4.1. Ice velocity

Winter ice velocity fields derived from InSAR for 2000/01 and for 2005/06 through 2008/09 were obtained from the NASA MEaSUREs dataset (http://nsidc.org/data/NSIDC-0481/versions/1/; Joughin and others, 2010b). Annual coverage for KNS was extended through 2012 using TerraSARX images but unfortunately TerraSARX did not cover NS. Estimated uncertainty of these SAR-derived velocities is 6%. For details on deriving these data from InSAR and TerraSAR images, the reader is referred to Joughin and others (2010a).

We extended the SAR dataset by deriving additional velocity fields for winters of 2013–15 for KNS and NS from Landsat 8 images. Velocities were determined by using feature tracking (Fahnestock and others, 2015); uncertainties are estimated at 3%. In addition, we used available Landsat 7 and earlier images to fill in gaps in the time series for both glaciers covering the period 1987–2012, albeit with some major gaps and with estimated uncertainties of 5%.

2.4.2. Glacier depths, ice flux and frontal ablation

We use KNS and NS bed depths derived by Morlighem and others (2014) to help calculate terminus ice fluxes and to investigate the potential control of glacier bed topography on future activity. The data were obtained from the Operation IceBridge (OIB) Earth Science DataSet at the National Snow and Ice Data Center (NSIDC) (Morlighem and others, 2015). These bed depths are based on CReSIS radio echo sounding data and the “conservation of mass” (or MC) method to evaluate glacier thickness (Morlighem and others, 2011). Estimates of error accompany the MC bed data and tend to be lowest along centerlines near the terminus (~20 m), and to increase towards glacier margins and further upglacier (50–100 m). However, we compared our bathymetry data obtained near the KNS and NS termini (discussed later), with the MC determined terminus depths, and discovered significant discrepancies with the MC data showing much shallower beds than our soundings. We therefore modified the MC terminus data by using our soundings as a control. This comparison also indicated that the published error estimates for the MC bed are too optimistic.

We evaluate frontal ablation, Q _fa , (which includes both calving and submarine melting) from the difference between ice flux arriving at the calving front, Q _i and the volume change at the terminus (advance/retreat):

(1) $$Q_{\,fa} = Q_i - {\rm d}V/{\rm d}t$$

where the volumetric rate of retreat is dV/dt (e.g., O'Neel and others, 2003).

To determine Q _i , we use “flux-gates” located ~3 and 4 km from the 2008 KNS and NS termini, respectively. The choice of gate location is dictated by the quality of bed and velocity data, and in the case of NS, 2014 terminus position. We use our three DEMs to establish flux-gate glacier surface elevations, Z _i (x,y) and then use dZ _i determined from ATM data to interpolate and adjust elevations for the years between our DEMs. We obtain flux-gate ice thickness H(x,y) by subtracting bed elevations Z _b(x,y), from Z _i (x,y). Ice flux through the gate, Q _g, is

(2) $$Q_{\rm g} = { \vint} H(y)\; U(y)\; {\rm d}y $$

where y is the distance along the gate, H is the ice thickness and U is ice velocity perpendicular to the gate obtained from the SAR and Landsat datasets. We assume that the surface velocity is entirely due to sliding (i.e., plug flow), a reasonable assumption for regions close to the terminus. For consistency, we utilize winter to early spring velocity fields derived from SAR and Landsat to determine gate-normal velocities. For years where only centerline velocities are available, we examine velocity fields derived from the year closest in time and then apply a linear adjustment based on comparison of centerline velocities to estimate velocities along the width of the flux gate.

Two further adjustments are needed to account for glacier thickness change, dH/dt, and for surface ablation, B, between our flux-gates and the glacier termini. Based on PROMICE data for this region (http://promice.org) we assume an average annual ablation of 5 m a⁻¹ in the terminus region and we determine dH/dt from analyses of the ATM record. We then integrate across the area below the flux gate to obtain the correction Q _b:

(3) $$Q_{\rm b} = \vint \limits_\alpha ^{} \left( {B - {\rm d}H/{\rm d}t} \right){\rm d}\alpha $$

where B is the SMB (negative for ablation) and α is the area between the gate and the terminus. Then

(4) $$Q_i = Q_{\rm g} + Q_{\rm b} $$

We then use our modified terminus bed model, the DEMs and the analysis of terminus change to compute dV/dt:

(5) $${\rm d}V/{\rm d}t = \; \mathop \vint \limits_A^{} H_i \; {\rm d}A$$

where A is the areal extent of retreat or advance and H _i is the ice thickness.

Uncertainties for Q _g accrue from our estimates of flux gate ice thickness, H(y), and gate-perpendicular ice velocities, U(y). We use ±30 m instead of the published bed error map (Morlighem and others, 2014). The uncertainty of our surface elevations varies from 4.5 m for 1985 to 2.3 m for 2008 and 2.9 m for 2014. For intermediate years, additional uncertainty is introduced by using ATM data to extrapolate surface elevations. We estimate the latter to be of the order of ±2 m. For surface velocities, uncertainties range from 3% to 6% depending on imagery used. Additional uncertainty is introduced by our assumption of plug flow, although we consider this to be of minor significance. A more serious consideration is the seasonal fluctuation in velocity and how representative are the velocities that we adopted for the annual average. We assign an additional 10% uncertainty to account for these uncertainties. Additional sources of uncertainty when computing Q _i come from Q _b, with uncertainties in the area α, the estimate of surface ablation (±1 m a⁻¹) and analyses of dH/dt. Uncertainties in digitized terminus margins are ±30 m for earlier Landsat and ±15 m for Landsat 7 and 8, resulting in uncertainties in α ranging from 1 to 2%. Sources of uncertainties in dV/dt include uncertainties in H and in assessing the area of retreat, A, (1–2%).

The DEM differencing captures the ice loss due to changes in ice-surface elevations (and this loss directly contributes to sea-level rise), but additional ice losses accrue from below-sea-level portions of the tidewater glacier termini due to terminus retreat (losses that do not contribute to sea-level rise). We exclude this ice loss from sea-level-rise because these regions were grounded and became inundated with sea water as the terminus retreated. We account for the loss of below-sea-level ice using a modified form of Eqn (5): by integrating the submerged ice thickness, H _s, instead of the total ice thickness, H _i , over the area of retreat. We determine H _s from our bathymetry and from Morlighem and others (2015) ice depths, where we assume 10% accuracy in terminus bed depths.

2.5. Runoff and SMB

We derive estimates of freshwater runoff and SMB for each glacier drainage using the HIRHAM5 regional climate model (RCM) data for 1989–2015. Divides between glacier drainages (Fig. 1) are adopted from van As and others (2014) and were determined using a Geographic Information Systems particle tracking tool and the 2006 velocity field (Joughin and others, 2010a). Details of applying the HIRHAM5 model to glaciers in the Godthåbsfjord region are discussed in Langen and others (2015). The HIRHAM5 for Greenland has a resolution of ~5.5 km. Analyses by Ettema and others (2009) and Lucas-Picher and others (2012) suggest that such resolution is necessary to resolve the coastal topography that impact variables such as near-surface air temperature and precipitation. However, the current HIRHAM5 model uses surface elevation of the ice sheet from Bamber and others (2001), which has high accuracy (within a few meters) in the flat ice-sheet interior, but suffers from significant error near the coast.

Langen and others (2015) compared HIRHAM5 with PROMICE data from stations located on glacier QS and concluded that the model may overestimate albedo, resulting in potentially underestimating melt below 1000 m by ~20% and thus overestimating SMB. Accumulation biases are within 10% compared with ice core-derived accumulation rates in the region (Lucas-Picher and others, 2012).

2.6. Ocean hydrography

We evaluate the potential influence of submarine melting on frontal ablation and terminus stability by analyzing hydrographic data in the vicinity of the calving front. Unfortunately, inner Kangersuneq is frequently choked with icebergs that can inhibit access by sea vessel; thus, data are sparse for the inner parts of this fjord system. However, we were able to obtain a detailed time series of hydrographic data within Kangersuneq that spans the period April 2008 to December 2014. Station locations varied because of ice conditions (see Mortensen and others (2013, their Fig. 1, for locations)); however the overwhelming majority were located just west of NS Bay with many between NS Bay and the KNS LIA sill. In addition, we were able to acquire hydrographic data from shipboard surveys within 12 km of KNS and within 5 km of the NS calving terminus during August 2011, and depth soundings close to the terminus of NS during a reconnaissance visit in July 2008. Additional spot soundings near the KNS terminus in 2010 and 2011 were obtained by using a helicopter and lowering a transducer through breaks in the ice mélange. These soundings near the terminus provide estimates of ice thickness at the glacier faces and help control MC-derived terminus glacier bed data. Accuracy of both the helicopter and shipboard soundings is estimated to be ±2 m.

Details from our August 2011 hydrographic transects and the 2008–14 time series for Kangersuneq will be subjects of separate papers. For this study, we draw upon hydrographic results that document inner fjord water depths, temperature and salinity. At each station, we lowered a calibrated SeaBird Electronics SeaCAT 19 “plus” conductivity-temperature-depth instrument and used standard procedures to measure temperature and salinity of the water column and averaged them into 1 m bins. For additional data on fjord and coastal temperature and salinity we draw upon Mortensen and others (2011, 2013) and Ribergaard (2014).

2.7. Fjord surface temperature - ice mélange proxy

An ice mélange is a mixture of icebergs and sea ice found in many proglacial fjords. Fjord constrictions and shallow sills often promote ice mélange accretion (Peters and others, 2015). In Kangersuneq, an ice mélange frequently extends to the LIA moraine ~22 km from the present KNS terminus (Landsat images cf. Supplementary Material Video), but in the winter the mélange can extend appreciably beyond the moraine, up to 40–60 km or more downfjord. In order to assess the impact of mélange variability, we used sea surface temperatures (SST) derived from the thermal bands of MODIS as a proxy for ice mélange conditions following the methods of Cassotto and others (2015). In general, the proxy interprets cold surface temperatures as cohesive ice mélange that exhibits rigid properties, while warmer temperatures represent increased mobility of a granular-like mélange. More than 5400 daily, 11 µm, level-2 Terra granules spanning >15 years were obtained from NASA's Ocean Data Processing System (Feldman and McClain, 2015). The granules were reprojected to 1 km grid spacing and filtered in time. The SST data product is based on brightness temperatures calculated from raw radiance values measured by the satellite. The emissivity of water is accounted for in the derivation; however, differences in the emissivities of water and ice are small (~0.01) and thus negligible. Therefore, the SST data product reflects the temperature in each pixel, which may contain seawater, sea ice, icebergs, or cloud tops. Herein, the record produced from the SST data product is referred to as the fjord surface temperature (FST) record to reflect variations in fjord surface – ice mélange conditions, and to avoid confusion with real water temperatures at the sea surface layer. The majority of noise in the record is derived from cloud top surfaces; however, to maintain a continuous record of fjord surface conditions, cloud temperatures were not removed. Instead, the data were filtered to reduce the high frequency, very cold temperatures (e.g. =−40°C) related to cloud tops. The filter ignores warm temperatures related to low strata or fog in the fjords; however, our assessment is based on winter mélange conditions when ambient air temperatures are typically low and below freezing in the arctic. Kangersuneq's fjord width ranges from 3 to 4 km, and the FST data product has a native resolution of 1 km. Inspection of time-lapse photographs and satellite images, when available, show little across-fjord variability in winter mélange conditions. Therefore, a 100 km long centerline profile was sampled and compared with variations in the terminus locations for KNS and NS.

2.8. Submarine melting

Submarine melting is a function of thermal forcing, TF, (dictated by ocean temperature, salinity and pressure) and subglacial discharge, Q _sg . Limited field data (e.g., Motyka and others, 2003, 2013) and modeling experiments (e.g., Jenkins, 2011) indicate that submarine melting, Q _m, is directly related to TF and to the 1/3 power of Q _sg . Thus, an increase or decrease in either quantity would increase or decrease Q _m. Unfortunately, our data are insufficient to accurately assess Q _m. Instead, we used parameters developed from models (e.g., Jenkins, 2011; Xu and others, 2013) to derive estimates of Q _m for NS and KNS in order to examine the sensitivity of Q _fa to changes in Q _m. For this exercise, we assumed Q _sg is equal to runoff as predicted by the HIRHAM5 model for KNS and NS. To determine TF we computed the average of temperature and salinity at depths of 120 m and 150 m determined from our hydrographic measurements in Kangersuneq. The data were then interpolated for each monthly Q _sg time increment between April 2008 and December 2014. We also made adjustments for the vertical area of the calving face for KNS and for NS as each retreated into deeper water.

3. RESULTS

3.1. Terminus change

The time series of mean terminus positions for KNS and NS show very different patterns of retreat (Fig. 2). KNS exhibits large seasonal variations (up to 2 km), but the annual minimum position has retreated only 0.8 km since 1985 (Fig. 2a). The summer position retreated 300 m between 1985 and 1999 and then remained relatively stationary until 2004 when it retreated ~300 m and another ~300 m by 2005. Since 2005 the summer position has oscillated, retreating some years and advancing in others. The net change in summer position between 1985 and 2015 (0.8 km) represents ~3.6% of the ~22 km retreat from its 1761 AD LIA maximum (Lea and others, 2014a).

Fig. 2. Plot of average terminus position vs time for KNS (a) and NS (b). Data points are color-coded with time, using the same color-coding in later figures of time series. Black diamond in (a) marks July 1985 terminus position of KNS, 300 m beyond its 1999 summer position. Dashed line in (a) traces summer retracted position. Inset in (b) extends record for NS back to 1985, which shows terminus was stable until 2005.

In contrast, NS was generally stable with minimal seasonal fluctuations before a large, episodic retreat began in 2004, with NS retreating a total of ~4 km by 2015 (Fig. 2b). The terminus remained at or near its LIA maximum position until 2004 (Fig. 2b). Prior to 2004, NS experienced seasonal oscillations of ~100–200 m, and did not appear to have any significant seasonal floating tongue. NS then underwent a retreat of ~600 m between 2004 and 2006. The glacier restabilized until the summer of 2010 when a second much stronger retreat began and continued into the summer of 2013, totaling ~3.3 km or ~1.1 km a⁻¹. Both pulses of retreat were marked by the development of large embayments into the terminus, particularly in the center and north side (cf. Fig. S1b and video in Supplementary Material). The most significant changes occurred in late 2010 when a large section of the northern terminus broke off and disappeared by the winter of 2011 (Fig. 2b and Fig. S1b and video in Supplementary Material). The NS terminus has experienced seasonal oscillations of ~500 m meters since the summer of 2013 and the retracted summer terminus retreated only slightly between 2013 and 2015 (~300 m).

Tidewater glaciers entering Kangersuneq are all presently grounded, although we attribute the large seasonal oscillations at KNS (and NS in later years) to the seasonal development of a floating tongue, similar to observations at Jakobshavn Isbrae (Amundson and others, 2008; Joughin and others, 2008; Cassotto and others, 2015). The early Landsat record shows this tongue sometimes coalesced with the AS terminus during the late 1980s (Fig. S1a and video in Supplementary Material). Advance typically begins in early fall, reaches a maximum in late winter to early spring and is followed by a rapid retreat to a late summer minimum (Fig. 2a). Noteworthy, the KNS floating tongue was significantly shorter during the winters of 2011 and 2015, and oscillations have become less pronounced in most recent years.

Most land-terminating glaciers in the region have remained at their LIA extent with total changes <100 m. QS is an exception, having sustained a 2 km post-LIA retreat by 1985 after the breakup and retreat of KNS (Weidick and others, 2012) but the glacier has been relatively stable since. Another exception is SS located north of NS (Fig. 1), which, based on air photos and satellite imagery, advanced slowly during the 20th century and 250 m between 1985 and 2009. The terminus showed minor seasonal oscillations (< 50 m) until 2013/14, when the summer position retreated ~80 m (cf. Fig. S2 in Supplementary Material).

3.2. Elevation and ice volume changes

All glaciers in the study region experienced significant drawdown between 1985 and 2014 (Figs 3 and 4, Table 2); however, the data show asynchronous patterns that emerge by region and glacier type. Table 2 also provides the uncertainty in volume change, σ_ΔV , and in ΔZ _ave , σ _we . The relatively small values for σ_ΔV as well as for σ _we are a direct consequence of using the correlated area analysis and a very large area of integration (A >> Ac).

Fig. 3. Total elevation change in meters between 1985 and 2008 DEMs (a) and between 2008 and 2014 DEMs (b) superimposed on shaded relief of 2008 DEM. Terrain contour intervals are 200 m. Red boundaries define the ice margin and glacier drainage divides: KSS, Kangaasaruup Sermia; IL, Isvand; KNS, Kangiata Nunaata Sermia; AS, Akullerssuup Sermia; NS, Narsap Sermia; QS, Qamanaarsuup Sermia; KS, Kangilinnguata Sermia; SS, Saqqap Sermia.

Fig. 4. ATM elevation change time series for (a) KNS, (b) NS and (c) SS referenced to the 1985 DEM as base. Color codes used to identify dates are similar to those used in Figure 2. Note Y-axis scale difference for (c) SS. HAE = “height above ellipsoid”. Panel (d) provides the elevation change per year between ATM dates averaged over the 100–200 m surface elevations at KNS (blue) and NS (red).

Table 2. Ice volume loss, ΔV, from different drainages as determined from DEM differencing over area A

Area-averaged water equivalent loss, ΔZ _ave also shown. Only areas common to all three DEMs were used to evaluate ΔŻ _ave . σ_ΔV and σ _we are uncertainties in volume change and area-averaged water equivalent loss, respectively.

Between 1985 and 2008, the greatest losses occurred predominately along southern glaciers (Fig. 3a, Table 2) with the KNS terminus thinning by 90 m or more. The ATM record (Fig. 4a) indicates most of this loss at KNS occurred between 2001 and 2008. In contrast, NS and northern land-terminating glaciers were the primary contributors to ice loss between 2008 and 2014, while losses diminished significantly along KNS and adjacent southern land-terminating glaciers (QS, KSS and IL). In total, the DEMs document 43.8 ± 0.2 km³ of ice loss between 1985 and 2014 over the surveyed region, equivalent to 0.10 mm eustatic sea-level rise. Below-sea-level ice losses from KNS and NS glacier retreat (Table 3), account for an additional 3.5 ± 0.3 km³ of ice loss for a total ice loss between 1985 and 2014 of 47.3 ± 0.4 km³. The largest of these below-sea-level ice losses is associated with the 2008–15 terminus retreat of NS and accounts for a large fraction of total ice loss from NS during this period (10.9 km³ vs 2.7 km³; Tables 2 and 3).

Table 3. Submerged ice volume loss, ΔV, between 1985 and 2008 and between 2008 and 2014 due to glacier retreat over area A. σ _ΔV are uncertainties in volume change

ATM data (Fig. 4) show that surface elevation changes have been quite variable at KNS. The glacier thinned a total of 10–20 m between 1985 and 1993, and then thickened to or above its 1985 elevation along most of its length by 1998. Thinning resumed after 1998 with greatest losses occurring between May 2001 and August 2008 (Figs 4a, d). KNS then thickened again between August 2008 and May 2009, although some of this change is a seasonal effect. KNS suffered some of its most significant ice losses between May 2010 and April 2011, dropping 30–40 m at the lower elevations (Fig. 4d), and a total of over 100 m compared with 1985. The glacier recovered somewhat in the succeeding year (April 2012) and elevations then remained constant into April 2014. Thus, the DEM computed ice losses at KNS for 2008–14 in Table 2 occurred primarily during 2010/11.

Figure 4b shows drawdown at NS was underway by 2008, 40 m since 1985 at lower elevations. Surface elevation then remained constant into May 2010 but thinning then continued unabated into April 2014 (Figs 4b, d). Drawdown at lower elevations compared with the 1985 surface reached a total 140–150 m with most of the thinning occurring after 2010.

Land-terminating glacier SS, subject to atmospheric changes but unaffected by terminus ocean processes, serves as a contrast to the tidewater glacier systems. Figure 4c shows that between 1985 and 2008 SS remained comparatively stable except for thickening at the lower elevations. However, SS began showing signs of thinning by 2012, by several meters compared with 2008.

3.3. Ice velocities

In general, winter centerline velocities at KNS were relatively constant between winter 2000/01 and February 2015 with some minor interannual oscillations (Fig. 5a). Velocities at the terminus ranged from ~5.5 km a⁻¹ in 2000/01 to a high of ~7 km a⁻¹ in February 2012 with a speedup of ~15% occurring between 2000/01 and 2006/07. Similar trends are apparent for the time series at 1.3 km upglacier (Fig. 5c): velocities ranged from 4 to 5 km a⁻¹ between 1987 and 2011, then increased to ~6 km a⁻¹ by 2012 but then declined back to 4–5 km a⁻¹ in 2013–15.

Fig. 5. Center line (mostly winter) velocities for KNS (a) and NS (b) derived from RADARSAT and Landsat. Color codes are same as used in Figures 2 and 4. The 2000/01 KNS and NS termini lie at the left end of panels (a) and (b). Centerline velocity time series (1987–2015) 1.3 km from the 2008 terminus of KNS and 5 km from the 2008 NS terminus are shown in (c) and (d), respectively. Symbols in (c) and (d) correspond to data source: + = Landsat 8, X = Landsat 7 and Y = Landsat 5. The black circles are from RADARSAT Mosaics (NS and KNS) and TerraSARX (KNS after 2009). Symbol color indicates date (cf. Figs 2, 4, and 5).

Velocity changes at NS have been more dramatic with speeds in the terminus region doubling between winters 2000/01 and 2005/06 from 1.5 to 3 km a⁻¹, and then nearly doubling again by February 2013 (Fig. 5b). Speeds accelerated all along the glacier and ranged from 5.5 km a⁻¹ at the terminus to 1.5 km a⁻¹ 20 km upglacier, similar in magnitude to KNS speeds. The time series at 5 km from the terminus (Fig. 5d) show velocity remained at ~1.5 km a⁻¹ from 1996 to the spring of 2002. Landsat imagery from spring 2002 until winter 2005/06 proved unsuitable for feature tracking, but images for 2005/06 confirmed an acceleration had taken place. The velocity time series from Landsat also show that there is considerable seasonal variation in near-terminus velocity at both KNS and NS during 2013 and 2014 (Figs 5c, d).

3.4. Fjord hydrography

Depth to the KNS LIA sill is ~160 m with the axial proglacial fjord deepening to over 400 m (Fig. 6a; cf also Fig. 9c). Depths in front of the 2011 KNS terminus range from 160 to 260 m. Depths over the NS LIA moraine are also ~−160 m then quickly deepen to −300 m outward from the moraine.

Fig. 6. Details of KNS (a) and NS (b) proglacial fjord and terminus on shaded relief images (August 2008 for KNS and July 1985 for NS). Land and ice contours (50 m) are in grey. Fluxgate locations are used for analyzing ice flux in Table 4. Red stars in both panels indicate shipboard soundings and conductivity-temperature-depth casts obtained August 2011. In panel (a) cyan crosses show locations of heli-soundings with depths above. LIA sill (black) and the 1985 terminus (red) are also shown. In panel (b), the blue contours are based on bathymetric soundings from July, 2008 (blue crosses show data points). Red dashed lines show NS terminus positions for August 02 2008 and August 2014.

Cross sections of temperature and salinity depth profiles acquired in August 2011 show highly stratified water columns near both glaciers (KNS: Fig. 7a; NS: Fig. 7b). Following water mass classifications of Mortensen and others (2011, 2013), we observed a ~50 m thick upper layer of cold (0–1°C), relatively fresh water (subglacial water) overlying an intermediate layer of warm (~3°C), more saline water (~33.3 PSU) extending from 50 m to ~200 m (sill region water), and then a third deeper layer of cooler (~2°C) saline water (~34 PSU) (basin water) extending to the bottom.

Fig. 7. Results of shipboard conductivity-temperature-depth casts obtained August 2011. (a) Cross section along fjord axis and over KNS LIA sill (Fig. 7a). The KNS terminus is to right. (b) Cross section across the mouth of NS Bay looking towards NS terminus (Fig. 7b). North is left.

The 2008–14 time series for temperature and salinity from depths of 120 and 150 m show considerable intra and interannual variation (Fig. 8). The data document the incursion of 3°C water beginning in fall of 2010. These elevated temperatures prevailed throughout 2011 into winter 2012. Temperatures then underwent seasonal oscillations, with coldest periods (~1.5°C) in late spring and warmest (~2.5–3°C) in late fall. In contrast, a colder regime prevailed prior to 2010, ranging from 0.7 to 2.5°C. Salinity fluctuated slightly, ranging from 32.6 to 33.8 PSU, with lower salinity roughly correlating with warmer temperatures (Fig. 8).

Fig. 8. Inner fjord water temperature (a) and salinity (b) obtained from Mortensen and others (2013) and Mortensen and others (unpublished data) at depths of 120 m (blue) and 150 m (red). Station locations varied but mostly occur between NS Bay and the KNS LIA sill.

3.5. KNS and NS Glacier bed topography

Bed geometries for the two tidewater glaciers are strikingly different (Figs 9a, b). Whereas KNS is currently retreating onto shallower bed topography (~−200 m), NS has retreated from its 160 m deep LIA moraine into a −300 m overdeepening. Upglacier, bed topography at both glaciers becomes increasingly shallow, rising to −75 m at KNS and −140 m at NS (Fig. 9c). Farther upglacier, the NS glacier bed dramatically deepens to −1000 m before emerging above sea level at ~40 km from the 2014 terminus position. At KNS, bed geometry reverses slope and deepens to −325 m at ~8 km from the 2014 terminus then emerges above sea level a few km further upglacier (Fig. 9c).

Fig. 9. Glacier bed topography for KNS (a) and NS (b) modified from Morlighem and others (2014; 2015) using near-terminus fjord bathymetry as a control. Note difference in axis scales between panels. Contouring of bed elevations above zero (red) is omitted in order to accentuate below-sea-level portions of glacier troughs. Same color code is used for both panels to accentuate differences in glacier beds. Panel (c) compares centerline bed depths for KNS (blue) vs NS (red). LIA moraines for the two glaciers are plotted coincidently. KNS began retreating from its LIA sill by 1761 AD; NS in 2004. Recent terminus positions are also shown.

3.6. Frontal ablation

Frontal ablation, Q _fa , at KNS was relatively constant between 2001 and 2013, averaging ~−4.4 ± 0.9 km³ a⁻¹, and then declined to ~−3.5 ± 0.6 km³ a⁻¹ in 2014 and 2015 (Table 4). Q _fa at NS increased from −1.5 ± 0.2 to −2.5 ± 0.3 km³ a⁻¹ between 2001 and 2006, and then nearly doubled to −4.2 ± 0.6 km³ a⁻¹ by 2013 (Table 4). Together, the total freshwater flux generated by Q _fa into Godthåbsfjord (Nuup Kangerlua) from KNS and NS increased from 6 to 8.5 km³ a⁻¹ (w.e.) between 2001 and 2013. Although our uncertainties are large, we note they are due in large part to systematic error associated with bed topography, which accounts for 50% of the total error. Thus, trends in fluxes are relatively well resolved.

3.7. HIRHAM5 meltwater runoff and SMB

There is considerable interannual variability in runoff at both KNS and NS (Fig. 10) with a significant peak in 2012 that coincides with a new melt record for the GrIS in July 2012. Other peaks occur in 2003, 2007 and 2010 at both KNS and NS. Total runoff into Godthåbsfjord from all glacier drainages averaged 16.9 km³ a⁻¹ during the period 1990–2004, increased to an average of 22.1 km³ a⁻¹ for 2004–12, and then declined to 13.8 km³ a⁻¹ in more recent years. At KNS, there is a significant correlation between the onset of spring runoff and the seasonal retreat of the KNS winter floating tongue and between the fall shutdown of runoff and the readvance of the tongue (Fig. 10a).

Fig. 10. Modeled monthly runoff predictions from HIRHAM5 are shown for (a) KNS and for (b) NS (inner left vertical axis). Terminus positions are indicated by black dots connected by dashed lines. In Panel (a), dash-dot line shows KNS summer retracted position; grey vertical bars show total annual runoff (outer left vertical axis). Vertical lines delineate correlations between onset and termination of runoff with advanced and retracted terminus positions.

Results for annual SMB (Fig. 11a) indicate that both KNS and NS had a relatively consistent positive balance regime that lasted from at least 1990 until ~2009 with SMB averaging 4.3 km³ a⁻¹ w.e. at KNS and ~1.7 km³ a⁻¹ w.e. at NS, although there is considerable interannual variability. SMB declined significantly at both glaciers during 2010–12, averaging 0.5 km³ a⁻¹ w.e. at KNS and −0.8 km³ a⁻¹ w.e. at NS before increasing to pre-2004 levels in 2013–15.

Fig. 11. (a) Annual SMB for individual drainages as determined from HIRHAM5. Red lines are results for tidewater glaciers: NS (filled circle), AS (open circle), and KNS (inverted triangle). Blue lines are land terminating glaciers: SS (filled circle), KS (open circle), QS (inverted triangle), KSS (open triangle) and IL (filled square). b) Comparison of cumulative SMB from HIRHAM5 (solid lines) with cumulative frontal ablation (dashed lines) (Table 4) at KNS (blue) and NS (red) for 2001–15. Note that Q _fa is shown as positive to facilitate comparison. Error bars for Q _fa were determined from propagation of uncertainties given in Table 4.

Table 4. Annual frontal ablation: Q _i , ice flux into the terminus_i; dV/dt, terminus volume change; Q _fa , frontal ablation (calving flux plus submarine melting); all in w.e. units

Average ice velocity into the terminus, Ū _i , in ice equivalent units.

^a Estimated using centerline velocities (Fig. 5).

Land-terminating glaciers all showed negative balance to varying degrees (Fig. 11a). These results qualitatively correspond with our geodetic assessment of ice loss for the land-terminating glaciers (Table 2). As at the tidewater glaciers, land-terminating glaciers also show a significant downturn in SMB during 2010–12. For SS, the HIRHAM5 SMB averaged −1.0 km³ a⁻¹ w.e. for the period 1990–2009, −5.7 km³ a⁻¹ w.e. during 2010–12, then −0.6 km³ a⁻¹ w.e. for 2013–15.

3.8. SMB vs frontal ablation

Here, we compare the balance between cumulative HIRHAM5 SMB and Q _fa at KNS and NS over our period of Q _fa data (2001–15) (Table 4; Fig. 11b). Cumulative SMB (C-SMB) and Q _fa (C-Q _fa ) were in balance until 2004, consistent with the observed relatively stable regimes at both KNS and NS (Fig. 2). The curves then diverge ~2005/06 with a stronger change occurring in 2011/12, with C-Q _fa significantly outpacing HIRHAM5 C-SMB. The changes correlate with retreats at both glaciers (Fig. 2), particularly for the second 2010 retreat at NS. The difference between C-Q _fa and C-SMB at NS reached 19 km³ in 2012 and continued to increase in following years. The difference at KNS was also greatest for 2012, and then narrowed slightly in succeeding years, when SMB reverted to a positive balance (Fig. 11a).

Frontal ablation at KNS did not fluctuate appreciably between 2001 and 2013, averaging ~4.4 ± 0.2 km³ a⁻¹ (w.e.) (Table 4). Based on these results and on the relative consistency of centerline velocity (Fig. 5c), we infer that this value is a reasonable estimate for KNS Q _fa extending back to 1990. Similarly, Q _fa averaged ~1.5 ± 0.2 km³ a⁻¹ w.e. between 2001 and 2004 at NS. Again, judging from the consistency of centerline velocity during the period 1996–2003 (Fig. 5d), we make the assumption that this value is reasonably representative of Q_fa at NS during 1990–2001. Based on these assumptions, there is a close balance between C-SMB and C-Q_fa at both glaciers between 1985 and 2001, which is consistent with the observed relative stability of both KNS and NS termini (Fig. 2).

3.9. Ice mélange variability

The time series of FST mélange conditions in Kangersuneq shows correlation between terminus retreat (Fig. 12a) and mélange variability (Fig. 12b). In general, retreat begins when FST rises above −5°C, and advances when FST drops below ~−5°C. The winter ice mélange appears as dark purple or black in the figure and frequently extends past Narsap Sermia's bay (white dashed line ~48 km). Winter ice mélange was quite pronounced into 2004, a period during which KNS and NS mean front positions did not significantly change. During the following winter, the mélange did not extend as far as NS bay or endure as long seasonally; NS initiated its retreat that fall and winter, while KNS reached a new minimum in summer 2005. The mélange was again spatially and temporally extensive during the winters of 2007, 2008, 2009 and neither glacier showed significant change. The mélange was at a minimum during the 2010 winter, and NS began its dramatic retreat that summer and fall. The mélange was persistent in front of KNS during the 2011, 2012, 2013 and 2014 winters, but did not extend as far downfjord as Narsap Sermia's bay; NS continued its large-scale retreat, while KNS showed some variations in summer terminus position.

Fig. 12. Panel (a) shows average terminus positions for KNS and NS. Panel (b) depicts results of FST data analysis derived from MODIS images. Panel (c) shows FST data (blue dots) for a position ~2 km from the KNS terminus (dashed cyan line in panel b). The black line in panel (c) is a 30 d running median value.

4. Submarine melting

The submarine melt rates for NS (Fig. 13) are given as an average over the area of the frontal glacier cross section in specific units (m d⁻¹); results for KNS (not shown) are slightly higher. Specific units allow direct comparison with near-terminus ice velocities and glacier length changes (Truffer and Motyka, 2016). Our estimate of peak Q _m during late summer and fall of 2008 and 2009 is ~3 m d⁻¹, which is similar to field measurements of glaciers north of Jakobshavn (Rignot and others, 2010), at Store Glacier (Xu and others, 2013), and at Kangerlussuup Sermia (Fried and others, 2015). Q _m at NS then increased to ~5–6 m d⁻¹ during summer-fall periods of 2010, 2011 and 2012, driven mainly by increased ocean temperature (Fig. 8) but also by increased Q _sg , particularly in 2012 (Fig. 10). Comparing Q _m with the average frontal ice velocity, Ū _i at NS (Table 4; Fig. 13) indicates that Q _m during the summers of 2008 and 2009 was about half the average terminus ice velocity, then increased and nearly equaled Ū _I in summers of 2010 and 2011.

Fig. 13. Estimates of submarine melting rates in m d⁻¹ for NS are in dark red using data from Figure 8. HIRHAM5 model runoff results were used as input for subglacial discharge (Q _sg). Parameters to obtain Q_m estimates were modified from Xu and others (2013) and Jenkins (2011). Black crosses and dashed line show Ū_i, average ice velocity (m d⁻¹) into the terminus (Table 4) for comparison. Grey bars depict periods of open mélange. Red line shows HIRHAM SMB trend (Fig. 10). Blue crosses and dashed line (L) show NS relative terminus position in km. Ellipse demarcates data associated with onset of major retreat.

5. DISCUSSION

5.1. Submarine melting

For submarine melting (Q _m) to plausibly impact terminus dynamics and initiate retreat, Q _m must, at least seasonally, be a significant fraction of frontal ablation, Q _fa . Submarine melting can also enhance and increase calving by simply undercutting the glacier face (Motyka and others, 2003; Bartholomaus and others, 2013; Fried and others, 2015), with one estimate suggesting an increase by a factor of 3–4 (O'Leary and Christoffersen, 2013). Indeed, the seasonal formation of embayments has been observed at many tidewater glaciers worldwide, invariably in close association with upwelling currents at the face of the glacier. The upwelling reflects buoyant convection of subglacial freshwater discharge from submarine conduits although the plume may not always reach the surface depending on the buoyancy of the entrained water (Jenkins, 2011; Sciascia and others, 2013). In Godthåbsfjord the plume has been observed to settle just below the summer surface layer (Mortensen and others, 2013; Bendtsen and others, 2015a, b) as well as upwelling at the face of KNS and NS. The rising plume entrains warm ambient seawater resulting in submarine melting that undercuts the terminus face, and causes the formation of deep concentric crevasses (Motyka and others, 2003, 2013; Mortensen and others, 2013; Bendtsen and others, 2015a). We have observed the formation of such features at KNS and NS during our field operations and from time-lapse camera photos. Landsat images indicate such embayments are common during summer at both NS and KNS. The development and expansion of such embayments are closely associated with retreats at both glaciers (Supplementary Material videos), indicating that submarine melting played a direct role in these retreats. For example, the increase in Q _m coincided directly with the inception of Narsap Sermia's major retreat that began in fall of 2010 (Fig. 13). At its peak during that summer, Q _m nearly equaled the average terminus ice velocity, U _i . Q _m increased during the next two summers, helping to sustain the continued retreat.

Looking at trends prior to 2008 at NS, runoff (or Q _sg ) gradually increased in the 3 years before the first retreat but then declined during the retreat of 2004–06 (Fig. 10). Thus, if an increase in Q _m was indeed involved in this first retreat it is more likely due to changes in ocean temperatures. The nearest record that overlaps the first retreat is ~200 km from NS, Fyllas Banke, on the continental shelf west of Godthåbsfjord (Fig. 17 in Ribergaard, 2014). Temperatures at depths of 50–150 m do show an increase of +1°C coinciding with the 2004–06 retreat. Warming of the West Greenland Current along southern and western Greenland coast, which flows past Fyllas Banke, has been implicated in triggering the retreat of Jakobshavn Isbrae (Holland and others, 2008; Motyka and others, 2011).

Another record from a station west of Nuuk and just inside Godthåbsfjord but still at considerable distance (150 km) from NS, spans the summer of 2005 to winter 2010 (Mortensen and others, 2011, their Fig. 5). The temperatures show seasonal oscillations at a depth of 240 m ranging from a winter low of −1 to ~4°C in summers. Most importantly, summer temperatures reached 4°C in 2005 then dropped in succeeding years, especially in 2009, before they warmed again to 4°C by 2010. However, flow of such waters into Godthåbsfjord is complex (Mortensen and others, 2011) because three prominent sills and tidal mixing modulate flow of ocean/fjord waters into the inner part of Kangersuneq (Mortensen and others, 2011, 2013, 2014), making any conclusions that increased submarine melting influenced the 2004 retreat equivocal.

At KNS, Q _m (≈3–7 m d⁻¹) is a comparatively lower fraction of Ū _I (≈14 m d⁻¹) than at NS, but still significant. Looking at summer positions, as with NS, an increase in ocean temperatures at the mouth of Godthåbsfjord coincides with the KNS 2004–06 retreat. Cooler inner fjord waters in 2008/09 coincide with a slight advance of KNS while warmer waters during 2010–12, plus a significant increase in runoff at KNS during 2010–12, coincide with a retreat.

5.2. Ice mélange, runoff and tidewater glacier stability

Velocity fields that are smooth and continuous across calving fronts and into ice mélanges are evidence of mélange rigidity (e.g. Joughin and others, 2008; Amundson and others, 2010; Foga and others, 2014; Cassotto and others, 2015). We see similar effects in our velocity data for KNS and NS. While back stress imposed by a rigid mélange is insufficient to directly influence glacier motion (e.g. Cook and others, 2014), only a small back stress is required to prevent detached icebergs from overturning (Amundson and others, 2010). The back stress can then indirectly affect glacier dynamics by inhibiting calving and permitting a floating tongue to develop (Cassotto and others, 2015). Landsat images (cf. Supplementary Material Video) show such a floating tongue seasonally developing at KNS.

The correlations of ice mélange extent and the meltwater runoff with terminus positions of KNS and NS termini on annual and interannual timescales indicate an interrelationship between the three. Seasonal variations in runoff directly impact the magnitude of subglacial discharge and thus Q _m. We hypothesize that such changes directly affect the floating tongue by progressively increasing submarine melting under the tongue as runoff increases through spring and summer (Fig. 13). Once runoff shuts down, Q _m becomes much reduced, Q _fa declines and the terminus advects a floating tongue once again. Mélange variability and runoff processes are thus likely interrelated and combine to produce the observed change in long-term stability at both glaciers as well as the seasonal cycle at KNS. Warm air temperatures trigger surface melt, which initiates runoff. They also melt sea ice, which weakens the fabric that binds iceberg clasts within the mélange.

5.3. What triggered the retreat of NS?

A steady position of a tidewater glacier's terminus requires a balance between ice supply to the front and ice loss by frontal ablation (Eqn (1) with dV/dt = 0). Several factors can initiate the rapid retreat of a tidewater glacier. Thinning of a glacier terminus by a reduction in SMB or other processes results in increased buoyancy, stimulating increases in mechanical calving, increasing Q _fa and destabilizing the terminus. Changes in Q _fa can also be initiated by changes in Q _m, via ocean temperature and/or Q _sg , and consequent enhanced calving and embayments, with accelerated ice flow quickly thinning the terminus region and rendering it buoyant. Lastly, changes in the proglacial mélange with feedbacks on calving can also affect Q _fa .

The ellipse in Figure 13 encircles various processes coincident with the inception of the 2010 retreat, a period that had the largest duration of relatively open water, a significant increase in estimated Q _m, a decline in SMB and a reduction in Ū _i . Not shown in Figure 13, however, is the 40 m thinning of NS (Fig. 4a), which likely occurred between 2004 and 2008. The thinning is a dynamic effect resulting from the acceleration of NS following the 2004 retreat (Figs 5b, d). This thinning made the terminus increasingly buoyant, particularly if the 2004 retreat caused the terminus to retract onto a reverse slope, thus increasing its vulnerability to further retreat. The glacier stabilized after 2008 and thinning and retreat did not resume again until 2010/11. The post-2010 retreat drove the terminus into even deeper waters (Fig. 9), increasing buoyancy and dynamic feedbacks. We conclude that the combination of the preceding processes initiated the retreat in late 2010. However, the evidence indicates that submarine melting was the most prominent of these processes with significant increases in Q _m coincident with the evolution of embayments and glacier retreat.

Deciphering what initiated the 2004 retreat is more problematic. Simulations of SMB using HIRHAM5 do not show any significant changes in SMB between 1960 and 2004 (Langen and others, 2015). The duration of open mélange was higher during 2004–06, which perhaps resulted in NS not fully advancing to its anchoring moraine during winter months. There are also indications of increased ocean temperatures near the mouth of Godthåbsfjord suggesting submarine melt may have had an influence, particularly in view of the large embayment that developed before this retreat. However, what triggered the initial 2004 retreat remains ambiguous because we lack sufficient ice-proximal ocean data for this time period.

5.4. Asynchronous behavior

What causes one tidewater glacier to retreat (or advance) and not a neighboring one depends on the size of the glacier and its geometry (e.g., bed depths and pinning points) and its past glacial history (Mercer, 1961; Meier and Post, 1987; Post and Motyka, 1995; Post and others, 2011; Truffer and Motyka, 2016), although other factors may also come into play. NS and KNS were both anchored on 160 m deep LIA fjord moraines and had similar outlet widths (~4.5 km) when each started their respective retreats into 300+ m over-deepened basins. However, KNS began retreating 240 years before NS began its retreat. Although the termini of the two glaciers are only ~40 km apart, a difference in climate regime may still be a factor. However, we lack data on atmospheric conditions in either region as the PROMICE data series is too short to make any conclusions. Assessments of DEMs (Table 2) and SMB comparisons (Fig. 11) show that land-terminating glaciers in the southern part of the study area all exhibited strong ice losses between 1985 and 2008. In contrast, SS and KS, on either side of NS, did not begin to show much change until after 2008, about the same time that NS began retreating. In the future the continued drawdown of NS may lead to divide migration and piracy of ice flow from SS and KS.

One significant difference between KNS and NS is their respective drainage areas: KNS is 2.4 times more extensive than NS (Table 2). But why this should have affected the timing of the two retreats is not clear. Another complicating factor is ocean environment. One might presume that given their proximity both glaciers would have been affected by the same ocean conditions. However, KNS lies directly in line with the axis of the principal fjord, whereas NS sits in a tributary valley. It could be that water and sea ice properties in the proximity of KNS and NS differ, driven by differences in external forces such as Coriolis effects, wind regimes and subglacial discharge.

Future changes at KNS and NS are intertwined with their respective channel geometries as well as with any changes in oceanic and atmospheric conditions. Figures 9b, c indicate that although the terminus of NS is currently located in an over-deepening, an additional retreat of a kilometer or so would bring the terminus up onto a ~150 m sill, potentially reducing frontal ablation and buoyancy, and resulting in a standstill. Should NS eventually retreat beyond this sill, the glacier will find itself in a deep trough (up to −1000 m). Buoyancy and glacier dynamics will then dominate the system, with feedbacks enhancing instability (e.g., Pfeffer, 2007; Vieli and Nick, 2011), potentially leading to a ~30 km retreat. In contrast, a retreat of KNS would put the terminus into progressively shallower bed topography (Figs 9a, c), which would likely stabilize the glacier for the near future. However, these predictions must be treated with caution given the large uncertainties associated with airborne radio echo sounding of glacier valleys.

6. CONCLUSIONS

We used multiple datasets to document changes to glaciers draining into inner Godthåbsfjord. Regional volume loss was 29.13 ± 0.13 km³ between 1985 and 2008 with KNS accounting for nearly half of this loss. An additional 14.67 ± 0.12 km³ was lost between 2008 and 2014, with 74% of the loss from NS. The total ice loss is equivalent to 0.10 mm eustatic sea-level rise. Including terminus retreat of the tidewater outlet glaciers (3.5 ± 0.3 km³) total ice loss between 1985 and 2014 was 47.3 ± 0.4 km³.

Tidewater glaciers KNS and NS have exhibited contrasting behavior. KNS has retreated 22 km from its 160 m deep LIA anchoring moraine while the terminus thinned 550 m since the LIA. The retreat began during the late 18th century with a brief standstill ~1920 (Lea and others, 2014a). KNS experienced considerable drawdown (100 m) and ice loss between 1985 and 2008 with the majority occurring between 2001 and 2008. Retreat of the KNS summer minimum positions has been relatively minor over the period 1985–2015: ~0.8 km. The glacier has retreated into a narrower and shallower part of its channel, constricting access to the calving front. Glacier velocities near the KNS terminus have remained relatively constant since 2001, ranging from 5 to 6 km a⁻¹. KNS appears to have stabilized, oscillating between small advances and retreats, and between thickening and thinning.

The formation and breakup of KNS's seasonal floating tongue appears to be caused by a combination of processes. In spring, increasing atmospheric temperatures initiate runoff, which increases the rate of submarine melting under the floating tongue, thus thinning and weakening it. The proglacial ice mélange is also weakened by the higher temperatures, and perhaps by submarine melting, reducing back-pressure on the terminus face. These processes eventually result in the collapse and disintegration of the seasonal tongue. In fall, the atmosphere cools below freezing, runoff shuts down, drastically reducing the rate of submarine melting and sea ice begins forming, strengthening the ice mélange, conditions, which lead to the advance of a floating tongue.

In contrast to KNS, NS did not advance a seasonal floating tongue and remained stable at its LIA maximum until 2004 when an initial retreat (~0.6 km) led to a doubling in speed from 1.5 to 3.0 km a⁻¹. NS then restabilized before undergoing a major calving retreat of 3.3 km between 2010 and 2013. The glacier now terminates in ~300 m deep water and is retreating at 150 m a⁻¹ (2013–15) with terminus velocities of ~5.5 km a⁻¹. An increase in submarine melting and variability in the strength of seasonal ice mélange may have triggered the initial retreat in 2004. A combination of processes likely triggered the 2010 retreat including dynamic thinning, longer duration of open mélange, a decline in SMB, but mainly an increase in submarine melting due to increased sea temperatures and increased subglacial discharge. Ice dynamics are likely now the dominant control as the glacier retreats into deeper water, the terminus thins and its buoyancy increases. NS may continue to retreat for another 1 km before stabilizing on a 150 m deep sill. However, if retreat continues beyond this sill, the glacier would likely destabilize and begin a 30-km long catastrophic calving retreat.