Retreat of the Greenland Ice Sheet leads to divergent patterns of reconfiguration at its freshwater and tidewater margins

Abstract

Greenland's marine- and land-terminating glaciers are retreating inland due to climate warming, reconfiguring the way the ice sheet interacts with its proglacial environment. Here we use three decades of satellite imagery to determine whether the ice-sheet margin is becoming more or less exposed to marine and lacustrine processes. During our 1990–2019 study period, we find that the length of ice-sheet perimeter in contact with the ocean shrank by 12.3 ± 3.8% (196.2 ± 10.4 km), due to the retreat of marine-terminating glaciers into narrower fjords. On the other hand, we find that the length of the ice-sheet perimeter in contact with freshwater lakes exhibited more divergent trends that is better explored at regional scales. The length of ice–lake boundaries increased in southwest, north and northwest Greenland but declined in southeast and central east Greenland. The magnitude of change we document during our study period leads us to conclude that the ice sheet is poised for further, substantial reconfiguration in the coming decades with consequences for the flux of fresh water, nutrients and primary productivity in Greenland's terrestrial and oceanic environment.

1. Introduction

Nearly 20% of observed global sea-level rise (2006–18) can be attributed to mass loss from the Greenland Ice Sheet, which has been in a state of negative mass balance for the past several decades (van den Broeke and others, 2016; Otosaka and others, 2023). In response to this mass loss, the margins of the ice sheet have undergone rapid and substantial retreat (Moon and others, 2020). The most dramatic rates of retreat have been observed at the Greenland Ice Sheet's marine-terminating outlet glaciers which have retreated a few tens to several hundreds of meters on average per year over the last four decades (~1980s–2020s) (Carr and others, 2017; King and others, 2020; Fahrner and others, 2021; Kochtitzky and Copland, 2022). Land-terminating glaciers have also receded, albeit at slower rates compared to their marine-terminating counterparts (Mernild and others, 2012; Mallalieu and others, 2021). These general patterns of retreat are reshaping the margins of the ice sheet with consequences for the flux of fresh water, sediment and nutrients into Greenland's terrestrial and oceanic environment (Storms and others, 2012; Meire and others, 2017; Stuart-Lee and others, 2021, 2023).

The character of reconfiguration at the ice-sheet margin in response to retreat mainly depends on the unique topographic setting of each glacier. Many glaciers have floating ice tongues or shelves which have decreased in area over the last few decades (Moon and Joughin, 2008; Box and Decker, 2011; Hill and others, 2018). One of the most famous examples is Petermann Glacier which lost 277 km² between end-of-summer periods in 2009 and 2010 (Box and Decker, 2011). Other marine-terminating glaciers (e.g. Rink Glacier) have retreated into narrower fjords (Moon and others, 2020), with at least three retreating so far inland that they have abandoned their fjords and now terminate on land (Williams and others, 2021; Kochtitzky and Copland, 2022). The combination of shrinking floating tongues and retreat of glaciers inland reduces the amount of ice in contact with the ocean which may subsequently alter circulation, nutrient availability and primary productivity in Greenland fjords (Meire and others, 2017; Stuart-Lee and others, 2021, 2023).

Land surface topography also determines how retreating land-terminating glaciers interact with their proglacial environment. Much of the meltwater produced on the ice sheet exits through subglacial portals, which drain into proglacial rivers that evacuate meltwater away from the ice sheet toward the ocean. However, if land-terminating glaciers retreat into topographic depressions (overdeepenings), meltwater may be temporarily stored in proglacial lakes (Carrivick and Tweed, 2013). Meltwater stored in these lakes can enhance the rate of mass loss through increased melting and calving at the glacier terminus, especially during periods of rapid water level fluctuation (Carrivick and Tweed, 2019; Mallalieu and others, 2020). An increase in ice-sheet area exposed to proglacial lakes could therefore result in an unanticipated acceleration in glacier mass loss since this process is currently not accounted for in ice-sheet models (Mallalieu and others, 2021).

A number of studies have documented multi-decadal changes occurring at the Greenland Ice Sheet's oceanic margins (Murray and others, 2015; Carr and others, 2017; Hill and others, 2018; King and others, 2020; Fahrner and others, 2021; Goliber and others, 2022). However, these studies tend to focus on calving front retreat and ice-flow velocity since they are often coupled and have implications for ice-sheet contributions to global sea-level rise. Fewer studies have quantified changing calving front widths. Those that have only investigated a regional subset of the largest marine-terminating glaciers (e.g. Seale and others, 2011; Catania and others, 2018; Fried and others, 2018) or only focused on one time period (e.g. Carrivick and others, 2022). To our knowledge, changes in the total length of the ice-sheet perimeter in contact with the ocean have yet to be reported.

