Introduction

In the past 50 years, mass loss from the Greenland Ice Sheet (GIS) has increased global sea level by 13.8 mm, with roughly half of this increase occurring within a single decade (2009-2018)1,2. Southeast Greenland is the second largest regional contributor to sea level rise (SLR) and has shown some of the most significant changes with mass loss increasing from 41 ± 17 Gt/yr to 286 ± 20 Gt/yr during the last 30 years. This has primarily been driven by increased melt and dynamic instabilities of marine-terminating glaciers2,3,4,5,6,7. With projected increases in Greenland’s mass loss expected to result in part from increased dynamic mass loss, accurately modelling future glacier dynamics is important for minimising the projection uncertainties that are correlated with the complexity of glacier behaviour8. Accurately capturing and reconstructing intricate and long-term glacial histories, through historical analogue data and archives9, is critical for determining if recent dynamic changes and/or collapses are directly linked to current warming or are the result of past ice flow reorganisations. For example, studies have highlighted that the current dynamic changes observed in the Amundsen Sea sector of Antarctica are the result of a long-term process that began in the 1940s10. It is therefore essential to capture long-term records of glacier evolution in order to produce robust numerical simulations that are capable of reproducing different glacier configurations, to calibrate critical ice flow parameters, and for model validation11.

One significant complexity not yet well represented in current models is related to dynamic instabilities and thinning, that lead to ice flow reconfigurations on the inter-catchment scale. This phenomenon, known as ice flow piracy, involves the redirection of ice flow from one glacier’s catchment basin to another12 and may significantly influence regional surface elevation changes and subsequent retreat rates. Glaciers that lose large parts of their drainage basins will experience reduced ice supply, impacting their dynamics and potentially contributing to changes in ice mass and retreat, while those that gain ice flow will experience increased flow and advance, both independently of the prevailing climatic conditions. Glacial piracy has been sparingly documented in Antarctica13,14, but it is expected to occur in complex glacial catchments elsewhere. Although it may play a critical role in the dynamics of Greenland’s ice sheets and glaciers, ice flow piracy has not been extensively documented here, and no comprehensive observation of this complex process has yet been reported.

Understanding the role of ice flow piracy is crucial for knowing how ice streams in complex regional catchments will respond to a changing climate, which is essential for estimating current and future mass loss. For example, changes in drainage basins resulting from ice flow piracy can impact the partitioning of mass balance into catchments in the input-output approaches, since this partitioning is based on discharge at the terminus and catchment-scale surface mass balance (SMB) estimates that are fixed in time1,15. Additionally, reference SMB, often averaged over the period 1970s–1990s, is sometimes used to scale ice discharge when ice velocity data are unavailable1,16. Hence, inaccurately defined catchments can directly impact these estimations and further propagate uncertainties. In terms of modelling, regional studies simulate dynamic evolution using catchments that are also fixed in time11,17,18,19,20,21,22, which can incorrectly specify the amount of ice influx, affecting the calculated SMB and ice discharge. Finally, understanding how ice flow piracy can occur is critical for understanding the physics of ice flow, as it is closely linked to changes in geometry, ice-bed interactions, thermal regimes, and subglacial hydrology23,24,25,26,27. The rearrangement of basins caused by shifts in glacier flow will alter glacier discharge and calving around Greenland, directly impacting the mass balance23 of individual glaciers. Specifically, accurately capturing the spatial distribution of ice discharge across the ice sheet is crucial, as it directly influences freshwater flux into the ocean, which can significantly impact ocean circulation28,29.

