Research ArticleGEOLOGY

Atmospheric forcing of rapid marine-terminating glacier retreat in the Canadian Arctic Archipelago

See allHide authors and affiliations

Science Advances  13 Mar 2019:
Vol. 5, no. 3, eaau8507
DOI: 10.1126/sciadv.aau8507


The Canadian Arctic Archipelago contains >300 glaciers that terminate in the ocean, but little is known about changes in their frontal positions in response to recent changes in the ocean-climate system. Here, we examine changes in glacier frontal positions since the 1950s and investigate the relative influence of oceanic temperature versus atmospheric temperature. Over 94% of glaciers retreated between 1958 and 2015, with a region-wide trend of gradual retreat before ~2000, followed by a fivefold increase in retreat rates up to 2015. Retreat patterns show no correlation with changes in subsurface ocean temperatures, in clear contrast to the dominance of ocean forcing in western Greenland and elsewhere. Rather, significant correlations with surface melt indicate that increased atmospheric temperature has been the primary driver of the acceleration in marine-terminating glacier frontal retreat in this region.


In the northern high latitudes, one of the regions thought to be highly vulnerable to environmental change is the Canadian Arctic Archipelago (CAA), which contains the greatest area of global land ice outside the ice sheets in Greenland and Antarctica (1). With a glaciated area of ~146,000 km2, the CAA has the potential to contribute 199 ± 30 mm to global mean sea level rise (2). Rates of ice mass loss from the CAA have increased sharply during the 21st century (36), from the Queen Elizabeth Islands (QEI) in the north [31 ± 8 gigatons (Gt) year−1 in 2004–2006, increasing to 92 ± 12 Gt year−1 in 2007–2009] (3) to Baffin and Bylot Islands (BBI) further south (from 11.1 ± 3.4 Gt year−1 in 1963–2006 to 23.8 ± 6.12 Gt year−1 in 2003–2011) (7). Although total ice mass loss is known, changes in individual glaciers have not been analyzed regionally.

There is increasing evidence in both polar regions that rates of frontal ablation of ice flowing into the ocean are strongly controlled by subsurface ocean temperatures (810). Unprecedented warming of the Arctic Ocean and surrounding water bodies (11, 12) therefore makes pan-Arctic marine-terminating glaciers vulnerable to rapid retreat (13, 14). Oceanic forcing has already been recognized as a key control on the fluctuations of outlet glaciers in western Greenland (1517), where rates of submarine melting are up to two orders of magnitude greater than surface melt rates (18). Warming of subpolar North Atlantic water entering Baffin Bay is considered to be triggering an increase in the rate of submarine melting of Greenland’s outlet glaciers (14). However, on the Canadian side of Baffin Bay, the role of oceanic forcing on marine-terminating glaciers has not been directly investigated and is largely unknown, although it has been hypothesized that ocean temperatures might be influencing glacier retreat (19).

With over 300 marine-terminating glaciers in the CAA, dynamic feedbacks associated with glacier frontal changes have clear potential to influence decadal-scale ice losses from the region. Recent studies have shown that increasing ice discharge to the ocean has contributed to recent mass loss from the CAA in some locations (6, 20). For example, ~62% of current ice discharge (i.e., via calving) from the QEI is channeled through Trinity and Wykeham glaciers, southeast Ellesmere Island, where flow speed has markedly accelerated over the past 10 to 15 years (20). In this area, at least half of the observed thinning rate at the glacier termini has been attributed to changes in ice dynamics rather than to surface mass balance (SMB) changes (20). Changes taking place at the terminus of marine-terminating glaciers form a crucial component of mass balance (21) and have a major influence on coastal ecosystems, and yet, long-term regional patterns and trends in the frontal changes of glaciers in the CAA have never been examined. Here, we present a comprehensive analysis of decadal changes in the frontal positions of 300 marine-terminating glaciers in the region. We analyze ocean temperature measurements, together with downscaled 1-km-resolution regional climate model output, to evaluate the relative importance of oceanographic and atmospheric influences on the observed changes in glacier frontal positions since the late 1950s.


Marine-terminating glacier frontal changes

We determined ~7 frontal positions for each glacier by digitizing images from a range of sources collected at approximately decadal intervals (see Materials and Methods, fig. S1, and table S1). Our comparison of historical aerial photography from 1958/1959 with Landsat satellite imagery in 2015 indicates that >94% of the 300 marine-terminating glaciers investigated have, overall, retreated during this period (Fig. 1). The mean frontal change rate was −9.3 m year−1 (±1.4 m year−1 SE) for the 211 glaciers in the QEI and −6.9 m year−1 (±0.7 m year−1 SE) for the 89 glaciers in BBI. Total area loss by frontal retreat since 1958 was 387 km2 in the QEI and 22.5 km2 in BBI. In both regions, the absolute area change associated with glacier retreat is related to both the width and length of the glacier (in absolute terms, larger glaciers retreat greater distances than smaller ones), so we normalized area changes according to the length of glacier basins in ~2000. Mean relative length changes since the late 1950s were −4.2% (±0.3% SE) in the QEI and −4.1% (±0.3% SE) in BBI (Fig. 1). Of the 211 marine-terminating glaciers in the QEI, 23 have previously been identified as surge type (22), so we excluded these from further analysis because their terminus behavior can be independent of climate forcing. In addition, our study does not include glaciers situated along the northern coast of Ellesmere Island, which are strongly influenced by the buttressing effect of multiyear sea ice and interactions with floating ice shelves and have been described elsewhere (23).