Studies have also documented changes occurring at Greenland's proglacial lake margins but long-term trends are inconclusive (Carrivick and Quincey, 2014; How and others, 2021). Carrivick and Quincey (2014), for example, identified a 20 ± 6.5% expansion in total lake surface area between 1987 and 2010 in western Greenland. But How and others (2021) found that the area of proglacial lakes did not significantly change between 1985–1988 and 2017 in the same region. Evidence for a progressive increase in proglacial lake area therefore remains disputed in southwest Greenland and long-term trends in other regions remain unknown. The criteria for identifying proglacial lakes also remain unsettled. Many studies define ‘ice-marginal’ lakes as those adjacent to the ice sheet (Carrivick and Quincey, 2014; How and others, 2021). For example, How and others (2021) identified lakes within <1 km from the ice-sheet margin. Yet, not all of these lakes are in contact with the ice sheet so do not have a direct impact of mass loss. Fewer studies have mapped lakes that are in direct contact with the ice sheet, and those that have either focused on a single time period (e.g. Carrivick and others, 2022) or a specific region (e.g. Mallalieu and others, 2021). No study has mapped the total length of the Greenland Ice Sheet perimeter in contact with the ocean and lakes for multiple time periods. Given the enhanced retreat that occurs at ice–water boundaries relative to ice–land boundaries, the lack of such information impedes any process-based analysis required to investigate current and future patterns of ice-sheet reconfiguration and associated impacts on mass balance and proglacial environments.

Here, we use satellite remote sensing to investigate how general patterns of ice-sheet retreat since the 1990s have altered the way the ice sheet interacts at its lacustrine and oceanic boundaries. We first delineate ice–lake and ice–ocean boundaries for three distinct time periods (1990–95, 2003–07 and 2019) using an edge detection method applied to Landsat 5, 7 and 8 imagery (see Section 2). We then assess how these boundaries have changed at ice sheet, regional and individual glacier scales. Finally, we discuss the implications of our findings for future reconfiguration of the Greenland Ice Sheet margin and potential consequences for ice-sheet mass loss.

2. Methods

2.1 Production of classified Greenland Ice Sheet image mosaic

We produced a cloud-free image mosaic of Greenland using Landsat 5™, Landsat 7 ETM+ and Landsat 8 OLI surface reflectance data for three distinct time periods (1990–95, 2003–07 and 2019) using Google Earth Engine (GEE). Different length time periods were required to ensure full image coverage of the ice sheet. We only included Landsat images acquired between 1 June and 31 August with <20% total cloud cover. We note that limiting the acquisition of images to the summertime reduces bias due to seasonal fluctuations in ice margins. The image mosaics for each period were produced by computing the median value of each pixel which removed artifacts due to clouds and contrails. Despite extending the earliest time period to 6 years, we were still not able to produce a complete ice-sheet mosaic for the 1990–95 period. We therefore masked the 1990–95 missing regions from the 2003–07 and 2019 mosaics during multi-decadal trend analysis (see Section 2.7).

To classify the image mosaics, we generated separate training datasets for each mosaic. Each training dataset contained ~1000 manually labeled data points, representing the surface reflectance values for water, ice and land in bands 1–7 for Landsat 5 and Landsat 7 and bands 2–6 for Landsat 8. Water, ice and land pixels were selected from all regions of Greenland to produce a training dataset that was representative of the ice sheet and surrounding margins. We used these labeled datasets to train three separate random forest classification models (one for each time period) with ten decision trees (which were found to be an optimal tradeoff between classification accuracy and computational efficiency). In total, 80% of the manually classified data was used for training the model while 20% was reserved for testing model accuracy. We achieve overall accuracies of 99.4% for the 1990–95 mosaic, 99.4% for the 2003–07 mosaic and 99.5% for the 2019 mosaic. We then applied the models to the corresponding image mosaic and exported the classified mosaics in GeoTIFF format for further analysis.

2.2 Pre-processing

After classification, our classified image mosaics contained ‘salt-and-pepper’ effects due to sparsely occurring, often singular, misclassified pixels. As a first step toward correcting these misclassified pixels, we applied a median filter with a 7 × 7 (210 m × 210 m) window. This step removed most of the ‘salt-and-pepper’ effects while maintaining distinct boundaries between classes. Next, we noticed that there were pixels classified as ice that were detached from the main ice sheet. To remove these pixels, we grouped all ice values together using connected-component labeling (8 pixels) and removed all groups that contained fewer than 60 000 pixels. We found that this threshold represented a good compromise between removing peripheral glaciers and ice caps while retaining the pixels that are connected to the ice sheet. We also note that our analysis is not particularly sensitive to this threshold because of the manual quality checks we perform later in the analysis.

2.3 Detection of boundaries between ice sheet and surface water