The Midgaard Glacier (MG), located in Sermilik Fjord (Fig. 1), Southeast Greenland (SE), has displayed some of the most remarkable mass loss observed in Greenland in the 20th and 21st centuries with a total retreat of 35.5 km since the end of the Little Ice Age (LIA – defined locally as 1900)5,6,30,31. Although covering less than 2% of Southeast Greenland’s glacierised area, the catchment is ranked the 5th largest contributor to mass loss in Greenland1. Since 1972, MG catchment has lost 42% of its total ice1, which accounts for 13% of Southeast Greenland’s mass loss and 3% of GIS’s mass loss between 1972 and 20181. On its own, MG catchment holds a potential of 0.8 mm Sea Level Equivalent (SLE)1; however, it is located in a topographically constrained and complex glaciological system (Fig. 1b) characterised by numerous tributaries that connect to neighbouring outlets. Thus, in the context of ice flow piracy, it is critical to consider the impact of neighbouring glacier catchments. Therefore, when including the neighbouring Fenris Glacier, Glacier de France, Knud Rasmussen Glacier, and Kiv Steenstrup Nodre Brae catchments (Supplementary Fig. 1), the combined potential SLR reaches up to 19.2 mm SLE1.

Fig. 1: Context of the study area and retreat in the region.
figure 1

a Inset showing the location of the study area in Greenland. b Map of the study area showing Sermilik Fjord and the location of the major glaciers in the study: Midgaard Glacier (MG), Franche Comté Glacier (FCG), Nigertiip Apusiia (NA), Kattilersopia (KA), and Pourquoi-Pas (PP). Hillshade is produced from the ArcticDEM64 mosaic. c Terminus positions from 1933 to 2022 based on historical analogue images and modern satellite imagery for MG. d Evolution of the medial moraine at the intersection between KA-PP-NA from 1933 to 2022, also based on historical and modern satellite imagery. e Terminus positions from 1933 to 2022 based on historical analogue images and modern satellite imagery for KA. (Background image c-e: LANDSAT 9 2022) f Relative glacier front position of MG (blue) and KA (orange) with the year 2011 highlighted in yellow to mark the end the single calving front named MG and the formation of two new fronts now called FCG (green) and NA (red).

Full size image

In this study, we present the longest record of glacier dynamic evolution in the region spanning a 90-year period from 1933 to 2023. We analyse a comprehensive set of data, including historical observational analogue aerial and satellite images, along with geomorphological evidence, historical maps, and modern remotely sensed data to shed light on the complex changes that have occurred in this region. By examining dynamic fluctuations, surface flow velocities, and successive flow reorganisations within the catchment that have occurred from the end of the LIA to the present, we were able to reconstruct the flow configuration evolution and examine factors influencing ice flow in this dynamically complex region.

Results

In this study, we employ a multi-data approach to document the collapse of Midgaard Glacier, and construct timeseries of resulting dynamic changes in its catchment as well as in the neighbouring Glacier de France catchment (Supplementary Fig. 1). The glaciers we investigate in the Midgaard Glacier catchment include Midgaard Glacier (MG), Pourquoi-Pas (PP), Franche Comté Glacier (FCG), Nigertiip Apusiia (NA), and Femstjernen glaciers. In the Glacier de France (GDF) catchment, we look mainly at the upper and lower reaches of Kattilersopia glacier (KA) (Fig. 1b). Note that in older studies, ‘Glacier de France’ referred to both the catchment and the glacier. However, to standardize terminology and use the official name32, ‘Kattilersopia’ will now be used for the glacier, and ‘Glacier de France’ for the catchment. Below, we first describe the observed changes on MG in terms of retreat, thinning, and flow, after which we lay out the effects of these changes on the neighbouring glaciers in the study area.

Collapse of Midgaard Glacier