Fig. 1 Spatial patterns of marine-terminating glacier front change from 1958/1959 to 2015 (as the percentage of glacier length in 2000).

In the QEI, 197 glaciers have retreated, and 12 have advanced since 1959. Five of those that have advanced are surge type. In the BBI region, 87 glaciers have retreated, and 2 have advanced since 1958. Named features are those referred to in text. GF, Grise Fiord weather station; PoW Icefield, Prince of Wales Icefield; Manson IF, Manson Icefield; T-W glaciers, Trinity and Wykeham glaciers.

Similar decadal patterns in frontal variability occurred throughout the study region: Glaciers showed slow retreat in the earlier years (1958 to late 1990s), and retreat rates increased considerably after ~2000 (Fig. 2, A and B). The mean frontal change rates in the QEI and BBI during the 1960s were −7.0 and −0.6 m year−1, respectively, increasing to −29.1 and −22.9 m year−1 since 2009/2010 (table S2). This trend is evident in both absolute and relative length changes and across all glacier types (e.g., different sizes, lengths, and frontal widths). The greatest absolute retreat rates occurred in the QEI, where mean glacier lengths and terminus widths are greater than those of glaciers further south. One notable anomaly in the BBI region is a period of increased retreat in the early 1980s (Fig. 2B). Another period of note is the earliest interval in the QEI region (1959–1975), when retreat rates were greater than in the two subsequent periods (Fig. 2A).

Fig. 2 Temporal trends in glacier length, surface ablation and runoff.

Marine-terminating glacier length changes (percentage per year) by time interval for (A) QEI and (B) BBI. Trends in mean runoff (millimeter w.e. year–1) and mean ablation area (as the percentage of glacier basin area) calculated from RACMO2.3 statistically downscaled to 1 km for (C) QEI and (D) BBI. Surge-type glaciers are excluded. Thick horizontal lines are mean values, and outer horizontal lines are ±1 SE.

Similar trends are observed across different latitudes, showing that there are no clear spatial variations in the temporal pattern of accelerating retreat (fig. S2). At all latitudes in the QEI, the lowest retreat rates occurred in the earlier periods, and rates were highest in the two periods since 2002. Although mean relative retreat rates are lower in the two most northerly latitudinal bands, they exhibit the same temporal trend. In the BBI region, the pattern is broadly similar, but during the period 1980–1985, there was a notable increase in retreat in each band north of 70°N. The highest regional mean retreat rate (−0.33% year−1) occurred in the most recent period, 2009–2015, in southeast Ellesmere Island (77°N to 78°N). The greatest absolute retreat rates in the CAA were observed on the coalescent Trinity and Wykeham glaciers (78.0°N, 78.6°W), where they reached −338 m year−1 between 1976 and 1984 and − 255 m year−1 between 2009 and 2015, with a total retreat of 4.8 km since 1959.

Ocean temperatures and regional frontal change

The coherency in glacier retreat rates across the region suggests a common external driver. Subsurface ocean temperatures on the eastern side of Baffin Bay have been identified as a significant control on frontal change rates of glaciers in western Greenland (1416, 18, 24), and so, we first consider the possible influence of ocean temperatures on glacier fluctuations on the Canadian side of Baffin Bay. Although ocean temperature measurements from locations close to glacier fronts in the CAA are generally sparse, there are sufficient data sources for regional-scale analyses. TOPAZ4 reanalysis data (gridded at 12.5 km from 1991 to 2015) (25) are available from the Copernicus Marine Environment Monitoring Service (CMEMS) and show spatial patterns in mean temperatures at standard depths (Fig. 3). Strong differences in temperature between the east and west sides of Baffin Bay are evident: Off southwest Greenland, warm waters (>1°C) are present at all depths, whereas waters at depths of 200 m and shallower surrounding the CAA are typically cold (<−1°C). This pattern is explained by the cyclonic circulation in northern Baffin Bay from the West Greenland Current System (26), which advects warm, saline Atlantic waters northward along western Greenland and carries fresh and cold Arctic waters southward along eastern Baffin Island. Additional cold water masses occur off northern Baffin Island as a result of wintertime freezing within the polynyas of northern Baffin Bay (in Smith, Jones, and Lancaster Sounds; locations marked on Fig. 1) (27). Cold waters in Smith Sound originate from water from the Pacific that enters via the Arctic Ocean and Nares Strait (26). In recent decades, however, changes have occurred in Baffin Bay, including warming in areas affected by the Atlantic inflow (12).