We detected boundaries between surface types by converting our classified mosaic into a binary mask where pixels classified as ice were assigned a value of 1 and all other pixels (i.e. land and water) were assigned a value of 0. We then used a Sobel filter to detect edges between these two classes. For each pixel classified as an ‘edge’ we determined whether the edge was an intersection between ice and water by inspecting the surface types of pixels either side of the edge using a 5 × 5 (150 m × 150 m) window. If the values in the 5 × 5 window surrounding the edge did not contain any land pixels, contained between 1/3 and 2/3 water pixels, and contained between 1/3 and 2/3 ice pixels, the edges were classified as boundaries between ice and water. We tested several different thresholds but found that the 1/3 and 2/3 thresholds reliably detected ice–water boundaries without being oversensitive to a few ice (or water) pixels. We vectorized the pixels labeled as ice–water edges to polylines and merged all features to produce a single polyline dataset for each time period.

To reduce computation time during this stage of the analysis, we only conducted our edge detection analysis within 10 km of the ice-sheet perimeter. To accomplish this, we produced a buffered mask using a polyline shapefile of the ice sheet's extent produced by PROMICE (Citterio and Ahlstrøm, 2013). Since the dataset distinguishes local ice masses from the main ice sheet, we first filtered the dataset to only include geometries of the main ice sheet. We then buffered the polyline by 10 km. This threshold was chosen since it ensured that we incorporated all potential boundaries of the ice sheet with proglacial lakes or the ocean.

2.4 Manual filtering

During our automated analysis, we aimed to reduce error of omission as much as possible so as not to miss any potential ice–lake or ice–ocean boundaries. As a result of these decisions, our error of commission was higher. But it was deemed more efficient and accurate to manually remove misclassified boundaries rather than manually attribute undetected boundaries. After vectorizing the boundaries between ice and water, we visually inspected every boundary and manually edited them, where necessary, in QGIS using the same true color Landsat image mosaics from GEE as a reference. Boundaries that were incorrect (e.g. between supraglacial lakes and ice) were removed and boundaries containing small errors were edited to match the satellite imagery. Although this manual approach was time-consuming, we argue that it is a prerequisite for producing a high-quality product. We also manually assigned a unique identifier to all boundaries so that we could track their dynamics between 1990–95 and 2019.

2.5 Classification of ice–lake vs ice–ocean boundaries

Since our method detects all boundaries between the ice sheet and surface water, we separated the boundaries into ‘ice–ocean’ and ‘ice–lake’ using the TermPicks dataset which contains the terminus positions for 278 marine-terminating glacier margins (Goliber and others, 2022). To do this, we first filtered the TermPicks dataset (Goliber and others, 2022) to only include terminus positions for time periods that correspond with our study period. We then buffered these terminus positions by 200 m and used a spatial join to detect intersections between the buffered TermPicks dataset and our ice–water boundaries. Boundaries that intersected with the buffered terminus positions data were assigned an ‘ocean’ attribute, while those which did not intersect were given a ‘lake’ attribute. During the manual editing stage (see the previous section), we identified 105 additional ice–ocean boundaries that were not contained in the TermPicks dataset. We therefore manually assigned these ‘ocean’ values.

2.6 Quantification of boundary length

In order to quantify the length of the ice sheet that interacts with surface water consistently we converted the detected ice–lake and ice–ocean boundaries into multipoint features with ten vertices per km. We note that ten vertices per km (i.e. 100 m node spacing) is similar to that recommended by Goliber and others (2022) for manual digitization and length assessment of calving front termini. It is also justified by the mean length of ice–ocean and ice–lake boundaries which is 4.2 ± 8.3 and 0.8 ± 0.9 km (mean ± std dev.), respectively. We then calculated the total length of the ice–ocean and ice–lake boundary for each time period. To eliminate distortions in length due to projection errors, we reprojected each boundary into a unique local orthographic coordinate system before computing the length of the boundary following Kochtitzky and Copland (2022). We note here that the change in length of 5.5% of ice–water boundaries (6.7% of total ice–water length) was not tracked because they were not contained in the 1990–95 mosaic. Regional statistics were based on the catchment areas defined by Mouginot and Rignot (2019). For the total perimeter, we used the same PROMICE perimeter polyline that was used to buffer our classified mosaics in our statistical analysis (Citterio and Ahlstrøm, 2013). For consistency, we also converted this polyline into multipoint features with ten vertices per km. We manually removed perimeter points (e.g. around nunataks) that were not directly at the outer margin of the ice sheet. Our ice-sheet perimeter is therefore 29 269 km in length.

2.7 Uncertainty