Based on manually digitised glacier front data (Fig. 1c; see Methods), we see that the retreat of MG’s calving front was marked by two phases of rapid retreat (Fig. 1f). The first phase occurred between 1933 and 1979 with a total retreat of 19 km. Bathymetric data from multibeam echo sounding reveal that during that time period, MG was dislodged from its LIA position on a shallow sill 200 m below sea level33 (Supplementary Fig. 1). The retreat rate decreased from 413 ± 22 m/yr before 1980 to 170 ± 16 m/yr until 2002. The second rapid retreat phase occurred between 2003–2011 with a retreat of 9 km, during which time the retreat rate peaked at 4021 ± 28 m/yr in 2008–2009. Overall, the total cumulative retreat of MG was 32.5 km over the entire study period. After 2011, MG had retreated into the FCG fjord in the west and NA fjord in the east producing two separate calving fronts (Fig. 1c), marking the end of the glacier known as MG. Between 2012 and 2021 NA retreated in total 1.5 km, at ~129 ± 4 m/yr and advanced 645 ± 28 m/yr between 2021–2022, whereas FCG retreated in total 1.1 km, at ~65 ± 3 m/yr between 2012–2022 (Fig. 1f).

Effects of Midgaard Glacier collapse on its catchment

Utilising historical aerial images and analysing the difference between derived Digital Elevation Models (DEMs) (Fig. 2 and Supplementary Fig. 2), we measure thinning rates at location 1, on Lower FCG (Supplementary Fig. 3) of 2.5 ± 0.5 m/yr from 1933 to 1981 (see Methods), which then rose to an average rate of 15.5 ± 2.5 m/yr in 2011–2014, peaking at 26.7 ± 5.0 m/yr in 2011, preceding the complete disappearance of MG. Later on at this location, thinning rates decreased to 4.3 ± 1.8 m/yr until 2021. According to historical DEM, MG could have thinned at most 746 ± 108 m during the entire study period between 1933 and 2021 (Supplementary Figs. 2c, 4b).

Fig. 2: Four decades of thinning and velocity data reveal remarkable responses to MG collapse.
figure 2

Elevation changes of the Midgaard Glacier (MG) and Glacier De France (GDF) catchments between 1981–2015 (a) and the MG catchment between 2011–2015 and the GDF between 2008–2013 (b). The glaciers in the MG and GDF catchments, and most notably the FCG, have experienced thinning and retreat since 1981, but due to the complex interplay of the streams, some areas have experienced slight thickening. Time series for the areas (average over 1kmx1km) marked by coloured circles are displayed in (c-h). Hillshades are produced from the ArcticDEM64 mosaic. c–h Time series showing ice speed from 1986 to 2023 (c–e), with shading indicating the uncertainty estimates by the method60 and elevation changes over different periods from 1981 to 2021 (f–h), the error bars in the plots represent the RSS (Root Sum of Squares) of the NMAD (Normalized Median Absolute Deviation) of each DEM involved. Time series are categorised by glacier branch, each with a distinct y-axis for enhanced clarity.

Full size image

Since 1980, MG and its contributing glaciers have also shown significant changes in velocity (Fig. 2c–e). Using velocity data derived from satellite sensors (see Methods), the speed at location 1 experienced an increase from an average of 1625 ± 60 m/yr between 1986–2002 to 2789 ± 45 m/yr between 2003–2013 (Fig. 2c). Then, the velocity decreased between 2014–2023, down to 2154 ± 54 m/yr in 2023, hence almost back to its pre-2000s conditions. During that last decade, we find the magnitude of the seasonal velocity fluctuation was 20% of the annual peak velocity.

Combining ice surface reconstruction from LIA trimlines5 and satellite mass balance calculations1,5, we estimate an average mass loss rate of 2 ± 0.6 Gt/yr in 1900–2003 for the MG catchment. Mass loss then increased by 85% to 3.7 ± 0.7 Gt/yr until 20191, with 80% of this loss attributable to ice discharge (Supplementary Fig. 5). In total, the MG catchment lost 247 ± 66 Gt between the end of the LIA to 2018, resulting in the equivalent of 0.7 mm global mean SLR.

Catchment reconfiguration and ice flow piracy

Between the end of LIA and 1933, historical mapping shows that the Femstjernen glacier (Fig. 1) flowed through both FCG and KA, crossing PP glacier. Based on visible surface lakes (Supplementary Fig. 2f) that form during periods of low flow, the NA glacier experienced minimal flow dynamics over this period and had a catchment of 152 km2 (only 4.5% of MG’s catchment) and was feeding into the main trunk of MG.