Fig. 3 Mean annual ocean temperatures at standard depths of 50, 100, 200, and 400 m.

Data are from CMEMS TOPAZ4 Arctic Ocean Physics reanalysis at 12.5-km grid spacing, modeled from all available satellite and in situ observations between 1991 and 2015. The black line is the 100-km buffer offshore from the CAA coastline.

To assess ocean temperatures close to marine-terminating glaciers in the CAA, we use in situ temperature measurements (1972–2015) obtained from the World Ocean Database 2013 (WOD13) (28) for standard depths between 50 and 500 m. These depths represent the range of bathymetry close to the coast, as calculated from the National Oceanic and Atmospheric Administration (NOAA) International Bathymetric Chart of the Arctic Ocean (see Materials and Methods). Although only 19 glaciers have ocean temperature measurements within 5 km of their fronts, these show comparable mean temperatures to those with measurements up to 20 km away (figs. S3 and S4). In total, 117 glaciers have ≥10 conductivity-temperature-depth (CTD) measurements within 20 km of their fronts, which provides a sample that is more than sufficient for statistical analysis. First, following methods outlined in (29), we compare mean overall glacier frontal changes (1958/1959 to 2015) with their adjacent mean ocean temperatures (1972 to 2015). From this, we find no statistically significant correlation at any depth (Fig. 4 and table S4). The mean ocean temperatures within 20 km of glacier fronts are below 0°C at all depths (fig. S4), and although these are above the freezing temperature of sea water (calculated from pressure and mean salinity), the melt induced by water at such low mean temperatures would be insufficient to drive retreat at the scale observed (29). For example, a study of glaciers in Svalbard (10) found that the relationship between frontal ablation rate and ambient subsurface ocean temperatures was strong and linear above 0°C, but little or no retreat occurred when ocean temperatures were below 0°C.

Fig. 4 Mean overall (1958–2015) glacier front changes (as the percentage of glacier length in 2000) by mean in situ ocean temperatures (1971–2015) at specific depths.

Results are for all glaciers with 10 or more ocean temperature measurements within 20 km of their fronts (glacier counts are in table S3, and correlation coefficients are in table S4). No significant correlation is observed between magnitude of frontal retreat and temperatures at any depth.

Spatial and temporal analysis of in situ temperature measurements shows that, north of 75°N, mean ocean temperatures up to 100 km from the coast have been below 0°C almost universally since the earliest records (Fig. 5 and table S5). Although heat content is low at these northern latitudes, a statistically significant increase in temperature occurred between the 1970s and 2000s at all depths (where data are available). Another study has shown that, on a local scale, there are variations to the regional trend, such as to the east of Devon Island and southeast Ellesmere Island, where temperatures at 200-m depth were lower in 2000–2010 relative to 1995–1999 (30). At latitudes south of 75°N, mean temperatures have been warmer than 0°C at 400- and 500-m depths since records began, and warming has occurred at all depths, except for a temperature decrease between 72°N and 74°N at 500-m depth (table S5). Although ocean warming has occurred throughout much of the region, the temperatures have been substantially below those seen to be driving melt in western Greenland and elsewhere [e.g., (16, 29)]. For example, the rapidly retreating Jakobshavn Isbrae is adjacent to the warm (mean, ~1.8°C) subsurface water of Disko Bay (16).

Fig. 5 Temporal trends from World Ocean Database (WOD13) data (up to 100 km offshore).

Mean ocean temperatures are separated by standard depths and degrees latitude for (A) QEI and (B) BBI. Error bars are ±1 SE.

The ice thickness and water depth at the grounding line also influence the degree of oceanic control on glacier retreat. Results from a recent study (31) show that the mean ice thickness at the grounding line of 35 tidewater glaciers in the QEI is estimated to be ~230 m and that only eight glaciers have a grounding line thickness that exceeds 300 m. In the BBI region, the mean ice thickness at the grounding lines is estimated to be even smaller (~100 m) (32). Therefore, the warmer water present at ocean depths deeper than 400 m will have little impact on the glaciers. By comparison, many outlet glaciers in Greenland extend into fjords to depths of 600 m or more, where they are exposed to deep warm Atlantic water (33, 34). Thus, while some large glaciers in the CAA may be vulnerable to melt in specific locations where they sit in deeper and warmer water, we find no evidence that ocean temperatures have driven the long-term regional patterns in the retreat of marine-terminating glaciers in the CAA.

Atmospheric temperatures and regional frontal change