Given the extensive manual correction applied to our dataset, we assume that uncertainties in boundary length are similar to those typically encountered when manually digitizing a satellite image. These uncertainties are usually considered to be ~1 pixel (±30 m in our case) for each segment in both the x and y directions (i.e. \(30\;\sqrt {2\;}\)). We therefore define the uncertainty in the length of each boundary as the root sum squared (RSS) of the uncertainty of each segment. Boundary length uncertainty is therefore proportional to the length of the boundary (or number of segments). Likewise, we define the uncertainty of regionally aggregated boundary lengths as the RSS of the individual boundary length uncertainties. This assumes that boundary length uncertainties are independent which is, at least partly, confirmed by low bias between our boundary lengths and manually traced ice–ocean boundary lengths provided by TermPicks dataset (Fig. S1).

3. Results

In the earliest period of our study (1990–95), we identify a total of 379 individual ice–ocean boundaries with an average length of 4.2 ± 8.3 km (mean ± std dev.). The ice–ocean boundaries have a total length of 1597.7 ± 5.4 km (mean ± uncertainty; Table 1), equivalent to 5.5 ± 0.1% of the perimeter of the Greenland Ice Sheet (modified from Citterio and Ahlstrøm, 2013). We identify more ice–lake boundaries (n = 659) but these have a smaller average length of 0.8 ± 0.9 km (mean ± std dev.) than the ice–ocean boundaries and therefore comprise a smaller proportion of the ice-sheet perimeter (546.4 ± 3.1 km or 1.9 ± 0.05%; Table 1).

Table 1. Summary of ice–ocean and ice–lake boundaries of the Greenland Ice Sheet for the 1990–95, 2003–07 and 2019 periods, discriminated by region

		1990–95		2003–07		2019
Region	Perimeterkm	Ice–ocean lengthkm	Ice–lake lengthkm	Ice–ocean lengthkm	Ice–lake lengthkm	Ice–ocean lengthkm	Ice–lake lengthkm
N	2548.9	204.5	48.8	184.1	59.4	180.2	67.3
NE	5641.0	321.2	103.5	269.4	107.8	201.7	111.6
CE	8223.2	269.9	32.1	217.8	31.8	229.1	20.4
SE	4045.8	217.4	21.5	206.5	13.1	235.6	11.4
SW	4131.4	51.6	237.0	51.5	220.5	41.5	199.2
CW	1296.9	91.2	77.3	86.4	55.3	97.1	62.4
NW	3381.8	442.0	26.2	422.9	43.5	416.2	45.9
Total Greenland Ice Sheet	29 269.0	1597.7	546.4	1438.6	531.3	1401.5	518.2

We identify ice–ocean and ice–lake boundaries in all regions of the ice sheet (Fig. 1), but their distribution varies substantially. By categorizing our findings into subregions defined by Mouginot and others (2019), we find that the northwest (NW), north (N) and central west (CW) sectors of the ice sheet have the largest proportion of ice–ocean boundaries (13.1 ± 1.0, 8.0 ± 0.7 and 7.0 ± 0.8% of the total perimeter of the ice-sheet margin, respectively; Fig. 1). In contrast, ice–ocean boundaries occupy a smaller proportion of the southwest (SW) and central east (CE) sector perimeters (1.2 ± 0.1 and 3.3 ± 0.2%, respectively). We find that regional patterns of ice–lake boundaries deviate from those of ice–ocean boundaries (Fig. 1). The largest proportions of ice–lake boundaries are found in CW and SW Greenland (6.0 ± 0.7 and 5.7 ± 0.4%, respectively) whereas smaller proportions are found in CE, SE, and NW Greenland (0.4 ± 0.02, 0.5 ± 0.04 and 0.8 ± 0.06%, respectively; Table 1).

Figure 1.

Location and trends in the length of ice–ocean and ice–lake boundaries for the Greenland Ice Sheet. Dots show location of ice–ocean (blue) and ice–lake (red) boundaries. Rings are proportional to length of ice–ocean and ice–lake boundaries in 2019. Bars show change for our three study periods: 1990–95, 2003–07 and 2019. Error bars are too small to show. Regional abbreviations correspond to: N, north; NE, northeast; CE, central east; SE, southeast; SW, southwest; CW, central west; NW, northwest Greenland.

[Figure omitted. See PDF]

We find that the length of the ice sheet in contact with the ocean decreased during our 1990–2019 study period (Fig. 1). Across the entire ice sheet, the ice–ocean boundary length decreased by a total of 196.2 ± 10.4 km, equivalent to a 12.3 ± 3.8% reduction of the ice-sheet margin in contact with the ocean. Regionally, the greatest decreases in ice–ocean boundary length occurred in NE and CE Greenland where −119.5 ± 3.1 and −40.8 ± 3.0 km of the ice–ocean boundary was lost over our study period, respectively. Substantial reductions have also occurred in N (−24.3 ± 2.6 km), SW (−10.0 ± 1.3 km) and NW (−25.8 ± 3.9 km) Greenland. In contrast, the ice–ocean boundary length increased in SE and CW Greenland (+18.3 ± 2.9 and +5.9 ± 1.8 km, respectively) during our study period. We explore the causes of these changes later in this section.