The significant post-1933 retreat and thinning of MG and thinning of FCG led to the first flow reversal that we observed in this area. This occurred sometime before the end of the 1970s, as deduced from the presence of a medial moraine in 1981 aerial imagery (Fig. 3). Due to this thinning and retreat, the ice from Femstjernen to KA transitioned to flow solely toward FCG, hence reversing flow and establishing a new ice divide (Figs. 2, 3). Between 1981 and 2021, the divide progressively migrated south-eastwards towards the KA fjord at an average speed of 133 ± 93 m/yr (Fig. 3). This first ice flow piracy event had a significant impact on the catchments’ areas. While the MG outlet gathered 62% of the total glacierised area in the region in 1933, this expanded to 77.5% by 1981 with the addition of NA (Fig. 4). During the same period the basin area of GDF decreased by 41%.

Fig. 3: Geomorphological evidence of ice piracy.
figure 3

a Overview showing the location of b, c, and d. b Evidence of past ice flow. The 1981 aerial image shows a medial moraine in the red box, indicating historic ice flow from Femstjernen to Kattilersopia (KA). The upper red arrow marks the probable origin of the medial moraine. c Evidence of ice flow redirection. The medial moraine loop indicates ice from Pourquoi-Pas (PP) is now directed towards Nigertiip Apusiia (NA), whereas it previously flowed towards KA. d Migration of the ice divide between Franche Comté Glacier (FCG) and KA. Visualization of flow reconfiguration and ice divide migration at the intersection between FCG and KA. Points indicate manual tracking of the ice divide along five profiles, showing spatial retreat over different years. The white dashed line marks the profile section shown in (e). e Elevation profile (m) along the white dashed line in (d), from A to B.

Full size image
Fig. 4: Evolution of ice flow and catchment configuration over the study period.
figure 4

a–c The pink arrows indicate the direction of ice flow while the red arrows indicate changes compared to the previous period. d–f Distribution of basins and their contributions to the Midgaard Glacier outlet or Kattilersopia outlet. Percentages indicate the catchments’ relative share of the total glacierised catchment at each timestep.

Full size image

Tracking of medial moraines revealed another significant ice flow piracy event at the intersection between KA and PP. Around 1990–1995, PP reduced its flow significantly towards lower KA, redirecting the flow almost entirely into NA. This change separated upper KA from its lower counterpart. Following the U-shaped medial moraine through time reveals a 7.7 km migration at an average rate of 240 m/yr between 1990 and 2022 (Fig. 1d). The most rapid flow migration of the medial moraine took place during 2005–2010, at a rate of 400 m/yr and covering 2 km. This is consistent with changes in surface elevations at location 4, which switched from thinning at 2.9 ± 0.5 m/yr in 1981–2008 to thickening at 1.8 ± 0.7 m/yr in 2011–2017 (Fig. 2b, g).

Satellite derived ice flow velocity measurements at NA (location 4) are also consistent with this ice flow piracy event. The ice velocity was stable from 1990 to 2003, with an average of 142 ± 61 m/yr (Fig.2d). Then, the velocity accelerated to 1076 ± 48 m/yr in 2009, followed by a decline between 2010–2021 to 161 ± 29 m/yr (Fig. 2d). We also measure remarkable acceleration of the entire NA branch after 2021, where the velocity in the vicinity of the ice front at location 3 increased from an average of 212 ± 18 m/yr in September 2021 up to 4642 ± 103 m/yr in May 2022 (Supplementary Figs. 3 and 6), before going back to its previous rate in July of the same year. This last flow reversal resulted in an additional expansion of the MG catchment by 87% of the total ice area. The flow switches and subsequent reorganisations has thus resulted in a net increase of 18% throughout the entire study period, despite significant retreat (Fig. 4 and Supplementary Fig. 7). At the same time, GDF catchment net area decreased by 66% between 1933 and 2022.