Atmospheric temperatures in the CAA have increased rapidly in the 21st century (5, 35, 36), but the relationship between rising temperatures and the retreat rates of marine-terminating glaciers has not previously been investigated. To quantify this, we evaluate modeled climate data from the Regional Atmospheric Climate Model 2.3 (RACMO2.3), statistically downscaled from 11- to 1-km resolution (5). The downscaled data provide a detailed reconstruction of glacier SMB components for the period 1958–2015, corrected for biases in elevation and ice albedo (5). The 1-km product corresponds well with altimetry records and shows a close correlation between surface meltwater production and measured mean annual near-surface temperatures (5). From the downscaled version of RACMO2.3, we calculate mean ablation rates [millimeter water equivalent (w.e.) per year], mean ablation area ratio (as the percentage of basin size), and mean runoff (millimeter w.e. per year) for the full drainage basins of all 300 marine-terminating glaciers (see Materials and Methods and fig. S5).

Correlations between overall relative changes in glacier length (1958–2015) and the mean surface ablation components over the same period are all statistically significant [Spearman’s rank (rs), two-tailed, P < 0.01; table S6]. Glaciers that have retreated more are found to drain basins that have (i) higher mean ablation rates, (ii) larger ablation areas, and (iii) greater runoff, compared to those that have retreated less (Fig. 6). The strongest correlation (rs = −0.506, P < 0.01) is between overall frontal change and mean runoff. Overall, retreat rates are lowest north of 78°N in Nares Strait (Fig. 1 and fig. S2A), where runoff is the least (fig. S6), near-surface air temperatures are the lowest (5), and sea ice concentration remains high (37). There has been an increase in runoff from marine-terminating glaciers throughout the CAA, with the highest runoff occurring most recently (figs. S6 and S7 and table S7). The mean runoff at all latitudes during this latter period is over 40% greater than the mean before 2001. Few weather stations are situated in the region, but available records all show a trend toward more positive anomalies in mean summer near-surface air temperature (fig. S8).

Fig. 6 Mean overall (1958–2015) glacier front changes (as the percentage of glacier length in 2000) by mean basin runoff and relative ablation area over the same time period.

Glacier counts are in table S3. Error bars are ±1 SE. Statistically significant (P < 0.01) correlations are present between retreat rate and the surface ablation variables (table S6).

The temporal trends in runoff and ablation area extent between 1958 and 2015 correspond closely with the trends in frontal retreat rates in both the QEI and BBI regions (Fig. 2 and table S7). During concurrent periods, retreat rates are lower during times of lower runoff and smaller ablation areas and greater when runoff and ablation areas are larger. The only period for which this is not the case is 2002–2009 in the QEI, when the runoff and ablation areas decreased slightly when frontal retreat rates were increasing (although the direction of change is still consistent; Fig. 2, A and C). In BBI, the correlation persists across all six time periods, including the anomalous period of greater glacier retreat in 1980–1985, when there was a distinct increase in runoff/ablation (Fig. 2, B and D). A recent study of surface mass loss from 1-km downscaled RACMO2.3 output shows that ice masses in the CAA experienced significant near-surface atmospheric warming (+1.1°C above the 1958–1995 mean) in the mid-1990s (5), which caused increases in meltwater runoff across the region. The region-wide acceleration in retreat rates since ~2000 corresponds well with the modeled and measured increase in mean annual near-surface atmospheric temperatures. No lag effects are observed between the surface ablation rates and the frontal changes at the scale of the multiyear periods used in this study. This indicates a rapid and near-synchronous response of marine-terminating glacier fronts to atmospheric temperature changes throughout this region.


Our analysis indicates a rapid increase in the rate of glacier retreat in the CAA since the start of the 21st century. The evidence does not support ocean temperatures as the driver of the regional trends in glacier retreat along the western side of either Baffin Bay or Nares Strait, despite ocean forcing being recognized as a key driver of glacier frontal retreat on the eastern side of Baffin Bay in west Greenland. Both ocean temperature measurements and reanalysis data show predominantly cold waters, both on a broad regional scale (up to 100 km from the coast) and closer to glacier fronts (where data are available within 20 and 5 km). In contrast to west Greenland, CAA glaciers are relatively thin at their grounding lines (31) and tend to occupy shallower water. They are, therefore, less vulnerable to deep ocean temperature changes, and the mean ocean temperatures at shallower depths have been too low to have driven the observed region-wide patterns of retreat (29).

Results from the 1-km downscaled RACMO2.3, however, demonstrate a significant correlation between marine-terminating glacier front retreat rates and SMB components (mean ablation rate, ablation area ratio, and mean runoff) over the period 1958–2015. The relationship is quantitatively robust and occurs across the region. Moreover, temporal variations in retreat rates and surface runoff rates are largely synchronous, suggesting that there is a rapid response of CAA marine-terminating glacier fronts to near-surface atmospheric warming. External forcings on the fluctuations of marine terminating glaciers have been noted in varying degrees in other locations on Earth, but the regional dominance of atmospheric temperatures in driving marine-terminating glacier retreat has not been observed elsewhere.

