1. INTRODUCTION
Glaciers along the margin of the Greenland ice sheet (GrIS) have been thinning and retreating over the past several decades (Warren, Reference Warren1991; Warren and Glasser, Reference Warren and Glasser1992; Moon and Joughin, Reference Moon and Joughin2008; Leclercq and others, Reference Leclercq2012) with tidewater outlet glaciers accounting for much of the ice loss: up to 58% before 2005 (e.g. Rignot and Kanagaratnam, Reference Rignot and Kanagaratnam2006; Van den Broeke and others, Reference Van den Broeke2009; Rignot and others, Reference Rignot, Velicogna, van den Broeke, Monaghan and Lenaerts2011; Shepherd and others, Reference Shepherd2012) and 32% between 2009 and 2012 (Enderlin and others, Reference Enderlin2014). 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 (Reference van As2014). 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, Reference Rignot and Kanagaratnam2006; Joughin and others, Reference Joughin, Smith, Howat, Scambos and Moon2010a). 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, Reference Weidick, Bennike, Citterio and Nørgaard-Pedersen2012; Lea and others, Reference Lea2014a). 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, Reference Weidick, Bennike, Citterio and Nørgaard-Pedersen2012; Lea and others, Reference Lea2014a). 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, Reference Post, O'Neel, Motyka and Streveler2011; McNabb and Hock, Reference McNabb and Hock2014; Truffer and Motyka, Reference Truffer and Motyka2016) as well as for other GrIS outlet glaciers (e.g., Carr and others, Reference Carr, Stokes and Vieli2013; Bartholomaus and others, Reference Bartholomaus2016).
A number of studies have addressed changing conditions that could impact Godthåbsfjord glaciers. Several have examined ocean conditions (Mortensen and others, Reference Mortensen, Lennert, Bendtsen and Rysgaard2011, Reference Mortensen2013, Reference Mortensen, Bendtsen, Lennert and Rysgaard2014; Bendtsen and others, Reference Bendtsen, Mortensen, Lennert and Rysgaard2015a) while van As and others (Reference van As2014) and Langen and others (Reference Langen2015) used regional climate models to estimate freshwater runoff into the fjord and surface mass balance (SMB). Weidick and others (Reference Weidick, Bennike, Citterio and Nørgaard-Pedersen2012) and Lea and others (Reference Lea2014a, Reference Leab) discussed the KNS LIA boundary and its post LIA retreat while Lea and others (Reference Lea2014b) 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, Reference Walter2010) and modeling (e.g., Pfeffer, Reference Pfeffer2007; Vieli and Nick, Reference Vieli and Nick2011); 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, Reference Straneo2013; Moon and others, Reference Moon, Joughin and Smith2015; Truffer and Motyka, Reference Truffer and Motyka2016). 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, Reference Langen2015). 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 km2 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, Reference Korona, Berthier, Bernard, Remy and Thouvenot2009). 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, Reference Noh and Howat2015) 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 (Reference Kjær2012), Kjeldsen and others (Reference Kjeldsen2015) and Korsgaard and others (Reference Korsgaard, Nuth, Khan, Kjeldsen, Bjørk, Schomacker and Kjær2016), based on DEM co-registration to ICESat data for this region.
Korona and others (Reference Korona, Berthier, Bernard, Remy and Thouvenot2009) 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 (Reference Noh and Howat2015) 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, Reference Hallet, Hunter and Bogen1996).
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 (Reference Rolstad, Haug and Denby2009) and outlined in Motyka and others (Reference Motyka, Fahnestock and Truffer2010). 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 km2 for both, which is considerably smaller than the area, A, for both the entire region of coverage (~2100 km2), 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−3.
2.2. Airborne laser altimetry
2.2.1. Airborne topographic mapper
Table 1 lists NASA ATM data for KNS and NS (Krabill, Reference Krabill2014). 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 (Reference Cassotto, Fahnestock, Amundson, Truffer and Joughin2015). 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, Reference Joughin, Smith, Howat and Scambos2010b). 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 (Reference Joughin, Smith, Howat, Scambos and Moon2010a).
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, Reference Fahnestock2015); 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 (Reference Morlighem, Rignot, Mouginot, Seroussi and Larour2014) 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, Reference Morlighem, Rignot, Mouginot, Seroussi and Larour2015). 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, Reference Morlighem2011). 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):
$$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, Reference O'Neel, Echelmeyer and Motyka2003).
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
$$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−1 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:
$$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
$$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:
$${\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, Reference Morlighem, Rignot, Mouginot, Seroussi and Larour2014). 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−1) 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 (Reference Morlighem, Rignot, Mouginot, Seroussi and Larour2015) 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 (Reference van As2014) and were determined using a Geographic Information Systems particle tracking tool and the 2006 velocity field (Joughin and others, Reference Joughin, Smith, Howat, Scambos and Moon2010a). Details of applying the HIRHAM5 model to glaciers in the Godthåbsfjord region are discussed in Langen and others (Reference Langen2015). The HIRHAM5 for Greenland has a resolution of ~5.5 km. Analyses by Ettema and others (Reference Ettema2009) and Lucas-Picher and others (Reference Lucas-Picher2012) 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 (Reference Bamber, Ekholm and Krabill2001), 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 (Reference Langen2015) 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, Reference Lucas-Picher2012).
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 (Reference Mortensen2013, 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 (Reference Mortensen, Lennert, Bendtsen and Rysgaard2011, Reference Mortensen2013) and Ribergaard (Reference Ribergaard2014).
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, Reference Peters2015). 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 (Reference Cassotto, Fahnestock, Amundson, Truffer and Joughin2015). 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, Reference Feldman and McClain2015). 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, Reference Motyka, Hunter, Echelmeyer and Connor2003, Reference Motyka, Dryer, Amundson, Truffer and Fahnestock2013) and modeling experiments (e.g., Jenkins, Reference Jenkins2011) 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, Reference Jenkins2011; Xu and others, Reference Xu, Rignot, Fenty, Menemenlis and Flexas2013) 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, Reference Lea2014a).
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−1. 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, Reference Amundson2008; Joughin and others, Reference Joughin2008; Cassotto and others, Reference Cassotto, Fahnestock, Amundson, Truffer and Joughin2015). 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, Reference Weidick, Bennike, Citterio and Nørgaard-Pedersen2012) 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 km3 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 km3 of ice loss for a total ice loss between 1985 and 2014 of 47.3 ± 0.4 km3. 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 km3 vs 2.7 km3; 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−1 in 2000/01 to a high of ~7 km a−1 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−1 between 1987 and 2011, then increased to ~6 km a−1 by 2012 but then declined back to 4–5 km a−1 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−1, and then nearly doubling again by February 2013 (Fig. 5b). Speeds accelerated all along the glacier and ranged from 5.5 km a−1 at the terminus to 1.5 km a−1 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−1 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 (Reference Mortensen, Lennert, Bendtsen and Rysgaard2011, Reference Mortensen2013), 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.
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 (Reference Mortensen2013) 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 km3 a−1, and then declined to ~−3.5 ± 0.6 km3 a−1 in 2014 and 2015 (Table 4). Q fa at NS increased from −1.5 ± 0.2 to −2.5 ± 0.3 km3 a−1 between 2001 and 2006, and then nearly doubled to −4.2 ± 0.6 km3 a−1 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 km3 a−1 (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 km3 a−1 during the period 1990–2004, increased to an average of 22.1 km3 a−1 for 2004–12, and then declined to 13.8 km3 a−1 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 km3 a−1 w.e. at KNS and ~1.7 km3 a−1 w.e. at NS, although there is considerable interannual variability. SMB declined significantly at both glaciers during 2010–12, averaging 0.5 km3 a−1 w.e. at KNS and −0.8 km3 a−1 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 terminusi; 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 km3 a−1 w.e. for the period 1990–2009, −5.7 km3 a−1 w.e. during 2010–12, then −0.6 km3 a−1 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 km3 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 km3 a−1 (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 km3 a−1 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 Qfa at NS during 1990–2001. Based on these assumptions, there is a close balance between C-SMB and C-Qfa 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−1); 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, Reference Truffer and Motyka2016). Our estimate of peak Q m during late summer and fall of 2008 and 2009 is ~3 m d−1, which is similar to field measurements of glaciers north of Jakobshavn (Rignot and others, Reference Rignot, Koppes and Velicogna2010), at Store Glacier (Xu and others, Reference Xu, Rignot, Fenty, Menemenlis and Flexas2013), and at Kangerlussuup Sermia (Fried and others, Reference Fried2015). Q m at NS then increased to ~5–6 m d−1 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−1 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 Qm estimates were modified from Xu and others (Reference Xu, Rignot, Fenty, Menemenlis and Flexas2013) and Jenkins (Reference Jenkins2011). Black crosses and dashed line show Ūi, average ice velocity (m d−1) 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, Reference Motyka, Hunter, Echelmeyer and Connor2003; Bartholomaus and others, Reference Bartholomaus, Larsen and O'Neel2013; Fried and others, Reference Fried2015), with one estimate suggesting an increase by a factor of 3–4 (O'Leary and Christoffersen, Reference O'Leary and Christoffersen2013). 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, Reference Jenkins2011; Sciascia and others, Reference Sciascia, Straneo, Cenedese and Heimbach2013). In Godthåbsfjord the plume has been observed to settle just below the summer surface layer (Mortensen and others, Reference Mortensen2013; Bendtsen and others, Reference Bendtsen, Mortensen, Lennert and Rysgaard2015a, Reference Bendtsen, Mortensen and Rysgaardb) 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, Reference Motyka, Hunter, Echelmeyer and Connor2003, Reference Motyka, Dryer, Amundson, Truffer and Fahnestock2013; Mortensen and others, Reference Mortensen2013; Bendtsen and others, Reference Bendtsen, Mortensen, Lennert and Rysgaard2015a). 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, Reference Ribergaard2014). 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, Reference Holland, Thomas, de Young, Ribergaard and Lyberth2008; Motyka and others, Reference Motyka2011).
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, Reference Mortensen, Lennert, Bendtsen and Rysgaard2011, 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, Reference Mortensen, Lennert, Bendtsen and Rysgaard2011) because three prominent sills and tidal mixing modulate flow of ocean/fjord waters into the inner part of Kangersuneq (Mortensen and others, Reference Mortensen, Lennert, Bendtsen and Rysgaard2011, Reference Mortensen2013, Reference Mortensen, Bendtsen, Lennert and Rysgaard2014), making any conclusions that increased submarine melting influenced the 2004 retreat equivocal.
At KNS, Q m (≈3–7 m d−1) is a comparatively lower fraction of Ū I (≈14 m d−1) 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, Reference Joughin2008; Amundson and others, Reference Amundson2010; Foga and others, Reference Foga, Stearns and Van der Veen2014; Cassotto and others, Reference Cassotto, Fahnestock, Amundson, Truffer and Joughin2015). 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, Reference Cook2014), only a small back stress is required to prevent detached icebergs from overturning (Amundson and others, Reference Amundson2010). The back stress can then indirectly affect glacier dynamics by inhibiting calving and permitting a floating tongue to develop (Cassotto and others, Reference Cassotto, Fahnestock, Amundson, Truffer and Joughin2015). 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, Reference Langen2015). 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, Reference Mercer1961; Meier and Post, Reference Meier and Post1987; Post and Motyka, Reference Post and Motyka1995; Post and others, Reference Post, O'Neel, Motyka and Streveler2011; Truffer and Motyka, Reference Truffer and Motyka2016), 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, Reference Pfeffer2007; Vieli and Nick, Reference Vieli and Nick2011), 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 km3 between 1985 and 2008 with KNS accounting for nearly half of this loss. An additional 14.67 ± 0.12 km3 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 km3) total ice loss between 1985 and 2014 was 47.3 ± 0.4 km3.
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, Reference Lea2014a). 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−1. 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−1. 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−1 (2013–15) with terminus velocities of ~5.5 km a−1. 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.
SUPPLEMENTARY MATERIAL
The supplementary material for this article can be found at http://dx.doi.org/10.1017/jog.2016.138
AUTHOR CONTRIBUTION STATEMENT
R. Motyka wrote most of the paper with the help of R. Cassotto and all coauthors. R. Motyka performed DEM differencing, analyses of ATM data and submarine melting, R. Cassotto derived terminus positions and analyzed FST and ice mélanges, M. Truffer performed ice flux calculations, N. J. Korsgaard developed the 1985 DEM, D. van As provided input on modeling SMB and runoff, P. Langen provided the HIRHAM5 results, M. Fahnestock analyzed Landsat images for velocity data, J. Mortensen and S. Rysgaard provided the hydrographic time series for Kangersuneq, I. Howat developed the 2014 DEM and K. Lennert provided the helicopter depth soundings.
ACKNOWLEDGEMENTS
We thank W. Dryer and D. Podrasky for assistance with fieldwork and L. Kenefic for assisting with digitizing glacier front positions. CH2 M HILL Polar Services and Air Greenland provided logistics support. The SPOT-5 images used for the 2008 DEM were provided by the SPIRIT program (Centre National d'Etudes Spatiales, France). The DigitalGlobe Worldview images used for the 2014 DEM were obtained from P. Morin. Terminus positions were derived from Landsat images courtesy of the U.S. Geological Survey. Funding was provided by the US National Science Foundation (NSF) Office of Polar Programs (OPP) grants NSF PLR-0909552 and NSF PLR-0909333. Cassotto is supported by NASA under the Earth and Space Science Fellowship Program (Grant NNX14AL29H). K. K. Kjeldsen acknowledges support from the Danish Council Research for Independent Research (grant no. DFF - 4090-00151). K. Kjær is thanked for his support during the earlier phases of this study. On-ice weather stations are operated by GEUS (Denmark) within the Programme for Monitoring of the Greenland Ice Sheet (PROMICE). J. Mortensen acknowledges support from IIKNN (Greenland), DEFROST project of the Nordic Centre of Excellence program “Interaction between Climate Change and the Cryosphere” and the Greenland Ecosystem Monitoring Programme (www.g-e-m.dk). S. Rysgaard was funded by the Canada Excellence Research Chair Programme. Additional funding was provided by the Geophysical Institute, University of Alaska Fairbanks, and Greenland Climate Research Centre. Scientific editor H. Fricker and reviewers H. Jiskoot and G. Cogley provided very constructive feedback that helped improve the paper.
Open access