The isolation of the lower section of KA as a consequence of two ice flow piracy events impacted the dynamic of its lower section. The KA ice front remained stable until 2002 and then started to retreat afterward, for a cumulative recession of 5.1 km between 2003–2022 (Fig. 1e,f). Surface elevation exhibited minor thinning of 2.1 ± 0.9 m/yr in 1981–2007 at location 7, which increased by 570% between 2008–2021, peaking at 32.7 ± 3.7 m/yr in 2018 (Fig. 2h). In terms of surface flow velocity, the lower section of KA (location 7) remained stable until 2010, with an ice front flow at 511 ± 29 m/yr (Fig. 2e). The velocity then quadrupled up to 2201 ± 34 m/yr in 2019 before decreasing back to 1137 ± 95 m/yr in 2023. This is correlated with significant dynamical changes for the two northern tributaries of lower KA (location 8 and 9; Supplementary Figs. 3, 6), which doubled their ice velocity between 1995–2000 and 2015–2020, hence supplying the lower section of KA with more ice. We posit that the enhanced thinning of KA permitted changes in ice dynamics of glacier tributaries, which are not observed close to the intersection between PP and KA where thinning is less pronounced. Finally, we are not able to measure any trends in speeds on the upper KA, where ice flow is often lower than 20 m/yr.

Discussion

Before the 20th century, during the LIA34, when surface air temperatures were relatively cold, Greenland’s glaciers grew considerably34, including MG which reached a 200-metre-deep fjord sill where the glacier stabilised until the 1930s33 (Fig. 1). Its stability at this shallow depth might have delayed its retreat compared to other glaciers in the region35,36. However, due to rising ocean and air temperatures (Supplementary Fig. 8f,g) between the 1920s–1940s, MG was dislodged from the sill after the 1930s, which was the starting point of a complex cascade of interconnected events. Between the 1930s and the 1980s, the rate of retreat of MG increased rapidly, which is consistent with an unstable retreat on a retrograde slope that reaches depths of 500 m37,38. The warmer climatic conditions during that period must have favoured the intrusion of warm Atlantic waters, originating from the Irminger current38,39. Recent CTD measurements have shown that warm Atlantic water (>3 °C at depth) could be found down to 200 metres deep in Southeast Greenland38. Ocean conditions likely played a key role in the chain of events40,41,42,43. Two multidecadal periods of warmer ocean and atmosphere conditions in the vicinity of Greenland stand out in the 20th century; the 1920s–1960s and the post-2000s separated by relatively colder sea surface temperatures in the 1960s–1990s (Supplementary Fig. 8f). The post-2000s warming is prominent and is also reflected in a regime shift in the sea ice conditions on the SE shelf44.