The total length of ice–lake boundaries exhibited more divergent trends during our 1990–2019 study period. Across the entire ice sheet, the length of ice–lake boundaries slightly decreased during our study period. However, it should be noted that the magnitude of change (−28.2 ± 4.4 km) is much less than the decrease in total ice–ocean boundary length. At the regional scale we find reductions in the ice–lake boundary length in SW (−37.8 ± 2.8 km), CW (−14.8 ± 1.6 km), CE (−11.7 ± 1.0 km) and SE (−10.1 ± 0.8 km) between 1990 and 2019. In contrast, we find increases in ice–lake boundary length in the N (+18.6 ± 1.4 km), NW (+19.6 ± 1.1 km) and NE (+8.1 ± 2.0 km) sectors of the ice sheet.

By tracking ice–ocean boundaries between 1990–95 and 2019, we find that 62.8 ± 16.8% decreased in length for a total of −344.7 ± 70.1 km. This was partly compensated by a +151.0 ± 48.0 km increase in ice–ocean boundary length over the same time period for a net loss of −193.3 ± 85.0 km. This implies that almost all (~99%) of the total reduction in ice–ocean boundary length that we observe (196.2 ± 10.4 km) is characterized by the retreat of marine-terminating glaciers into narrower fjords, as opposed to their retreat onto land. To illustrate, we describe some representative case studies of marine-terminating glaciers retreating into narrower fjords. Apuseeq Anittangasikkaajuk (or Storbræ), for example, is known to be one of the fastest retreating glaciers in CE Greenland during the satellite era (Jiskoot and others, 2012). In 1990–95, this glacier had a 3.5 ± 0.3 km wide calving front which decreased to 1.9 ± 0.2 km by 2019 in response to persistent thinning during the early 2000s and retreat into a narrower part of the fjord (Jiskoot and others, 2012; Ultee and Bassis, 2020) (Fig. 2a). Similarly, Hagen Bræ has exhibited some of the largest rates of retreat of glaciers in N Greenland. Hagen Bræ is a surge-type glacier with a floating ice tongue (Rignot and others, 2001). At the beginning of our study period, the floating tongue was grounded on a small island in the middle of the fjord which acted as a pinning point (Joughin and others, 2010). After subsequent thinning in the early 2000s, the ice tongue fractured and lost contact with the island, leading to a dramatic loss of floating ice (Fig. 2b; Murray and others, 2015). We find that the length of ice tongue in contact with the ocean decreased from 16.5 ± 0.5 km in 1990–95 to 11.3 ± 0.5 km in 2019.

Figure 2.

Reconfiguration of the Greenland Ice Sheet's ice–ocean interface during the periods 1990–95, 2003–07 and 2019. North is up and locations can be found in Figure 1. (a) The calving front of Apuseeq Anittangasikkaajuk (CE) shrank from 3.5 ± 0.3 km in 1990–95 to 1.9 ± 0.2 km in 2019 as it retreated into a narrower portion of the fjord. (b) The floating ice tongue of Hagen Bræ (N) fractured during our study period, reducing in length from 16.5 ± 0.5 km in 1990–95 to 11.3 ± 0.5 km in 2019. (c) Morell Glacier (NW) split from one to three calving fronts causing an increase in total length from 3.6 ± 0.3 km 1990–95 to 6.3 ± 0.3 km in 2019. (d) Nigertiip Apusiia (SE) also split from one to many calving fronts causing an increase in total length from 3.8 ± 0.3 km in 1990–95 to 4.0 ± 0.3 km in 2019.

[Figure omitted. See PDF]

A process that complicates the general reduction in ice–ocean boundary is the splitting of single marine-terminating glacier calving fronts into multiple calving fronts. For example, Morell Glacier in NW Greenland, split from a single calving front with a length of 3.6 ± 0.3 km in 1990–95 into three distinct calving fronts with a combined length of 6.3 ± 0.3 km in 2019 (Fig. 2c). Likewise, Nigertiip Apusiia (or Midgård Glacier) in SE Greenland has been characterized by a persistent negative mass balance in the 21st century, retreating at rates of >500 m a⁻¹ between 2000 and 2010 (Howat and Eddy, 2011; Walsh and others, 2012; Williams and others, 2021). As Nigertiip Apusiia retreated, it split into four distinct calving fronts in 2010–11 which have a greater combined length than the original single calving front (Fig. 2d). We find 23 marine-terminating glaciers that split between 1990–95 and 2019 (Table S1). After splitting, the length of the ice–ocean boundary for most (74%) of these glaciers increased. Overall, the splitting of marine-terminating glaciers caused a +16.0 ± 2.2 km increase in the total ice–ocean boundary and may have therefore obscured more general patterns of ice–ocean boundary loss during our study period.