Glacier frontal retreat rates have followed a similar trend to the rates of surface mass loss in the CAA, which have accelerated since ~2000 (1, 5). We suggest that the mechanisms by which increased surface air temperatures drive rapid retreat of marine-terminating glaciers in this region are likely to be much the same as those that drive retreat of the land-terminating glaciers. That is, increased melt and subsequent meltwater runoff reduces mass delivery to glacier termini through reductions in glacier flow, with the consequence that glaciers thin and retreat in their ablation areas (38). Other mechanisms observed on marine-terminating glaciers in Greenland, such as increased surface melt amplifying the delivery of ocean heat to the terminus via the creation of an upwelling plume at the grounding line (39), are less likely to be a primary mechanism behind regional glacier retreat in the CAA. This is because most of CAA glaciers terminate in shallow water, making them less vulnerable to intrusion of warmer deep water. However, this mechanism could be important on the local scale for the few CAA glaciers that terminate in deep water, although further research is required to examine this.

It is now widely acknowledged that ocean temperature increase has been the dominant driver of glacier retreat in other polar regions in recent years, particularly along the western Antarctic Peninsula, around the West Antarctic Ice Sheet, and in western Greenland. In contrast, we show that, in the CAA, the substantial rise in atmospheric temperature in the 21st century has outweighed any regional impact of changing ocean temperature on marine-terminating glacier frontal behavior. It follows that ocean temperature cannot be assumed to be the primary driver of marine-terminating glacier retreat in all polar regions and that studies of local processes are needed to understand the impacts of climate change on glacier behavior.


Glacier frontal positions data sources