Both warm periods in the 1920s–1960s and the post-2000s have likely led to enhanced submarine melting, rapid ice front retreat, increased ice velocity and dynamic thinning from MG and subsequently FCG glacier (Fig. 2c and Supplementary Fig. 4). This intensive thinning isolated the upper section of KA, located at higher elevation, cutting off its ice supply from Femstjernen (Fig. 4), which triggered the first ice flow reconfiguration of the region in the early 70 s (Fig. 4). This disruption in ice inflow at upper KA led to significant ice thinning at higher elevations (Supplementary Fig. 4c, d), initiating a series of subsequent dynamic changes at lower elevations (Fig. 4) after the early 1970s. The lack of ice supply to KA upper section, from upstream ice flow of Femstjernen (Fig. 4) along with increased thinning (Supplementary Fig. 2 and Supplementary Fig. 4) has weakened the ice tongue close to the PP intersection. This allowed the flow of PP through KA, starting from the early 1990s and likely caused NA to move faster and thicken after 2004 (Fig. 2 and Supplementary Fig. 3). This change in flow pattern also likely initiated the retreat at the KA front from 2003 (Fig.1f), as the ice supply at upper KAinto the lower part of KA had been significantly reduced. Despite the reduction in ice flow velocity downstream of the PP-NA-KA intersection, KA’s front ice flow velocity remained constant until 2010 and even quadrupled from 2010 to 2020 (Fig. 2e), likely partly due to ice supply coming from the concomitant acceleration of the Northern tributaries of lower KA (Supplementary Fig. 3d) and accelerated undercutting. However, the increase in warm ocean temperatures39 post-2000s may have also increased ice front undercutting, which could have intensified the retreat of the front of KA observed in 2003–2022. Therefore, the reorganisation of flow isolated some basins, reducing glacier ice input and accelerating their retreat. This, combined with previously mentioned factors, likely played a crucial role in the retreat acceleration of KA. Indeed, more than half of KA’s total observed retreat occurred between 2015 and 2022, and its discharge doubled from 2009 to 20181. At the NA glacier front (location 3), it is important to note that the glacier shrank by 3.20 ± 0.20 km2 due to a front retreat that occurred between 2012 and 2021. However, from 2021 to 2022, the front advanced by 670 ± 28 m (Fig. 1c,d) leading to an increase in the glacier surface area by 1.03 ± 0.01 km2, which is likely related to the widespread increase in speed of this glacier after 2021. Specifically, the surge-type behaviour measured in May 2022 is a strong candidate to explain the frontal advance of NA.

The initial retreat since the 1930s and the intensified events post-2000s align with warmer atmosphere and ocean temperatures, as illustrated in Supplementary Fig. 8 summarising key events in MG and GDF catchments. However, assessing the synchronicity of all changes with climate warming is challenging, as these changes trigger further cascading events, highlighting the intricate dependency and prolonged reaction times within the system. Also, directly correlating climate with glacier events is difficult due to the system’s complexity and limited understanding and data availability of the various physical processes involved. The interplay between climate and glacier response is not fully understood, as events are intertwined between climate and glacier response. Thus, current warming alone cannot fully explain recent changes in glacier dynamics.

Although our results suggest that changes in flow are initially driven by significant variations in glacier geometries, such as thinning, other studies indicate that factors like substrate deformability and subglacial hydrology also play a role in altering ice flow patterns23. Gaining a better understanding of geology, subglacial hydrology and basal motion would improve our ability to link climate drivers with observed glacier dynamics23. Nonetheless, this study has demonstrated that current rapid ice reconfiguration in glacial outlet dynamics can be attributed to previous ice flow piracy in flow and long-term changes in the overall dynamic regime. Given the current projected increase in melt rates8, which will further alter glacier geometry, through thinning, as well as subglacial hydrology45,46, are factors that are likely to contribute to future potential reconfigurations of ice sheet flow.

The case of MG is a prime example of such a cascade of events, where post-LIA retreat triggered substantial surface lowering, allowing major catchment piracy throughout the following nine decades. By examining the history of these dynamic changes, we can directly measure their impact on ice flow and the retreat of glacier fronts, such as KA and NA glaciers. Capturing these complex changes in models is essential for accurately replicating and projecting future ice sheet dynamics.

This study highlights that future glacial dynamics, especially in complex terrain, are controlled by the interplay between thinning and the surrounding terrain. Studying mountainous regions, as in our study, is highly complex due to intricate topography and poorly defined bedrock constraints, creating significant uncertainty. If bed topography is complex, a combination of extreme thinning and strong glacier dynamics can lead to ice flow piracy, resulting in drastic changes in ice dynamics for neighbouring outlets. Currently, thinning rates can vary significantly among tributaries in SE Greenland (e.g., Supplementary Fig. 3), implying the potential for future switches that would further alter local ice flow patterns. Local observations of flow changes are scarce, and the phenomenon is poorly understood.