At the ice sheet's land-terminating glaciers, we find that retreat can both increase and decrease the exposure to fresh water. Between 1990–95 and 2019, the ice-sheet margin just south of Sermeq Kujalleq retreated by 1.0–2.5 km inland, causing the expansion of several proglacial lakes which increased the length of ice–lake boundary from 3.6 ± 0.3 to 5.1 ± 0.3 km (Fig. 3a). Likewise, an unnamed glacier 60 km west of C.H. Ostenfeld Glacier in N Greenland was land-terminating in the period 1990–95 with a small ice–lake boundary of 0.3 ± 0.1 km (Fig. 3b). Over the last few decades, the valley has flooded with meltwater. By 2007 the meltwater-filled lake reached the terminus of the two glaciers that terminate in the fjord and, by 2019, the length of ice–lake boundary increased to 3.6 ± 0.3 km which is equivalent to a 5.3% increase in the total ice–lake boundary in the region (Fig. 3b). In contrast, Akuliarutsip Sermia in CW Greenland has experienced a decrease in ice–lake boundary (Fig. 3c). Retreat from an island into a narrower part of the basin has reduced the length of the ice in contact with lake from 5.1 ± 0.3 to 3.9 ± 0.3 km even though the overall area of the lake has actually increased from 39.8 to 42.2 km².

Figure 3.

Reconfiguration of Greenland's land-terminating glaciers during the periods 1990–95 and 2019. North is up and locations can be found in Figure 1. (a) As the ice-sheet margin just south of Sermeq Kujalleq retreated by 1.0–2.5 km, the total length of ice–lake boundary increased from 3.6 ± 0.3 to 5.1 ± 0.3 km. (b) An unnamed glacier to the west of C.H. Ostenfeld Glacier has experienced a large increase in freshwater interaction (from 0.3 ± 0.1 km in 1990–95 to 3.6 ± 0.3 km in 2019) due to the filling of the large basin. (c) Retreat of Akuliarutsip Sermia (CW) on the other hand has reduced boundary length interaction with the lake from 5.1 ± 0.3 km in 1990–95 to 3.9 ± 0.3 km in 2019. (d) Ice-sheet margin just south of Sermeq Kujalleq in 2019. (e) An unnamed glacier to the west of C.H. Ostenfeld Glacier in 2019. (f) Akuliarutsip Sermia in 2019.

[Figure omitted. See PDF]

We find one case where ice-sheet retreat caused a change in boundary type. In 1990, Thrym Glacier, SE Greenland, had a single marine-terminating calving front in Skjoldungen fjord (Fig. 4a). Inland of the calving front, an ice-marginal lake formed where Jomfruen Glacier met Thrym Glacier due to the branching of the fjord to the north and west. However, after retreat of Thrym Glacier into the northwestern branch of the fjord, this lake drained and the Jomfruen Glacier became exposed to the ocean through a narrow channel (Fig. 4b). Retreat of the Thrym Glacier therefore led to the replacement of two ice–lake boundaries (totaling 4.8 ± 0.3 km in length) with a 1.4 ± 0.2 km ice–ocean boundary.

Figure 4.

Retreat of Thrym Glacier, SE Greenland, leads to the replacement of two ice–lake boundaries with an ice–ocean boundary. Location can be found in Figure 1. (a) Thrym Glacier terminates in the ocean in 1990–95 and Jomfruen Glacier, to the north, terminates into a proglacial lake. (b) Jomfruen Glacier terminates into the ocean in 2019 due to the opening of a channel caused by the retreat of Thrym Glacier.

[Figure omitted. See PDF]

4. Discussion

Over the next few 100 years, the Greenland Ice Sheet is expected to shrink and become increasingly grounded above sea level (Aschwanden and others, 2019). Therefore the mass balance of the ice sheet will become increasingly dominated by processes occurring at the ice–atmosphere interface. In contrast, the relative importance of ice–ocean processes on mass balance will wane as marine-terminating glaciers recede into narrower fjords and eventually onto land. Our study demonstrates for the first time that this process is detectable in satellite imagery spanning just 30 years. By combining Landsat remote-sensing imagery with a novel edge detection algorithm, we demonstrate that the overall length of the ice–ocean boundary has declined by 12.3 ± 3.8% between 1990 and 2019. Most regions of the ice sheet (e.g. N, NE, CE, SW, NW) are experiencing reductions in ice–ocean boundary length due to the retreat of marine-terminating glaciers into narrower fjords, rather than their retreat onto land (Figs 2a, b). However there are two regions of the ice sheet where the length of the ice–ocean boundary is increasing (CW and SE). This can partly be explained by the splitting of single marine-terminating glacier calving fronts into multiple calving fronts (Figs 2c, d).