The baseline image on which our study is based is the Natural Resources Canada (NRCan) image mosaic of the CAA, compiled from Landsat 7 scenes between 1999 and 2002, and orthorectified using ground control and conjugate points (40). The seamless mosaic is considered accurate to ~20 m, or less than one image pixel (30 m) (40). We downloaded the image tiles for our region of interest from the Canadian Open Government portal ( The glacier basin outlines were downloaded from the Global Land Ice Measurements from Space (GLIMS) database (41), with manual checking and correction of gross errors (e.g., coregistration of the basin outlines to the NRCan Landsat 7 mosaic). The glacier frontal positions for the period 1999–2002 were digitized from the mosaic and joined to the glacier basin outlines to produce closed polygons, with relevant attributes (such as glacier length, type, and ID) from the GLIMS database. Marine-terminating glaciers were identified as those that terminated in the ocean at the time of the Landsat 7 mosaic (1999–2002). In many cases, it was not possible to identify from the images the extent to which the glacier descended below sea level, and so, we used the general term “marine-terminating” rather than “tidewater” glacier. In total, we mapped all 300 marine-terminating glaciers in the CAA with basin size greater than 2 km2 (in 1999–2002), excluding those on the northern coast of Ellesmere Island (Fig. 1).

Glacier frontal positions were manually digitized, using Esri ArcGIS10.5 software from a range of sources (table S1). Aerial photographs, from survey missions between 1956 and 1961, were the first observations that covered the whole region. These were accessed from the NRCan National Air Photo Library. The glacier frontal positions were mapped from the aerial photographs by on-screen digitizing from scans of the photographs georeferenced to the Landsat 7 mosaic base image. Where original photographs were not available, the glacier fronts were digitized from scans of National Topographic System map sheets, which were compiled in the 1960s from the aerial photographs. An average of seven frontal positions, at approximately decadal intervals, were mapped per glacier. Post-1960, most glacier frontal positions were digitized from Landsat satellite images. These were accessed and downloaded using the U.S. Geological Survey Earth Explorer portal ( Image resolutions ranged from 5 m (aerial photographs) to 60 m (Landsat Multispectral Scanner), with all post-2000 images (Landsat Enhanced Thematic Mapper Plus and Landsat 8) at 15-m resolution. High relative positional accuracy was maintained by ensuring that the frontal positions were georegistered to the geospatially accurate Landsat 7 mosaic. The digitization methodology follows standard procedures when mapping glacier changes in previous work (42). Uncertainties in the glacier change results stem from the reliability of the source material, measurement error, and methods used to calculate change, as documented in Cook et al. (42). Similarly, in this study, it is estimated that the digitization uncertainties are less than the largest image pixel size (60 m). Uncertainty in the data analysis methods is described in the following section.

All glacier frontal data were assigned the relevant source material attributes including date, year, and source ID. In total, 2039 frontal positions were digitized for 211 glaciers in the QEI and 89 in BBI (fig. S1).

Glacier frontal change analysis

Change rate calculations. All frontal positions digitized for each basin were joined to the glacier basin outline, resulting in individual area polygons for each time period. Glacier area measurements take into account uneven changes along the ice front and offer a more precise assessment of change than can be obtained using along-flow sample lines that intersect the glacier fronts. Differencing of polygon areas from one epoch to another provides values of absolute area change. We used the established “box method” (43), with some minor modifications. Using ArcMap, a box polygon was digitized for each glacier, intersecting the glacier front at its maximum width and extending upflow to an undefined distance, ensuring that all decadal positions were contained within the box (fig. S1). Each frontal position was then clipped to the box and projected into equal-area projection (WGS84 Lambert Azimuthal Equal Area). The area of each clipped polygon (in square kilometers) enabled area change calculations between time periods, and the mean changes in length were calculated by dividing the area changes by the width of the glacier box. The rates of change (in meters per year) were then calculated by dividing the changes in length by the number of years that elapsed between the frontal positions. The magnitude of absolute glacier retreat is related to glacier size (larger glaciers tend to change more in absolute terms), so to enable comparisons between basins, we normalized changes according to glacier length. The relative length changes (as percentage per year) were calculated relative to the glacier length in 1999–2002 (using the length attribute reported in GLIMS).

Statistical analyses of regional and temporal glacier changes include mean absolute (meters per year) and relative (percentage per year) rates of frontal change for glaciers in the QEI and BBI over each time interval (Fig. 2 and table S2) and relative rates of frontal changes over each time interval at 1° latitudes (fig. S2). For all statistical analyses, mean values are considered a suitable measure, with SE bars to show the uncertainty around the estimate of the mean, considered the most valuable measure based on the large sample size. The CAA were split into the two broad regions of QEI and BBI for the following reasons: (i) clarity in data analysis, (ii) the notably different glacier basin sizes and frontal widths between the two regions, and (iii) the differences in image availability, which results in differing start/end years for the time periods.

Grounding line ice thickness. Ice thickness at the grounding lines of tidewater glaciers in the QEI was calculated by Van Wychen et al. (31) using data from a range of sources. Ice thickness values for the glaciers in the present study were obtained from Van Wychen et al. [Table S1 in (31)].

Ocean temperature data and methods

Ocean temperature data. 1) TOPAZ4 Arctic Ocean reanalysis data ( The TOPAZ4 Arctic Ocean Physics reanalysis data provide physical ocean variables for the time period 1991–2015. These data are available from the CMEMS (product identifier: ARCTIC_REANALYSIS_PHYS_002_003). The TOPAZ reanalysis system assimilates all available satellite and in situ temperature-salinity profile observations, resulting in modeled data as a 12.5-km resolution grid (25). We downloaded the mean monthly temperature and salinity NetCDF files at the subsurface standard depths available: 50, 100, 200, 400, and 700 m. Mean annual temperatures were calculated from the grid values at each depth. The quality information document ( gives the estimated accuracy for modeled temperatures as statistical values for bias and root mean square (RMS) differences. The forecast biases are −0.35°C at 5- to 100-m depth, −0.2°C at 100- to 300-m depth, and 0.5°C at 300- to 800-m depth. Forecast RMS differences are 0.6°C at 5- to 100-m depth, 0.8°C at 100- to 300-m depth, and 1.0°C at 300- to 800-m depth. These accuracy estimates are for the full domain, but there is variability in the system accuracy, e.g., in some regions with few temperature profile data, the reanalysis Root Mean Square error (RMSE) can be up to 1.0°C. Our statistical analysis does not include TOPAZ4 reanalysis data, and therefore, these errors do not affect our results. Rather, the data are used to illustrate the regional overview of mean annual ocean temperature at each depth (Fig. 3).

2) World Ocean Database in situ measurements ( The World Ocean Database (WOD13) (28) is a comprehensive set of historical ocean profile data collected between 1972 and 2015. Data are provided by the NOAA National Centers for Environmental Information. We selected in situ data from CTD profiles for the region surrounding the CAA. We downloaded temperatures at the following standard depths available: 50, 100, 150, 200, 300, 400, and 500 m. WOD13 water properties at standard depths were measured and reported either directly or by interpolation, as set out by commonly used criteria (28). The data variables include the measurement station ID, latitude/longitude, date, and year, in addition to temperature and salinity values. To download the data for our region of interest, we used the Ocean Data View software provided by NOAA. Error estimates described in the WOD13 report (28) state that the accuracy of CTD measurements varies with instrument design, typically from 0.005° to 0.001°C for temperature. The overall quality of CTD measurements are also affected by calibration drift (i.e., prolonged use of the CTD instrument in the sea), data streaming problems, varying speeds of the CTD, and rapid changes in the ocean environment. The measurements are assigned a flag in quality control checks and measurements that are considered to fail quality control are excluded from WOD13. Other flags range from good quality (1), to potentially bad quality but based on visual inspection found no reason not to be used in scientific research (4) (28). For our study, we included all measurements available in WOD13, and our representation of error in the statistical analysis methods is SE, considered the most valuable owing to the large sample size.

Bathymetric data ( The standard depths at which the ocean temperature measurements were chosen were based on the range of depths close to the glacier fronts. These were established from the NOAA International Bathymetric Chart of the Arctic Ocean v3, 500-m grid (fig. S3). The grid is compiled from all available multibeam, dense single beam, and land data and built using a gridding algorithm at 500-m resolution (44). It is currently recognized as the best Arctic Ocean Digital Bathymetric Model available. We downloaded and clipped the dataset to our region of interest and removed land data before undertaking statistical analysis of depths close to glacier fronts.

Buffers extending from the marine-terminating glacier fronts (at their positions ca. 2000) were created at 5- and 20-km radii (fig. S3). The mean bathymetry within the 20-km buffer for all 300 glaciers is 166 m of which 40 glaciers have mean bathymetry deeper than 400 m. Within 5 km of all glacier fronts, the mean bathymetry is 90 m, and 11 glaciers have a mean bathymetry deeper than 400 m within 5 km of their fronts. Of the 117 glaciers with ≥10 ocean temperature measurements within 20 km, 36 have a mean depth deeper than 400 m (mean, 477 m). When all glaciers are separated by 1° latitudes, the mean depths within 20 km are shallower than 400 m in all latitudinal degree bands, and within 5 km, mean depths are shallower than 200 m.

Ocean temperature statistics. Regional modeled mean ocean temperatures were mapped in ArcMap using TOPAZ4 gridded data for the region surrounding the CAA and west Greenland (Fig 3). Ocean temperature analysis close to glacier fronts was carried out using in situ data from the WOD13 within the buffers at 20- and 5-km radii from each glacier front (fig. S3). In total, 117 glaciers have ≥10 CTD measurements within 20 km of their fronts, and 19 glaciers have ≥5 CTD measurements within 5 km of their fronts. Measurements within both radii have similar mean temperatures at each depth (fig. S4). Glaciers with ocean measurements within 20 km provided a wide-ranging sample for statistical analysis. Our investigations of glacier front change and ocean temperature correlations include glacier relative change value bins and mean ocean temperatures at each depth (glacier counts are shown in table S3). Spatiotemporal trends are based on WOD13 ocean temperature data within the 100-km buffer zone (Fig. 5). The region size was broken down by 1° latitudinal bands for regional mean ocean temperature analysis (fig. S3).

Atmospheric temperature data and methods

RACMO2.3 downscaled to 1-km data source. A detailed reconstruction of SMB between 1958 and 2015 is available from the RACMO2.3 statistically downscaled to 1 km. The model is a downscaled version of the 11-km RACMO2.3, corrected for biases in elevation and bare ice albedo, and provides regional SMB data at a scale suitable for measuring mean values within individual glacier basins in the CAA. The SMB output is in millimeter w.e. year−1.

The methods used in producing the model, along with accuracy assessments based on in situ mass balance measurements, are reported in Noël et al. (5). These are summarized as follows. The SMB values from both the 11- and 1-km versions of RACMO2.3 were assessed against 4356 stake measurements within the CAA. Correcting for elevation and bare ice albedo bias during the downscaling process significantly improved the accuracy. Assessment with observational data showed that, relative to the 11-km model, the 1-km downscaled version of RACMO2.3 has a much improved mean SMB bias (40 mm w.e.), RMSE (380 mm w.e.), and R2 (0.74). The product uncertainty for SMB was estimated at 11.0 and 4.5 Gt year−1 for Northern CAA (NCAA = QEI) and Southern CAA (SCAA = BBI), respectively. To obtain these estimates, the 4356 in situ SMB measurements from 1961 to 2016, e.g., from Sverdrup Glacier on Devon Island (yellow dots in fig. S5), were binned in 250 mm w.e. intervals for which the mean bias was estimated (i.e., modeled minus measured SMB), which was then applied to total ablation and accumulation areas across the region (5). The RACMO2.3 product was further evaluated using independently derived mass changes from ICESat and CryoSat-2 altimetry measurements. The downscaled dataset shows good agreement with the altimetry records, and both winter accumulation and summer ablation are accurately resolved. In addition, the SMB product falls well within ICESat/CryoSat-2 measurement uncertainty. Ice mass loss estimates from the downscaled SMB product are similar to estimates from other mass loss studies (3, 6, 32). Mass loss in 2004–2009 was estimated at 33.8 ± 11.5 Gt year−1 for the NCAA and 24 ± 4.5 Gt year−1 for the SCAA (5). For comparison, a study by Gardner et al. (3) estimated mass loss from GRACE (Gravity Recovery and Climate Experiment) to be 39 ± 9 and 24 ± 7 Gt year−1 for NCAA and SCAA, respectively.

SMB calculations per glacier basin. For our study, the SMB calculations per glacier basin were made from the downscaled RACMO2.3 by converting the gridded values to points for analysis in ArcGIS. The 1-km resolution data enabled statistical analysis of ablation components within each glacier basin (fig. S5).

The mean ablation rates per glacier basin were calculated from negative (i.e., <0) SMB values, in mm w.e. year−1. Net ablation occurs in the region of the glacier basin where ice loss through surface melt, sublimation, and runoff is greater than ice gain through precipitation and riming. The ablation area ratio (in percentage) was therefore calculated as the area where SMB < 0 mm w.e. year−1, divided by the glacier basin area. This is calculated from the gridded data, i.e., the number of negative SMB points divided by the total number of grid points that constitute each glacier basin. The ablation area ratio enables analysis of patterns of change across all glaciers by normalization of basin areas. Mean annual runoff (in mm w.e. year−1) represents basin-wide nonrefrozen rain and meltwater induced by the surface energy balance, and we calculated this from the sum of the runoff component values divided by the number of grid points per basin. It should be noted that summer runoff consist partly of rainfall, which does not contribute to mass loss, but the flux is small [Figure 8 in (5)].

We calculated mean values of these components across the whole study period (1958–2015) and, for each time period, used in the glacier front change study. Mean runoff patterns show similar temporal trends for marine-terminating glaciers at all degrees latitude, with a statistically significant increase in recent time periods (fig. S6). This is reflected in the spatiotemporal patterns of runoff from ice masses throughout the CAA (fig. S7).

Automatic weather station data. Near-surface (2 m) air temperature data are available from weather stations operated by Environment Canada from 1948 to the present ( There are six weather stations located throughout the CAA (locations shown in Fig. 1), and although some are situated at a considerable distance from the marine-terminating glaciers in this study, the records provide a useful indication of long-term temperature trends in the wider region. Temperature measurements are available at the following weather stations in the QEI: Eureka (1947–2014), Resolute (1947–2014), and Grise Fiord (1984–2007); and in BBI: Pond Inlet (1975–2007), Cape Hooper (1957–2007), and Cape Dyer (1959–1993). Mean summer air temperatures, for the months of June to August (JJA), were calculated from monthly mean temperatures derived from daily averages. Annual anomalies in mean summer temperatures were calculated relative to the mean temperature across the complete time period available for each weather station (fig. S8).


Supplementary material for this article is available at

Fig. S1. Sample glaciers on southeast Devon Island, illustrating the database compilation and analysis methods.

Fig. S2. Mean relative change per year (as the percentage of glacier length in ~2000), separated by latitudinal degrees.

Fig. S3. Ocean temperature in situ measurement locations and marine-terminating glacier buffer limits.

Fig. S4. Mean ocean temperatures from WOD13 measurements available within glacier front buffer regions.

Fig. S5. Modeled climate data accuracy assessment on Sverdrup Glacier, on north Devon Island (75.5°N, 83°W).

Fig. S6. Change in surface meltwater runoff over time (mm w.e. a−1).

Fig. S7. Mean runoff calculated from RACMO2.3 statistically downscaled to 1 km resolution for glaciers, ice caps and ice fields in the CAA.

Fig. S8. Temperature trends from six weather stations in the CAA.

Table S1. Image sources from which CAA glacier frontal positions were mapped.

Table S2. Mean absolute and relative length change rates in each time period, for glaciers in the CAA.

Table S3. Glacier counts for overall (1958–2015) glacier length change bins (see Figs. 4 and 6).

Table S4. Spearman’s rank (rs) correlations between mean relative change rates and mean ocean temperature at each depth.

Table S5. Ocean temperature measurements descriptive statistics.

Table S6. Spearman’s rank (rs) correlations between overall (1958–2015) glacier front change rates (as the percentage of basin length) and mean SMB components.

Table S7. Mean glacier runoff (mm w.e. a−1) and mean frontal change rates (% a−1) for glaciers separated by latitudinal degrees and time periods for BBI and QEI.

Data file S1. Marine terminating glaciers (>2 km2 in size) on Baffin and Bylot Islands (BBI): data used in statistical analysis.

Data file S2. Marine terminating glaciers (>2 km2 in size) on the Queen Elizabeth Islands (QEI): data used in statistical analysis.

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.


Acknowledgments: We would like to thank NRCan, GLIMS, U.S. Geological Survey, CMEMS, the National Centers for Environmental Information, the Government of Canada, and the World Glacier Monitoring Service for access to datasets. We would also like to thank P. Whitehouse for processing the CMEMS data. We thank two anonymous reviewers for their insightful suggestions that improved the manuscript. Funding: A.J.C. was supported by a Leverhulme Early Career Fellowship and the Department of Geography, Durham University. Additional funding was provided by the Natural Sciences and Engineering Research Council of Canada, University of Ottawa, and Canada Foundation for Innovation. B.P.Y.N. and M.R.v.d.B. acknowledge funding from the Netherlands Earth System Science Centre (NESSSC) and the Polar Program of the Netherlands Organisation for Scientific Research (NWO/NPP). Author contributions: A.J.C. led the analysis, writing, and compilation of graphics and tables. L.C. contributed image data and expertise, and B.P.Y.N. and M.R.v.d.B. contributed atmospheric climate model data and analysis. C.R.S., M.J.B., M.J.S., and R.G.B. designed the study and contributed knowledge. All of the authors contributed to the writing of the manuscript. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors. Data are publically available from the Polar Data Catalogue (

Stay Connected to Science Advances

Navigate This Article