The dataset presented here is a comprehensive collection on this rarely observed phenomenon, providing a stepping stone for improved modelling. To explain these changes, it is therefore essential to assimilate and reproduce in models these observations, in order to understand the underlying physics. Currently, most ice-sheet-wide models typically have a resolution >1 km at the ice front8,47, which is too coarse to accurately replicate dynamic changes observed in this study, where variations can occur over a few hundred metres. Higher resolution simulations are therefore needed to better replicate ice flow rearrangement at finer scales23,48 and the resulting mass loss evolution through discharge and calving, as even minor local changes can have significant regional effects, influencing the overall mass balance23,24,27,49.

More broadly, our study also shows the importance of closely monitoring glacier changes in glacier flow patterns at ice divides, particularly in the context of ice coring, where it is usually assumed that the flow at ice divides is close to zero23. We have shown that understanding ice flow piracy is crucial for accurately predicting and modelling glacial behaviour, particularly in the context of climate change and its impact on Greenland’s ice sheet. Local changes can have a regional impact by reconfiguring basins, which directly affects glacier discharge. Our work sheds light on previously unknown histories of glacial complexity and emphasises the significance of knowing the history and causes of long-term ice dynamics when modelling and interpreting contemporary glacial changes. Therefore, expanding this analysis at the ice sheet scale should be performed and continued in the future.

Methods

Historical aerial images and historical map—7th Thule expedition Knud Rasmussen 1933

The historical map from Knud Rasmussen’s 7th Thule expedition50 in 1933 was georeferenced for digitising the ice, trimlines, glacier fronts, glacier flow direction, and glacier moraines. However, in recent decades, a vast collection of the original aerial photographs has been scanned using a photogrammetry-grade scanner30 and has been utilised in various scientific studies5,30.

Processing of the historical images from the 7th Thule expedition Knud Rasmussen 1933

The oblique images have a focal length of 127 mm and were processed in Agisoft Metashape v. 1.7.4, using structure-from-motion (SfM) photogrammetry. As the images are of varying quality and with inconsistent overlap, we used a selection of the best images covering Midgaard Glacier. First, we standardized the internal geometry of the images by manually detecting the image fiducials and aligning them by extracting up to 40,000 tie points per image. The 3D model was georeferenced by manually placing 23 ground control points (GCPs) in the aerial images. The GCPs used are stable bedrock features with specified (x,y,z) locations extracted from the Arctic DEM51. We performed a dense multi-view stereo (MVS) reconstruction and grid the point cloud into DEMs and lastly, the orthophoto mosaic was produced. We co-registered our 1933 DEM to the Arctic DEM and assessed its accuracy relative to the Arctic DEM, by comparing elevations over stable bedrock. We found our 1933 DEM to have an accuracy of 108 m, which is likely due to a combination of varying image quality, insufficient overlap, and unknown camera parameters. Therefore, the 1933 DEM is used as a rough approximate measure of elevation change since 1933 (Supplementary Fig. 2).

Historical satellite imagery: 1960s declassified Corona photographs and Sojuzkarta KFA-1000 images

We used declassified Corona satellite images from 1962, 1964, and 1966, at 50, 40, and 5 m resolution, respectively. We used the KFA-1000 images52 from 1992, digitised at a resolution of 0.012 mm, providing a ground resolution of 3  m52. We used them to digitise the extent of ice, lateral elevation, glacier fronts, and glacier moraines.

To georeference these images, we employed ArcMap 10.8. In this study, we simply identified GCPs by using visible bedrock with unique geomorphological features also visible on a georeferenced 25 m orthomosaic from 198130 to define (x, y) coordinates. We identified 617 GPS points in four scenes in total representing 1234 manual placements, and we concentrated most GCPS around the glacierised areas so as to constrain the ice surface reconstruction in these areas. We then performed a spline transformation with an average residual error of less than a metre for all georeferencing.

Tracking calving front and medial moraine positions