The reconfiguration of the ice–ocean boundary has direct implications for the habitats of pagophilic (ice-dependent) marine species as well as the flux of fresh water, sediment and nutrients into Greenland's oceanic environment (Stuart-Lee and others, 2021). Fjords with ice tongues provide calm, stable environments that concentrate plankton and fish so therefore provide vital foraging habitats for marine mammals (Lomac-MacNair and others, 2018). The loss of stable floating ice tongues that we (e.g. Fig. 2b) and others (Box and Decker, 2011; Münchow and others, 2014; Murray and others, 2015) have observed over the past few decades suggests the potential for complete ice-tongue breakup in the near future. The transition to vertical calving fronts will enhance the production of ice melange (a conglomeration of sea ice, icebergs, bergy bits and brash ice) which may unfavorably influence the habitat use and distribution of some marine mammals (Lomac-MacNair and others, 2018). Having said that, recent research demonstrates that a subpopulation of polar bears in southeast Greenland use ice melange as a platform for hunting during the sea-ice-free season (Laidre and others, 2022), indicating that more ice melange may actually benefit some marine mammals. Likewise, marine-terminating glaciers that retreat onto land will still provide sediment-laden fresh water to the fjord environment. However, there is evidence to suggest that fresh water injected into the fjord at depth has important effects on the ecosystem (Stuart-Lee and others, 2021, 2023). Subsurface plumes entrain large volumes of ambient water, ensuring a continuous resupply of intermediate depth waters, zooplankton and nutrients to the surface where they are more accessible to marine mammals and seabirds (Lydersen and others, 2014). The reduction of subglacial discharge as marine-terminating glaciers retreat onto land will therefore have long-lasting consequences for the complex and productive ecosystems of Greenland.

Rapid reconfiguration of the Greenland Ice Sheet margin is not only occurring at the ice–ocean interface. We find that the length of the ice–lake boundary declined at the ice-sheet scale (−28.2 ± 4.4 km), but that regional trends were more diverging during our study period. In N and NW Greenland, the length of ice–lake boundaries has increased between 1990 and 2019 whereas the length of ice–lake boundaries has declined in the SW and CW sectors. Changes in ice–lake boundary will likely influence rates of mass loss and retreat of the ice-sheet margin. Previous studies in Alaska, the Himalaya, New Zealand and Patagonia have found that mass loss of glaciers in contact with lakes is larger than that of similar-sized land-terminating glaciers and retreat rates are faster (Larsen and others, 2007; Willis and others, 2012; Maurer and others, 2019; Tsutaki and others, 2019; Sutherland and others, 2020). A recent study in western Greenland demonstrated that rates of retreat were also much higher for glaciers that were in contact with lakes than glaciers that terminated on land (Mallalieu and others, 2021). As well as mass loss, these changes are also likely to have lasting impacts on downstream environments and ecosystems. Ice-marginal lakes can impact fjord circulation patterns when they episodically drain (Stuart-Lee and others, 2021). Glacial lake outburst floods also enhance sediment transport rates, which have the potential to alter light attenuation, primary productivity and nutrient availability downstream (Marín and others, 2013; Meerhoff and others, 2019). The continued evolution of ice-marginal lakes therefore has wide ranging consequences on Greenland's terrestrial and coastal environment.

Our findings emphasize that the relative importance of oceanic and freshwater processes on ice-sheet mass change are constantly evolving. The trends we observe during our relatively short study period are, however, unlikely to be linear. Many glaciers (e.g. Humboldt Glacier, Jakobshavn Isbræ and Petermann Glacier) are poised to retreat into wider, deeper channels (Morlighem and others, 2014; Aschwanden and others, 2019; Williams and others, 2021). In doing so, the calving fronts of these glaciers would widen and deepen, which would, at least temporarily, enhance submarine melting and iceberg calving. Accurately forecasting frontal ablation at the ice–ocean interface and the associated retreat of Greenland's marine-terminating glaciers will therefore remain a critical component of sea-level projections in the near future.

Our study extends previous research on Greenland Ice Sheet ice-marginal lakes, although we note that the definition of ‘ice-marginal lakes’ has evolved. For example, Carrivick and Quincey (2014) mapped ice-marginal lakes in western Greenland regardless of how far they were from the ice-sheet margin. How and others (2021) mapped ice-marginal lakes for the entire of Greenland but focused on those that were within <1 km of the ice-sheet margin. Both these studies focused on characterizing changes in the area of these lakes. More recently, Mallalieu and others (2021) defined ice-marginal lakes as those in contact with the ice-sheet margin by intersecting the lakes identified by Carrivick and Quincey (2014) with an ice-sheet perimeter. But while they investigated both lake and ocean boundaries between 1987 and 2015, they focused on a 5000 km length of the western ice-sheet margin. Carrivick and others (2022) extended the approach of Mallalieu and others (2021) to the entire Greenland Ice Sheet by intersecting ice-marginal lakes identified by How and others (2021) with the ice-sheet perimeter. However, they only investigated one time period (2017). We extend Mallalieu and others (2021) and Carrivick and others (2022) by investigating both lake and ocean boundaries for the entire Greenland Ice Sheet for multiple time periods between 1990 and 2019.