We used historical aerial and satellite images, along with modern digital imagery, spanning from 1933 to 2022, to map the calving fronts of MG and KA glaciers. This almost 90-year record time series was generated by combining manual tracing from the historical dataset with selected modern images (Supplementary Table 1) and completed with glacier fronts from the Termpicks dataset26. Our image selection focused on a specific period between July 1st and October 31st, corresponding to the summer season, to maintain consistency. The ice front positions of MG and KA glaciers, as well as the medial moraine separating KA and PP, were manually digitised using QGIS. In total 171 frontal positions were digitised and 43 medial moraine locations. Relative errors of digitisation (Supplementary Table 1) were estimated by following the methods from Bevan et al.53. However, we encountered limitations during our analysis. Some images taken between August and September had snow cover and could not be digitised. Furthermore, the Landsat-7 satellite’s scan line corrector failure posed another challenge54.

We determined glacier retreat using the centreline method, employing the intersections between digitised fronts and topographic centreline profiles as point indicators for ice front positions55. This approach is less labour-intensive, allowing measurements even if part of the front is obscured. While this method captures variability at a single point, the results were sufficient for the objectives of our study56 as we are interested in large-scale changes.

DEM post-processing

We used DEMs sourced from various datasets (Supplementary Table 1) spanning the period from 1981 to 2021, specifically constrained between 1st July and 31st October whenever possible. To ensure accuracy and consistency, we co-registered57,58 the DEMs and each individual ArcticDEM strip59, using a regular grid with a resolution of 40 × 40 m, to the Worldview DEM ArcticDEM Mosaic V351 using Shean et al.58 In addition to the DEMs, we incorporated altimetry data (Supplementary Table 1) and compared the differences between each DEM and the corresponding altimetry data over stable bedrock, with a mean accuracy of 2.7 m. To assess the absolute model accuracy for each DEM, we calculated the Normalized Median Absolute Deviation (NMAD) of the elevation differences over stable bedrock59. To evaluate changes in ice dynamics since 1981, we computed elevation changes (dh/dt) for all possible pairs and extracted elevation profiles along the centre lines of MG, KA, and PP-NA, different close-up transects (Supplementary Fig. 4), as well as averaging thinning rate over a 1 km × 1 km square at different key locations (Fig. 2 and Supplementary Fig. 3). For each pair we determined the uncertainty as the RSS (Root Sum of Squares) of their NMADs.

Ice velocity mapping

To generate surface flow velocity maps, we use data from four satellite sensors between 2013 and 2023, we analysed surface displacements to monitor glacier velocity. Three of these sensors employ optical imaging: ESA’s Sentinel-2 (S2) and NASA’s Landsat-7 (L7) and Landsat-8 (L8), and one utilises synthetic aperture radar technology: ESA’s S1. By leveraging persistent surface features or speckles, we mapped ice displacements between consecutive images. To calculate the normalised cross-correlations between primary and secondary image chips, we used repeat cycles shorter than 30 days for L7/L8 and S2, and 12 days for S160,61,62. We calibrated our displacement maps by referencing ice velocity products from previous surveys61. The calibrated maps were resampled to a 150-metre posting in the north-polar stereographic projection (EPSG:3413). To complement the time series, we included historical measurements from various sources such as ERS-1/2, RADARSAT-1, ALOS PALSAR, Envisat ASAR, Landsat-4 to 7, and TerraSAR-X.

Drainage basins

Drainage basins and sub-basins were generated using flow accumulation generated in ArcGIS 10.3 from the ArcticDEM51 ice surface mosaic. We relied on supraglacial catchments as the best available bed elevation model is largely interpolated and does not provide the necessary spatial scale63. The analysis defined supraglacial topographic divides for the entire area and assigns incoming flow to each pixel, enabling the determination of a catchment for each pixel in the DEM. Four points were placed at all major ice intersections and outlets with high aggregated flow. This method divided the major catchment into several sub-catchments (see Supplementary Fig. 1), reflecting individual components of the modern flow organisation and historical catchment switches. These sub-catchments were then manually assigned to the correct outlets for past flow configurations.