Even though our study region and time periods are slightly different, it is possible to make some comparisons between the findings of Mallalieu and others (2021), Carrivick and others (2022) and our study. Mallalieu and others (2021), for example, identified 374 ice–lake boundaries with a total length of 434 km in CW and SW Greenland in 2015: equivalent to 8.8% of their ice-sheet perimeter. Carrivick and others (2022), on the other hand, found more than double the number and length of ice–lake boundaries in CW and SW Greenland in 2017 (n = 1048, total length of 929 km): equivalent to 16% of their ice-sheet perimeter. It is possible that some of these differences can be explained by the fact that Carrivick and others (2022) include ‘inland’ ice-marginal lakes, such as those around nunataks, whereas Mallalieu and others (2021) focus on the outer perimeter of the ice sheet. However, most of the ice-marginal lakes in both datasets are located at the outer margin of the ice sheet. Instead, the differences between the two studies emphasize the challenges with mapping ice-marginal lakes across the Greenland Ice Sheet.

We find 329 ice–lake boundaries in CW and SW Greenland in 2019 with a total length of 262 km: equivalent to 4.8% of the ice-sheet perimeter. Likewise, we find that ~2% of the ice sheet is in contact with ice-marginal lakes in 2019 whereas Carrivick and others (2022) find 10% in 2017. Again, at least some of these differences can be explained by the different treatments of inland ice-marginal lakes around nunataks. But it is also likely that these differences can be explained by the different approaches to mapping ice-marginal lakes. We define ice–lake boundaries as pixels classified as ice that are directly adjacent to pixels classified as water (i.e. lakes). While ice–lake boundaries in Carrivick and others (2022) are lake perimeter points within 200 m of the ice-sheet perimeter (personal communication from Jonathan Carrivick, 2022). Our strict definition of ice–water boundaries makes it more likely that our approach misses some boundaries. Even with 30 m resolution Landsat imagery, it is difficult to determine whether small glaciers (i.e. widths <1 km) are in contact with water. This is made more challenging by the fact that some of the ice-sheet margin, especially in SW Greenland, is debris-covered, and some of the lakes, especially in N Greenland, are ice-covered almost year-round. We mitigate these factors to some extent by conducting rigorous manual quality control on our dataset. Despite this, our findings should be considered a lower bound on the number and length of ice-marginal lakes boundaries. However, given the differences in ice–lake boundary lengths between Carrivick and others (2022) and Mallalieu and others (2021) in CW and SW Greenland, it is also possible that Carrivick and others (2022) overestimate the length of ice sheet in contact with lakes. In the context of the three studies discussed here, a conservative estimate on the total length of ice–lake boundaries in CW and SW Greenland Ice Sheet (combined) at the outer margin of the ice sheet would be 5–9% (in 2019). For the Greenland Ice Sheet as a whole, we expect that the ice–lake boundary is no more than 2–5% of the outer margin of the ice sheet. In other words, the length of ice–lake boundaries is likely smaller than the ice–ocean boundary length which both our study and Carrivick and others (2022) found to be ~5%.

5. Conclusions

In this study, we used three decades of Landsat imagery to determine the length of the Greenland Ice Sheet margin that interacts with both freshwater lakes and the ocean. To do this, we first developed a semi-automated approach to detect boundaries between the ice sheet and water. We then manually checked and edited many of the boundaries to improve the quality of the final dataset. We found that, during our 1990–2019 study period, the length of ice-sheet perimeter in contact with the ocean decreased by 196.2 ± 10.4 km (12.3 ± 3.8%). This was mainly due to the retreat of marine-terminating glaciers into narrower fjords since we did not find much evidence of marine-terminating glaciers retreating on land. Interestingly, the ice–ocean boundary length increased in two regions of the ice sheet (SE and CW) which can be partly attributed to the splitting of marine-terminating glacier calving fronts into multiple calving fronts. The exposure of the ice sheet to ice-marginal lakes also exhibited complex trends during our study period. The length of ice–lake boundaries, for example, increased in SW, N and NW Greenland but declined in SE and CE Greenland. Overall, our large-scale mapping indicates that substantial reconfiguration of the ice-sheet margin boundary has occurred over the past 30 years. These changes will have direct consequences on the rate and character of future ice-sheet retreat. They will also have implications for the habitats of pagophilic species, such as seals and polar bears, as well as the flux of fresh water, sediment and nutrients into Greenland's oceanic environment.