A Greenland-wide empirical reconstruction of paleo ice sheet retreat informed by ice extent markers: PaleoGrIS version 1.0

The Greenland Ice Sheet is a large contributor to global sea level rise, and current mass losses are projected to accelerate. However, model projections of future ice sheet evolution are limited by the fact that the ice sheet is not in equilibrium with present-day climate but is still adjusting to past changes that occurred over thousands of years. While the influence of such committed adjustments on future ice sheet evolution remains unquantified, it could be addressed by calibrating numerical ice sheet models over larger timescales and, importantly, against empirical data on ice margin positions. To enable such paleo data–model interactions, we need Greenland-wide empirical reconstructions of past ice sheet extent that combine geomorphological and geochronological evidence. Despite an increasing number of field studies producing new chronologies, such a reconstruction is currently lacking in Greenland. Furthermore, a time slice reconstruction can help to (i) answer open questions regarding the rate and pattern of ice margin evolution in Greenland since the glacial maximum, (ii) develop a standardised record of empirical data, and (iii) identify new sites for future field campaigns. Based on these motivations, we here present PaleoGrIS 1.0, a new Greenland-wide isochrone reconstruction of ice sheet extent evolution through the Late Glacial and early- to mid-Holocene informed by both geomorphological and geochronological markers. Our isochrones have a temporal resolution of 500 years and span ∼ 7.5 kyr from approximately 14 to 6.5 kyr BP. We describe the resulting reconstruction of the shrinking ice sheet and conduct a series of ice-sheet-wide and regional analyses to quantify retreat rates, areal extent change, and their variability across space and time. During the Late Glacial and early- to mid-Holocene, we find the Greenland Ice Sheet has lost about one-third of its areal extent (0.89 million km²). Between ∼ 14 and ∼ 8.5 kyr BP, it experienced a near-constant rate of areal extent loss of 170 ± 27 km² yr⁻¹. We find that the ice-sheet-scale pattern of margin retreat is well correlated to atmospheric and oceanic temperature variations, which implies a high sensitivity of the ice sheet to deglacial warming. However, during the Holocene, we observe inertia in the ice sheet system that likely caused a centennial- to millennial-scale time lag in ice extent response. At the regional scale, we observe highly heterogeneous deglacial responses in ice extent evident in both the magnitude and rate of retreat. We hypothesise that non-climatic factors, such as the asymmetrical nature of continental shelves and onshore bed topographies, play important roles in determining the regional- to valley-scale dynamics. PaleoGrIS 1.0 is an open-access database designed to be used by both the empirical and numerical modelling communities. It should prove a useful basis for improved future versions of the reconstruction when new geomorphological and geochronological data become available.

1 Introduction

The Greenland Ice Sheet holds an estimated ∼ 7.4 m of global sea level equivalent (Morlighem et al., 2017). It is currently experiencing mass loss at a rate of ∼ 150 Gt yr⁻¹ (over the 1992–2018 period; The IMBIE Team, 2019; Otosaka et al., 2023), making it a large contributor to current global mean sea level rise. In a warming world, these losses are projected to accelerate (Meredith et al., 2019), with the latest Greenland Ice Sheet Model Intercomparison Project (ISMIP6) estimating a contribution to global mean sea level rise of 32 ± 17 and 90 ± 50 mm by the year 2100 for the Representative Concentration Pathway (RCP)2.6 and RCP8.5 greenhouse gas concentrations scenarios, respectively (Goelzer et al., 2020). The robustness of such projections is limited by the short time span of instrumental evidence (tens to hundreds of years) used for verification and testing (Larsen et al., 2015). This is problematic because ongoing and future ice sheet mass losses are not only exclusively a consequence of contemporary and future climate but also represent adjustments to climate variations stretching further back in time (hundreds to thousands of years; Rogozhina et al., 2011; Yang et al., 2022). As few models include the committed response to Late Quaternary environmental changes, current projection simulations might effectively be starting from a possibly unrealistic state, whose influence is thus far unquantified. Consequently, improved empirical reconstructions of Late Quaternary ice sheet evolution may be useful for refining future predictions of mass loss and global sea level rise by enabling the calibration of ice sheet models over larger timescales encompassing thousands of years (e.g. Lecavalier et al., 2014; Albrecht et al., 2020).

The response of the ice sheet to climate change during the Late Quaternary likely contributes to major imbalances in its current state (Calov and Hutter, 1996), yet many key questions pertaining to this time remain unresolved. The position, rate, and pattern of ice margin retreat during the Late Glacial and early- to mid-Holocene periods (between ∼ 14 and ∼ 6 kyr BP) remains heterogeneously constrained in numerous regions. During that time, the ice sheet was responding to rapid and high-amplitude fluctuations in atmospheric and oceanic boundary conditions (e.g. the Bølling–Allerød interstadial; Buizert et al., 2018b), while also adjusting to fast rates of relative sea level change, following the demise of other ice sheets. Consequent changes in ice sheet margin, volume, and surface mass balance were likely associated with modifications to the configurations of ice divides and to the position and discharge regimes of major ice streams (e.g. Franke et al., 2022). A key question concerns how far the ice sheet retreated behind its contemporary margin in response to the warming of the Holocene thermal maximum (10–5 kyr BP; Weidick et al., 2004; Cartapanis et al., 2022), a possible analogue to atmospheric warming expected in future decades (Funder et al., 2011). Calibrating model simulations against known Holocene margin positions around the full ice sheet perimeter might help answer this latter question.

The Greenland Ice Sheet has been subject to an increasing number of field studies reconstructing the deglacial evolution of the former ice sheet margin, both on land and offshore. These typically use geomorphological, sedimentological, and geochronological dating analyses (e.g. Ten Brink and Weidick, 1974; Hjort, 1979; England, 1985; Bennike and Weidick, 2001; Hughes et al., 2012; Briner et al., 2014; Larsen et al., 2014; Young et al., 2021; Garcia-Oteyza et al., 2022; Sbarra et al., 2022). This growing library of geological evidence enables an improved understanding of the ice sheet response to deglacial climate change at both the regional and continental scale (Sinclair et al., 2016). However, to our knowledge, it appears that a Greenland-Ice-Sheet-wide reconstruction with an open-access, reproducible database of geomorphological and geochronological evidence is lacking. Such a product would prove useful to both the modelling and empirical glaciology communities. It would enable ice-sheet-scale analyses, help keep a standardised record of dating investigations, facilitate targeting new study regions for future fieldwork, and provide a reconstruction for calibrating ice sheet model simulations. We note that ice-sheet-wide summary datasets of accumulated evidence have played an essential and enabling role in advancing the understanding of other paleo ice sheets and in answering key glaciological questions on former ice dynamics (e.g. Dyke, 2004; Bentley et al., 2014; Hughes et al., 2016; Davies et al., 2020; Dalton et al., 2020; Clark et al., 2022).

Based on motivations outlined above, we present an ice-sheet-wide isochrone reconstruction of the Greenland Ice Sheet's Late Glacial and early- to mid-Holocene evolution informed by an ice-sheet-scale dataset of geomorphological and geochronological ice extent markers compiled in an open-access database. This reconstruction, hereafter referred to as PaleoGrIS 1.0, has a temporal resolution of 500 years, and spans 7500 years, between ∼ 14 and ∼ 6.5 kyr BP, the period bracketing most terrestrial landform evidence of ice margin change. This reconstruction represents a first attempt at combining both geochronological and geomorphological records to estimate the pattern and timing of the ice sheet deglacial retreat. In this work, we describe our methodology in detail to make the reconstruction transparent and reproducible, with the hope that improved versions of PaleoGrIS can be produced as more empirical evidence arises. We use our isochrone reconstruction to conduct a series of ice-sheet-wide and regional analyses to quantify the natural variation in retreat rates and areal extent change and to address scientific questions on the nature and style of retreat dynamics in Greenland during the Late Glacial to mid-Holocene period. We envisage that modellers can use this reconstruction to quantitatively compare deglacial simulations of the ice sheet with a vast quantity of empirical evidence. For example, the reconstruction (made available in both shapefile and NetCDF formats) could be used to quantitatively score model simulations run as part of ensemble experiments and/or to narrow down parameter spaces when conducting model sensitivity analyses, by testing the fit with margin extent and retreat pace (e.g. Patton et al., 2017; Ely et al., 2019; Pittard et al., 2022).

Figure 1 Greenland geographical context, present-day ice sheet thickness (BedMachine v4; Morlighem et al., 2017), and location of 18 regions (boxes outlined in light green), where detailed ice marginal landform mapping (displayed in purple) was conducted as part of this study. Note that few areas were mapped in detail in the southeast due to relative lack of ice-free land. A simplified pattern of contemporary ocean circulation is represented, after Yang et al. (2016; Fig. 1), where red and orange arrows symbolise warmer Atlantic-origin water, while blue arrows represent colder Arctic-origin water. NAC, IC, EGC, WGC, BC, and LC stand for the North Atlantic, Irminger, east Greenland, west Greenland, Baffin, and Labrador currents, respectively. The digital elevation model (DEM) displaying both topographic and bathymetric data is from the General Bathymetric Chart of the Oceans (GEBCO) 2022 release (450 m resolution). Geographical extent of Greenland periphery ice caps and glaciers (white polygons) is from BedMachine v4. The locations of main Greenland ice cores (further mentioned in text) are shown by black stars. All data (here and in subsequent figures) are displayed projected to the WGS 84/NSIDC Sea Ice Polar Stereographic north coordinate reference system.

2 Methods

2.1 Geomorphological reconstruction of ice retreat patterns

Across Greenland's exposed land area there has been a rich but sporadic field collection of landform and sedimentary evidence for former ice margin positions. These studies have been based on field investigation and mapping from aerial photographs or satellite images (e.g. Weidick, 1968, 1971; Ten Brink and Weidick, 1974; van Tatenhove et al., 1996; Roberts et al., 2009; Levy et al., 2012; Young et al., 2013b) and, more recently, using digital elevation models (DEMs; e.g. Carrivick et al., 2017; Pearce et al., 2018). Such investigations have typically mapped individual moraines, meltwater channels, and trimlines to define former ice margin positions, which were then used to build local reconstructions of ice margin retreat. This information covers a small fraction (< 10 %) of the land area, leaving most of Greenland's ice marginal history unexplored. The release of the 2 m ArcticDEM (Porter et al., 2018) for the whole of Greenland provides a consistent dataset that has the potential to revolutionise our understanding of Holocene ice retreat from the Greenland coast to the present-day ice sheet margin position. Motivated by the availability of this new resource, we devised a mapping scheme and protocol that could capture the first-order pattern of retreat for the whole of Greenland's terrestrial ice sheet periphery, an area of 430 500 km². We did not identify and map every ice marginal landform from the ArcticDEM, as this would represent a decade-long task. Instead, we sampled the area in sufficient detail to provide the first landform-based reconstruction.

2.1.1 Geomorphological mapping

We focussed our investigation on landforms indicating the position and shape of former Greenland Ice Sheet margins, paying less attention to peripheral ice caps and mountain glaciers. To capture high-resolution details of ice marginal retreat across the whole area and with a consistent approach is a challenge. To accomplish this, we adopted a sampling approach in which, for some regions that we called “swathes”, we investigated and mapped ice marginal landforms in detail (e.g. Figs. 1–3). For the intervening areas between swathes, we only identified the most prominent landforms that would permit us to connect paleo ice margins across these areas, joining up swathe to swathe (see Sect. 2.1.2). A total of 18 swathes were chosen, based on presenting especially dense and well-preserved ice marginal landforms (Fig. 1), to be positioned approximately evenly around the ice sheet perimeter and which covered ∼ 60 % of Greenland's ice-free periphery. Although it was beyond the scope of this continental-scale investigation to map all ice marginal landforms identifiable from the ArcticDEM, such an endeavour would be a valuable future goal. We suggest that the main focus should be to pursue mapping efforts in between the 18 swathes, which could contribute to updating future versions of the PaleoGrIS database.

All geomorphological mapping was conducted using remotely sensed data. Identification of ice marginal landforms was carried out using the 2 m ArcticDEM (Porter et al., 2018) and the 30 m ALOS WORLD 3D (AW3D30) DEM (AW3D30, https://www.Eorc.jaxa.jp/ALOS/en/aw3d30/, last access: 1 December 2022). As recommended by Smith and Clark (2005), optimised orthogonal hillshaded relief models of 315 and 45° azimuth angles and 45° inclination were toggled between to minimise azimuth biasing and better identify landforms of interest. Moreover, to aid us in distinguishing glaciogenic sediments from bedrock features, Google Earth Pro software was used for consulting three-dimensional visualisations of satellite imagery (as recommended by Chandler et al., 2018). All mapped landforms were digitised manually as shapefiles in the WGS 84/NSIDC Sea Ice Polar Stereographic north (EPSG code 3413) reference coordinate system. To remain consistent throughout, and to avoid introducing bias from other investigators' landform interpretations, our mapping was conducted without input from previous mapping investigations. Furthermore, mapping was not conducted at a fixed spatial scale but by zooming in or out to enable better visualisation of the spatially variable landform details.

The main indicators of former ice extent were terminal and lateral moraine ridges, more expansive moraine complexes, hummocks and hummocky ridges, lateral meltwater channels, and trimlines. These landforms were interpreted and identified based on their morphology and texture, their position in relation to the wider topographic setting, and following criteria for landform interpretation as detailed in Benn and Evans (2014), Chandler et al. (2018), Barr and Clark (2009), Rootes and Clark (2020), and Leger et al. (2020). The most numerous landforms were moraine ridges, which are typically discerned on DEM hillshades as arcuate, steep-sided, and sometimes sinuous ridges with positive relief and often displaying sharp crests (e.g. Fig. 2). Given our aim of covering the entire landmass in a single pass with a small number of investigators, we did not individually map and classify landforms into their respective types (e.g. moraine ridge and meltwater channel). Instead, we captured and summarised information from these landforms by digitising lines all grouped into a single layer (shapefile) called “ice marginal landforms” (Figs. 2, 3). Due to a lack of pre-Last-Glacial-Maximum (LGM) ice margin chronologies across Greenland, and following the assumption that glaciogenic deposits relating to previous glaciations were overridden during the last glaciation (Funder et al., 2011) and thus less likely to be preserved (with exceptions, see Mejdahl and Funder, 1994; Kelly et al., 1999), we adopt the assumption that our mapped ice marginal landforms were deposited during the last deglaciation between ∼ 17 kyr BP and the present.

Figure 2 Comparison between DEM hillshade (light azimuth 315°; incline 45°) of the 2 m spatial resolution ArcticDEM (Porter et al., 2018) in the left-hand panel, and the same area with our ice marginal landform mapping superimposed (right-hand panel). The area presented here is part of J. C. Christensen Land in north Greenland, as highlighted by the red box in the inset. This region displays remarkable preservation and spatial density of terminal and lateral moraine rides, moraine complexes, lateral meltwater channels, and trimlines.

2.1.2 Establishing an internally consistent map of retreat pattern

The high number (n=194 302) of ice marginal landforms identified inhibited reconstruction at an ice-sheet-wide scale. Thus, the landform record underwent two stages of simplification to create a pattern of retreat, following the method employed by Clark et al. (2012). First, ice marginal fragments were drawn to summarise the collective pattern in which ice marginal landforms could be reasonably linked in close proximity (thick black lines; Fig. 4). As a rule, ice marginal fragments were mostly confined within the same valley, unless the landform evidence was overwhelming for expanding beyond this scale. A further interpretive, and more speculative, step involved joining up these ice marginal fragments with ice marginal connectors (thick yellow lines; Fig. 4). These lines attempt to connect nearby ice marginal fragments interpreted as being approximately time-synchronous. Such connectors are guided by less prominent or less dense spreads of landform evidence and by considering relationships between topography and plausible ice dynamics. As a result, the 18 swathes featuring more detailed mapping (Fig. 1) enable us to connect landforms with both ice marginal fragments and connectors, while inter-swathe areas are mostly dominated by ice marginal connectors. To interpolate ice marginal connectors across offshore areas, submarine topographic data were obtained from the 15 arcsec spatial resolution General Bathymetric Chart of the Oceans (GEBCO) 2022 release (https://www.gebco.net/data_and_products/gridded_bathymetry_data/, last access: 1 October 2022). Together, the ice marginal fragments and connectors depict the direction and relative age of ice marginal recession in undated time steps; this is the first-order pattern of ice marginal retreat (e.g. Fig. 4). The underlying assumption behind our first-order retreat pattern map is that retreat was generally considered to be monotonic, unless there was geomorphological evidence for time-transgressive margins (e.g. cross-cutting moraines).

Figure 3 Examples of mapped ice marginal landforms in two distinct regions displaying high density and preservation of ice marginal glacial deposits. The mapped landforms are displayed overlaying topographic data from the AW3D30 DEM (30 m resolution). The left-hand panel presents mapping in a region of north Greenland (also referred to as Wulff Land), while the right-hand panel focuses on the deglaciated region directly south of Disko Bay in central west Greenland.

Figure 4 Visual and cartographic description of methodology followed to produce an ice-sheet-wide retreat pattern from DEM and imagery data. Our procedure involves three main steps of data generalisation to incrementally summarise the raw topographic and imagery data to near-continuous ice marginal connectors interpreted as representing time-synchronous former ice sheet margins. The latter product is then used to inform our isochrone time slice reconstruction, which involves adding geochronological information from our compilation of TCN and radiocarbon-derived event ages.

2.2 Compiling geochronological evidence

We attempt to compile published ages relating to the Greenland Ice Sheet grounded margin retreat from ∼ 14 kyr BP through to the present day for the entirety of the Greenland domain (Fig. 1). We include ages from both terrestrial cosmogenic nuclide (TCN) surface exposure dating and radiocarbon dating methods. Other dating methods (e.g. luminescence and lichenometry) were not included in this compilation, as they have been applied much more sporadically for the specific objective of dating the former position of the ice sheet margin in Greenland.

2.2.1 Terrestrial cosmogenic nuclide surface exposure ages

Compilation of ages

The geochronological component of the PaleoGrIS 1.0 database compiles published TCN surface exposure ages, and associated metadata, produced with the aim of dating the deglacial and Holocene ice extent fluctuations in the Greenland Ice Sheet. Such ages estimate the timing of moraine emplacement and stabilisation, deposition of glacial erratics, or the abandonment of ice from exposed bedrock surfaces, which are all relevant for constraining the timing of ice retreat. These exposure ages were extracted predominantly from the maintained ICE-D Greenland online database (ICE-D Greenland, https://version2.ice-d.org/greenland/publications, last access: 1 December 2022). In addition, a complementary review of the existing literature was conducted to compile relevant studies not currently included in ICE-D Greenland. Exposure ages dating the former margin evolution of periphery ice caps and mountain glaciers were intentionally excluded from our collection. Any new data added to ICE-D and any new study published later than our census date of 21 October 2022 are not included in the PaleoGrIS 1.0 database.

Age calibration and re-calculation

All ages were re-calculated using the online calculator formerly known as the CRONUS-Earth online calculator version 3 (Balco et al., 2008). We do not apply an erosion rate correction to our exposure age calculations. We use the LSDn scaling scheme (Lifton et al., 2014), and the ¹⁰Be west Greenland production rate (Young et al., 2013a) obtained from the ICE-D online calibration database (http://calibration.ice-d.org/, last access: 1 December 2022). Furthermore, no corrections were applied for post-exposure isostatic and/or tectonic uplift or subsidence, which, given the young nature of the ages (< 14 kyr BP) and the relatively low magnitude of surface elevation change (50–120 m) related to glacial isostatic adjustment during the Holocene, is not thought to cause age offsets greater than analytical uncertainties in most Greenland regions (Jones et al., 2019). To remain consistent across the domain, no corrections were applied for post-exposure vegetation or snow cover. Final exposure ages are reported in calendar years (or thousands of years, kyr) before present (BP), with “present” defined as the year of sampling. Nuclide concentrations and sample metadata were retrieved from ICE-D and/or from original publications when needed. When multiple ages were described by original publications as dating a single time-synchronous margin, these were grouped and summarised by a “summary event age” (Fig. 5) and uncertainty calculated using the population arithmetic mean and the 1σ standard deviation, respectively. No grouping of ages was applied if not clearly described in source study as dating the same event. While TCN exposure ages positioned on mountains can act as former ice thickness indicators (e.g. elevation dipstick models), they are not used directly by our isochrone reconstruction but are included in the database, as they may prove useful for comparing with modelled ice sheet thickness predictions, for instance.

Age filtering and quality control

Stratigraphical and/or statistical outliers were removed from event–age calculations only when considered such in original publications. The only exception to this is when sample coordinates provided in original publications were either missing or erroneous (e.g. plotting offshore or within the ice sheet). In this case, the age was considered an outlier and was not included in summary event age calculations. For a few more recent investigations (e.g. Søndergaard et al., 2020), in situ cosmogenic ¹⁴C was also measured alongside ¹⁰Be or ²⁶Al. When ¹⁴C-derived exposure ages display younger exposure ages and are described by authors as presenting less nuclide inheritance than other nuclides, they were given priority in our summary event age calculations. Our database also includes the details of whether paired nuclide analysis on a given sample (e.g. ²⁶Al $/$ ¹⁰Be) suggests a complex exposure/burial history and thus guides us to exposure ages that are likely too old to be considered in our Holocene retreat reconstruction. Furthermore, to help combine the geomorphological and geochronological evidence in a time slice reconstruction of ice sheet evolution, TCN exposure ages compiled within PaleoGrIS 1.0 are provided a quality control rating and classified in three categories described as high-, medium-, and low-confidence ages. The criteria list followed to apply this quality control is presented in Table 1 and is derived from the investigations of Hughes et al. (2016), Small et al. (2017), and Davies et al. (2020).

Table 1 Quality control assessment criteria list for TCN and radiocarbon summary event ages compiled in the PaleoGrIS 1.0 database. Criteria are adapted from Hughes et al. (2016), Small et al. (2017), and Davies et al. (2020) to fit the Greenland-specific context of numerical dating of former ice sheet margin retreat. Note that AMS is for accelerator mass spectrometer; n/a is for not applicable.

2.2.2 Radiocarbon ages

Compilation of ages

The PaleoGrIS 1.0 database also features a collection of Greenland-wide radiocarbon ages and associated metadata, which provide minimum-limiting age estimates of organic deposition in ice-free conditions following ice retreat during the deglacial and Holocene periods. This collection of radiocarbon ages was assembled by consulting former ice-sheet-scale reviews by Bennike and Björck (2002), Dyke (2004), Sinclair et al. (2016), and Dalton et al. (2020). Other more regional reviews, including those by Rinterknecht et al. (2014), Dyke et al. (2014), Larsen et al. (2014), and Young et al. (2021), were also examined. Furthermore, a review of the existing literature was conducted with the aim of finding other relevant studies not included in the above. Any new radiocarbon date published after 21 October 2022, our census date, is not included in the PaleoGrIS 1.0 database.

Age calibration and re-calculation

All radiocarbon ages were consistently recalibrated using the IntCal20 curve (Reimer et al., 2020) for terrestrial samples, the Marine20 (Heaton et al., 2020) curve for marine samples, and the CALIB 8.2 online calibration software (CALIB, http://calib.org/calib/calib.html, last access: 1 December 2022). We report final calibrated ages and uncertainty as the midpoint ± half of the calibrated age range at 95 % probability (2σ). Calibrated ages are reported in calendar years (or kyr) before present (BP), with “present” defined as year 1950.

For marine samples, we apply a marine reservoir age correction protocol that attempts to consistently account for both spatial heterogeneity in the reservoir effect itself and for variability between ΔR calibration sites. To do so, we calibrate all marine ages against the Marine20 curve using regional ΔR values obtained from the maintained online Marine Reservoir Correction Database (http://calib.org/marine/, last access: 1 December 2022) (Reimer and Reimer, 2001). For each sample and location, a final ΔR value is obtained by computing the weighted mean of the 10 nearest available ΔR calibration sites, as determined directly from the correction database. The reported final uncertainty following this calculation is the maximum of the standard deviation of ΔR and the weighted uncertainty in mean of ΔR (Bevington, 1969). Following this protocol, ΔR values in our database range from −113 to 73 years, while ΔR uncertainties range from 36 to 150 years. We note that these ΔR values mostly overlap with the newest Greenland-specific marine reservoir age correction assessment of Pearce et al. (2023), which was submitted after our compiled ages were re-calculated). We acknowledge that, for polar latitudes (> 50° N), calibrating marine radiocarbon ages against the Marine20 curve may be problematic, due to greater variability in ocean ventilation and air–sea gas exchange caused by fluctuations in sea ice extent and wind strength, leading to increased and more time-variable marine reservoir effects (Butzin et al., 2005; Heaton et al., 2022). However, this is more likely to be problematic during glacial periods (Reimer et al., 2020). For polar samples dating to the Holocene (11.5–0 kyr BP), Heaton et al. (2022) recommend calibrating directly against Marine20. Since the PaleoGrIS 1.0 reconstruction spans the Late Glacial and early- to mid-Holocene (14–6.5 kyr BP) period and because 90 % of the calibrated radiocarbon ages compiled in our database are younger than 11 kyr BP, we choose to treat all samples the same, for consistency.

Age filtering and quality control

Previous reviews and published studies were systematically filtered so that only radiocarbon dates produced with the aim to constrain the evolution of the ice sheet margin through time were compiled. Following this logic, only dated events that present a clear stratigraphic link to the evolution of the former ice sheet extent are incorporated in the database. In Greenland, a common example of such dates are multiple radiocarbon dates down a core of lacustrine sediments featuring a clear sedimentological transition from subglacial (e.g. till) to proglacial (e.g. silts and clays) deposits. In the latter case, and for any deglacial chronologies presenting multiple radiocarbon ages from a single core or single location, the oldest age of the sequence was considered the closest age estimate of ice-free conditions and organic deposition following ice retreat and was thus retained as the summary event age (following Hughes et al., 2016; Small et al., 2017; Dalton et al., 2020) (Fig. 5). The data were also filtered so that radiocarbon ages dating the evolution of Greenland periphery ice caps and mountain glaciers were not included in this collection. Stratigraphical and/or statistical outliers were not included in the summary event age compilation when described as such in original publications. Samples featuring missing or erroneous geographical coordinates were also considered outliers. Radiocarbon ages compiled within PaleoGrIS 1.0 are given a quality control rating and are classified in three categories described as high-, medium-, and low-confidence ages. The criteria list followed to apply this quality control consistently through the dataset is displayed in Table 1.

Figure 5 Spatial distribution of radiocarbon-derived (a) and TCN-exposure-age-derived (b) summary event ages (after filtering and statistics) used to produce the PaleoGrIS 1.0 isochrone reconstruction. Topographic and bathymetric data (GEBCO 2022 release) are overlaid with 1000 m interval contour lines. The dashed line shows our division scheme of the ice sheet in seven main drainage regions that we refer to throughout the paper and in associated quantitative analyses, after the drainage basin sampling scheme of Rignot and Mouginot (2012), labelled in panel (b). Panel (b) also features lines denoting our estimations of maximum and minimum scenarios of the last full glacial Greenland Ice Sheet extent (grounded ice only), based on a review of the literature (see Sect. 2.6), and which informs the mapping of our outermost isochrones.

Table 2 List of publications and associated numbers of compiled ages in PaleoGrIS 1.0 for TCN surface exposure dating.

Table 3 List of publications and associated number of compiled ages in PaleoGrIS 1.0 for radiocarbon dating.

2.2.3 The PaleoGrIS 1.0 geochronological database format

TCN and radiocarbon ages compiled were entered into two respective Excel (.xlsx) spreadsheets made available in the PaleoGrIS 1.0 database (data available from https://doi.org/10.17632/nh57cz4gys.1, Leger et al., 2024). Both spreadsheets (one for TCN and one for radiocarbon) document sample information and source publication details (also in Tables 2 and 3), metadata relevant to age calibration and re-calculations, all event identification and summary event ages, and age quality control attributions. A subset of these details including sample names, locations, source publications, summary event ages and uncertainties were used to generate point shapefiles for use with any geographic information system software (e.g. ArcMap or QGIS). Further details concerning the several shapefiles published alongside this work are described in the ReadMe files provided in the online database.

2.3 Producing Greenland-Ice-Sheet-wide isochrones

Within the context of reconstructing past ice sheet extent, isochrones are defined as time-stamped and spatially continuous margins highlighting changes in the former spatial extent of an entire ice sheet over time. Reconstructing such continuous perimeters is a challenge, given a fragmentary evidence base. Geochronological and geomorphological data available to empirically reconstruct isochrones are either point data (ages) or fragments of line data (e.g. moraine mapping). Such data are characterised by highly variable temporal and spatial densities and are associated with uncertainties of their own. Therefore, the task of drawing time-stamped and spatially continuous ice sheet perimeters involves interpolating between empirical evidence and extrapolating across blank areas of the map (e.g. Stroeven et al., 2016; Hughes et al., 2016). The following paragraphs describe the methods followed to produce isochrones that combine geomorphological and geochronological evidence, while separately accounting for temporal and spatial uncertainties in empirical data.

2.3.1 Isochrone time span, temporal resolution, and uncertainty

For PaleoGrIS 1.0, we map isochrones that delineate grounded ice only (and not floating ice fronts) (Fig. 6). Since Greenland ice-free land areas were deglaciated mostly during the Late Glacial and early- to mid-Holocene periods, there is a higher density of terrestrial TCN and radiocarbon dates displaying ages between 12 and 6.5 kyr BP (Fig. 7). Following the Holocene thermal maximum, between 6.5 and 2 kyr BP, the number of available dates drop significantly (Fig. 7). This coincides with the ice sheet reaching its minimum Holocene extent (between ∼ 7 and ∼ 4.5 kyr BP), which, in most Greenland regions, was either as extensive as or more retreated than today's margin position (Larsen et al., 2015; Briner et al., 2016). After ∼ 5–4 kyr BP, the onset of neoglacial cooling caused ice sheet readvances culminating in the Little Ice Age, which explains the small relative increase in the number of dates in our compilation after 2 kyr BP (Fig. 7). Consequently, the time period featuring enough terrestrial geochronological evidence for reconstructing past retreat at the ice sheet scale currently spans only 6–7 kyr. Recent methodological and technological improvements to TCN exposure and radiocarbon dating now enable the production of Holocene ages with analytical uncertainties that can be inferior to 500 and 200 years (at 1σ), respectively. For these reasons, we chose to produce ice-sheet-wide isochrones at 500-year temporal resolution between 12 and 6.5 kyr BP.

Between 14 and 12 kyr BP, a period characterised by a predominantly marine-terminating ice sheet constrained by fewer empirical data, we chose to draw isochrones at a 1000-year temporal resolution. Prior to this, between the ice sheet's last full glacial extent and 14 kyr BP, the ice sheet is thought to have been mostly grounded on the presently submerged continental shelf (Funder et al., 2011). A few offshore sampling studies provide geochronological constraints for the approximate location of the grounded ice margin during this time period (e.g. Smith and Licht, 2000; Kuijpers et al., 2003; Ó Cofaigh et al., 2013). However, we believe such studies are currently too scarce and spatially scattered to enable tracing ice-sheet-wide isochrones between the last full glacial extent and 14 kyr BP and hope such improvements can be made with future versions of this reconstruction as more data from less studied regions arise. The PaleoGrIS 1.0 reconstruction therefore features 14 isochrones (or time slices) between 14 and 6.5 kyr BP.

In an attempt to account for uncertainties inherent to TCN exposure and radiocarbon dating, we chose to allocate a time range to individual isochrones and thus to each time slice of our reconstruction (Fig. 6). For instance, the youngest and innermost isochrone of PaleoGrIS 1.0 is here referred to as the 7–6.5 kyr BP isochrone. This means that we estimate the former margin to have been near the location of that isochrone line at any time between ∼ 7 and ∼ 6.5 kyr BP. This line was thus drawn to connect, as much as possible, the landforms and summary event ages comprised between 7.0 and 6.5 kyr BP when rounded to the nearest 100 years. This approach is different from previous time slice reconstructions that attributed a single timestamp to isochrones and choosing to represent the dating uncertainty spatially by differentiating minimum, maximum, and/or optimum positions for isochrones for a given time slice (e.g. Dyke and Prest, 1987; Hughes et al., 2016; Davies et al., 2020; Dalton et al., 2020; Clark et al., 2022). Contrastingly, the PaleoGrIS 1.0 approach aims to more clearly distinguish and separate the temporal from the spatial uncertainties inherent to isochrone reconstructions. Isochrone temporal uncertainty is exclusively associated with the analytical and calculation/calibration uncertainties of the numerical dates (TCN or radiocarbon), while isochrone spatial uncertainty instead results from spatially variable density of geochronological and geomorphological evidence. Here, the latter is treated by attributing various confidence levels along a single isochrone line (more details in Sect. 2.3.5). Therefore, our approach allows comparing numerical ice sheet model outputs to a single isochrone, while enabling model time to vary within our isochrone temporal error range and thus account for analytical uncertainty in geochronological data. To help modellers use the PaleoGrIS 1.0 isochrones in their model–data comparison procedures, our online database includes details (ReadMe files) on how various data formats could be used, depending on the model resolution and the type of experiment conducted.

2.3.2 Rules followed when drawing isochrones

To draw spatially continuous isochrones as consistently as possible around the ice sheet's periphery, the following workflow and set of rules were applied.

• Using ArcMap 10.7.1 software, we displayed our landform mapping database accompanied by our synthesised pattern of retreat map.

• All TCN and radiocarbon summary event ages were displayed rounded to 100 years and feature a traffic light colour code relating to one of the three confidence level categories (Table 1).

• All information was displayed on a rendition of topography from the ALOS World 3D 30 m spatial resolution DEM and bathymetry from the GEBCO 2022 release.

• The location of contemporary grounded ice was displayed at all times using the raster mask of the IceBridge BedMachine Greenland version 4 dataset (Morlighem et al., 2017).

• Isochrones were interpreted and mapped, working clockwise around the full Greenland perimeter, and sequentially following chronological order (starting from oldest). The process was conducted iteratively and with numerous adjustments, as the position of an individual isochrone might depend upon the preceding and succeeding ones. Multiple authors re-interpreted and contributed to the final isochrones to try and ensure a consensus view of possible alternate behaviours.

• When the ice sheet was more extensive than present everywhere, isochrones were drawn around the whole of Greenland. However, during younger isochrone time slices (i.e. 10–6.5 kyr BP), paleo ice sheet margins were likely similar to present-day positions, or in a more retreated position, for certain regions (Briner et al., 2010; Larsen et al., 2015). In such cases, isochrones were not drawn around the full ice sheet perimeter but were instead interrupted where they meet the present-day margin. We do not attempt to map out the extent of ice where it retreated inside of the present-day position, as the position of the ice margin is largely undefined/unknown, even within regions where it is known to have occurred.

The precise method employed when drawing isochrone lines varied, depending on the nature and spatiotemporal density of empirical data available. When both mapped ice marginal landforms and TCN or radiocarbon event ages coincided, the line was drawn through both sets of information. When only reliable event ages were available, the line was drawn to connect them while considering the topography, bathymetry, and spatial configuration of the modern ice sheet margin (e.g. Fig. 6). On the other hand, when only mapped ice marginal landforms were available, the line was drawn following the landforms and topography/bathymetry only. In the absence of any chronological constraints, the retreat pattern was assumed to be monotonic in nature, i.e. producing a decreasing extent that is spatially and temporally consistent. The resulting ice margin isochrones are thus much more complex than a product of lines joining up geochronological point data (Fig. 6). Instead, the PaleoGrIS 1.0 isochrones represent a qualitative reconciliation and interpretation of topographic, geomorphological, and geochronological data. Such a heuristic approach sometimes relies on soft knowledge of typical interactions between ice sheet margins and topography, including the lobate behaviour of outlet glaciers and the dynamics of ice flow around present-day Greenland. The isochrones therefore stand as an informed interpretation of successive ice margin positions, well constrained in some places and in others much less so, and which are likely to require adjustment once further landform or geochronological data become available.

Figure 6 Example maps showing the details of the PaleoGrIS 1.0 isochrones in a land-dominated region of central west Greenland (a) and a fjord-dominated region of central east Greenland (Scoresby Sund (b)). The two examples highlight our choice of distinguishing isochrone temporal uncertainty (shown by time ranges labelled in bold black) from isochrone spatial uncertainty (shown by the four colour-coded confidence levels). The figure also displays the location of local event ages compiled in our database which is colour-coded here based on our quality control assessment (Table 1). All data are displayed with the AW3D30 DEM overlaid (30 m resolution).

2.3.3 Addressing periphery ice caps and glaciers

During the Late Glacial and early- to mid-Holocene, as the Greenland Ice Sheet margin receded, peripheral ice caps and mountain glaciers became separated from the main ice sheet body. Ideally, a full reconstruction of the ice sheet Late Glacial and Holocene evolution would also include the changing extent of these peripheral ice bodies. Future versions of PaleoGrIS could usefully be extended to include the development of such peripheral ice caps, but for this first version, we chose to exclusively reconstruct the retreat pattern of the main body of the ice sheet. Peripheral ice caps were thus either (i) included within the perimeter of the reconstructed Greenland Ice Sheet margin, when the latter was considered extensive enough, or (ii) excluded from our reconstruction, after complete separation from the ice sheet was estimated to have occurred.

2.3.4 Connection with the Innuitian Ice Sheet

Ellesmere Island lies just 100 km to the NW of Greenland (Fig. 1), and along with much of the rest of the Canadian Arctic Archipelago, it was covered by the Innuitian Ice Sheet, which is thought to have been connected to the Laurentide Ice Sheet until approximately 10.5–9.5 kyr BP (England et al., 2006). Given such a close proximity, it is deemed likely that the Greenland and Innuitian ice sheets coalesced over Nares Strait during glacial maxima (Dyke, 2004; Sinclair et al., 2016; Funder et al., 2011; Georgiadis et al., 2018; Dalton et al., 2020). We chose to make no inferences regarding the pattern of retreat associated with the margins of the Innuitian Ice Sheet. We thus draw isochrones that are consistent with the two ice sheets merging, but we interrupt them towards the point of deepest bathymetry between Greenland and Ellesmere Island throughout the length of Nares Strait. The same method was applied to close polygon shapefiles over Nares Strait when measuring the ice sheet's areal extent and create time slice maps.

2.3.5 Isochrone confidence levels

Geochronological and geomorphological evidence are spatiotemporally heterogeneous, causing each isochrone to feature spatially variable levels of reliability around its perimeter (e.g. Fig. 6). Given this uneven distribution of information, our mapping procedure was adapted accordingly. We split each isochrone line into four confidence level categories referred to as high, medium, low, or very low confidence isochrones, by applying the following set of rules.

• High-confidence isochrone sections. These are drawn when empirical evidence features both mapped ice marginal landforms and TCN or radiocarbon event ages and deemed reliable, i.e. graded high or medium confidence.

• Medium-confidence isochrone sections. These are drawn when empirical evidence does not feature mapped ice marginal landforms but features reliable event ages, i.e. graded high or medium confidence.

• Low-confidence isochrone sections. These are drawn when empirical evidence features only ice marginal landforms or only low reliability dates, i.e. ages graded low confidence.

• Very low-confidence isochrone sections. These are drawn in the absence of any mapped geomorphological or any geochronological evidence but relying on topographic and/or bathymetric information exclusively.

2.3.6 Ages that could not reasonably be reconciled to isochrones

When building isochrones by integrating geochronological evidence from numerous published sources and locations, contradictory evidence can become challenging to resolve. For instance, this can occur when nearby event ages display high age variability while presenting similar levels of reliability but also when they are found reversed relative to the general direction of presumed ice retreat, and/or relative to stratigraphic order of events in nearby valleys/regions.

For TCN exposure dating, misleading ages can result from nuclide inheritance causing apparent ages that are too old, post-depositional disturbance causing apparent ages that are too young, or laboratory contamination/errors potentially causing both (Dunai, 2010). Such issues are in most cases reported in original publications, thus enabling us to easily identify outliers. In Greenland more specifically, Late Glacial and Holocene studies generally report nuclide inheritance to be the dominating cause of exposure age scatter (Larsen et al., 2021). For radiocarbon ages, delayed organic growth following ice retreat or contamination by younger organics (e.g. for bulk samples) can cause ages to be too young. This challenge is common in the Arctic, where post-deglaciation biodiversity establishment is relatively slow. For marine ages, an erroneous estimation of the marine reservoir effect can also cause calibrated ages to misrepresent the true deposition age. Moreover, in Greenland, numerous radiocarbon ages were produced to date the marine-to-freshwater transition of basins following isostatic uplift caused by ice retreat (Weidick et al., 2004). Delay in such uplift can increase the potential for radiocarbon ages underestimating the timing of deglaciation. The misfit between our isochrone reconstruction and a specific event age may also be due to oversights in our age-filtering exercise.

In these challenging cases, drawing isochrones requires subjective decisions to either favour and/or ignore certain seemingly awkward or anomalous event ages when weighed against our interpretation of the most representative timing of local deglaciation. Such interpretations were either based on a review of the region-specific literature or on assessing the distance separating relevant isochrones in adjacent regions displaying more consistent geochronological constraints. Summary event ages deemed challenging and ignored were identified and gathered after the mapping of isochrones in a separate shapefile available in the online database (https://doi.org/10.17632/nh57cz4gys.1, Leger et al., 2024). We acknowledge that decisions to ignore event ages are necessarily subjective (see Sect. 2.3.2).

2.3.7 Limitations and uncertainties in isochrone reconstruction

In regions of Greenland with a high quality and density of geochronological and geomorphological constraints (e.g. the Sisimiut and Disko Bay regions), former ice margins can be reconstructed with reasonably high levels of confidence. Contrastingly, our reconstruction also features crude interpolations of ice margins over vast areas, due to the low density or absence of empirical evidence in numerous locations. Thus, the heterogeneous nature of our reconstruction's uncertainties, as depicted by our four different confidence levels (Sect. 2.3.5), calls for caution. One must stress that the PaleoGrIS 1.0 isochrone reconstruction is not intended for use at the valley or regional scale, for which local investigations are likely more accurate (e.g. Pearce et al., 2018), but rather at continental scale. Caution should also be applied concerning the temporal resolution of our reconstruction, which aims to provide estimates of ice sheet margin positions every 1000–500 years. In fact, our reconstruction should be regarded as averaged or net retreat over such time intervals, and any short-lived readvances or dynamic ice sheet margin response (e.g. tens to hundreds of years) are not captured in our reconstruction. In most locations, the PaleoGrIS 1.0 isochrones should not be considered precise enough to accurately depict the margin evolution of individual outlet glaciers, to predict the formation of small proglacial lakes and spillways, or to reconstruct past meltwater pathways, for instance. The reliability of the reconstruction is, of course, dependent on the underpinning evidence, an examination of which should guide a user away from more speculative areas. The continental-scale mapping approach can sometimes cause resulting ice configurations to be glaciologically unrealistic and inconsistent with a complex landscape at the valley scale (< 5 km spatial resolution), despite our efforts to consider the influence of local topography/bathymetry as much as possible.

Over the past 30 years, several offshore investigations have established the approximate location of the ice margin during initial deglaciation from its local full glacial position (∼ 17–12 kyr BP) when the ice sheet was mostly marine-terminating (e.g. Smith and Licht, 2000; Nørgaard-Pedersen et al., 2008; Evans et al., 2009; Ó Cofaigh et al., 2004, 2013; Rasmussen et al., 2022; Hansen et al., 2022; Lloyd et al., 2023). Such studies remain scarce, and we find former grounding line positions remain largely understudied and undated. Readers should thus be aware that our oldest isochrones (i.e. between ∼ 14 and ∼ 12 kyr BP) are associated with crude interpolations across large offshore areas, thus presenting higher levels of uncertainty than we now have for the terrestrial areas.

We believe there is much scope for reducing isochrone uncertainty in future versions of the PaleoGrIS reconstruction, as more mapping and dating are conducted from less studied regions. We make suggestions as to which regions would most benefit from this attention in a discussion section later in the paper.

2.4 Estimating areal extent change in the Greenland Ice Sheet

To evaluate the former areal extent of the ice sheet (i.e. two-dimensional surface area) between ∼ 14 and ∼ 6.5 kyr BP, we produced a series of polygon shapefiles covering the area delineated by each isochrone perimeter and measured their areal extent using the ArcMap geometry calculator. To compare variations regionally, we divided the ice sheet into seven major hydrological basins, following the ice divide sampling of Rignot and Mouginot (2012) and the IMBIE Team (2019) (Fig. 5). Region-dividing polylines were extended towards the outermost isochrones, following topographic/bathymetric highs, while remaining as perpendicular to isochrones as possible. Therefore, for the exclusive purpose of subdividing the ice sheet into regions, the catchment areas were crudely assumed to remain similar from ∼ 14 kyr BP onwards, as their paleo configurations remain unknown. Between the 14–13 and 10.5–10 kyr BP time slices, we consider our measurements as absolute estimates of areal extent. However, for time slices between 10–9.5 and 7–6.5 kyr BP, our measurements should be regarded as maximum-limiting estimates of ice sheet areal extent. That is because empirical evidence suggests the ice sheet was as extensive as or smaller than today in several regions during that time interval. In these cases, we chose not to guess the extent of retreat behind the current margin and merged polygons representing more extensive margins with the contemporary ice sheet extent in data-free regions. Present-day ice sheet areal extent was computed from BedMachine Greenland version 4 (Morlighem et al., 2017) after removing peripheral glaciers.

We chose not to convert our ice sheet areal extent reconstruction into a volume and mass estimation. The Greenland Ice Sheet was characterised by very cold and arid conditions during extensive advance and, conversely, significant increases in accumulation occurred during deglacial margin retreat. The relationship between areal extent and ice thickness is thus complex and not necessarily positively correlated across the ice sheet (Cuffey and Clow, 1997). Obtaining a realistic deglacial volume reconstruction of the Greenland Ice Sheet requires extensive modelling experiments that take climate, surface mass balance, glacial isostatic adjustment, relative sea level change, calving, and basal sliding (e.g. Bradley et al., 2018) into account and that conduct quantitative model–data comparisons. Our team is conducting such an experiment with the Parallel Ice Sheet Model (Winkelmann et al., 2011) and will thus present volume reconstructions in a separate publication.

2.5 Assessing rates of retreat and their variation

To quantitatively estimate the former retreat rates of outlet glaciers, 72 transects emanating from present-day outlets were drawn. Transects were traced from the outermost (14–13 kyr BP) to the innermost (7–6.5 kyr BP) PaleoGrIS 1.0 isochrone margins, following our estimate of the former glacier front central position at each time step and guided by topography. We thus follow the assumption that the approximate centre of former glacier termini (or grounding lines) was located near the point of lowest topography. Hence, such transects should not be interpreted as supraglacial flowlines, as these would have likely evolved through time with glacier catchments and ice divides potentially readjusting during deglaciation. The transects are made available in a polyline shapefile (see the online database). Present-day outlet glacier names (preferably new Greenlandic, otherwise Danish) were obtained from the database of Bjørk et al. (2015), while distances along transects were measured for each isochrone spacing using ArcMap. As each isochrone is associated with a temporal uncertainty (i.e. 1000 or 500 years), both minimum and maximum retreat rates were calculated for each time step. The reported retreat rate for each step is taken as the midpoint ± half the range. For the purpose of evaluating former retreat rates only, the positions of isochrones are regarded as definite, and the reported error is exclusively temporal. For each transect, an overall retreat rate (± temporal uncertainty) is also computed using the total distance divided by the minimum and maximum time span between outermost and innermost isochrones.

2.6 Where and when was Greenland's last full glacial extent?

Although we mostly focus on the terrestrial deglaciation from ∼ 14 kyr BP onwards, it was important for positioning this outermost isochrone to have some knowledge of greater ice extents and indeed the maximum achieved extent in the last glacial. A review of Greenland's potential last full glacial extent (thought to occur around 18–16 kyr BP; Simpson et al., 2009) was conducted. We updated the full glacial margin drawn by Funder et al. (2011) after consulting publications contributing new empirical knowledge to that specific question (Möller et al., 2010; Ó Cofaigh et al., 2013; Arndt et al., 2017; Laberg et al., 2017; Jennings et al., 2017; Newton et al., 2017; Seidenkrantz et al., 2019; Sbarra et al., 2022; Couette et al., 2022; Rasmussen et al., 2022; Hansen et al., 2022). In several regions, a debate prevails regarding whether grounded ice reached the outer continental shelf when in the last full glacial configuration. We find that in all studied regions, apart from offshore the Northeast Greenland Ice Stream (NEGIS; Rasmussen et al., 2022), more recent investigations tend to suggest a more extensive full glacial configuration than previously proposed, with grounded ice often argued to have reached the mid- to outer shelf. For instance, in central west Greenland, Ó Cofaigh et al. (2013) used a series of dated marine sediment cores and bathymetric subglacial landform mapping to show that the Uummannaq and Jakobshavn Isbræ ice streams remained grounded near the continental shelf edge until ∼ 15 kyr BP. A similar interpretation was, for instance, proposed by Hansen et al. (2022), who argued that grounded ice in the Westwind Trough region (northeast Greenland) was located towards the outer shelf during last full glacial extent and prior to 13.5 kyr BP. However, near the same site (79.5° N), Rasmussen et al. (2022) suggest that grounded ice did not reach the shelf edge during the last glacial. Such contrasts demonstrate that the question remains open in this region. Vast areas of the Greenland continental shelf are still data scarce for this purpose. We revise the limit of Funder et al. (2011) by considering two scenarios, a minimum and maximum last full glacial extent, with the aim of highlighting spatial uncertainties in debated regions (Fig. 5). These two full glacial extent scenarios are used to inform the mapping of our outermost isochrones and are included here as shapefiles (see online database). The maximum extent scenario displays an ice sheet that reaches the shelf edge around the entire ice sheet perimeter. The minimum extent scenario modifies the latter by following the outline of Funder et al. (2011) in the central east and northeast regions, except towards the Westwind Trough, where it accounts for new data by Hansen et al. (2022). The minimum extent scenario also displays a more retreated full glacial ice sheet in the northwest region, where little data constrain whether a mid-shelf or outer shelf position was reached. We suppose that the grounding line was likely highly dynamic and that any maximum achieved last glacial extent configuration between these lines is feasible.

Figure 7 Histogram displaying the temporal distribution of all summary event ages calculated in the PaleoGrIS 1.0 database. Colour-coded relative proportions of the three quality control categories are shown for each 1000-year bin. A light-blue polygon highlights the time period covered by the PaleoGrIS 1.0 isochrone reconstruction. This Late Glacial to mid-Holocene period (14–6.5 kyr BP) features the large majority (> 90 %) of compiled event ages dating the deglacial evolution of former Greenland Ice Sheet margins.

3 Results

3.1 The PaleoGrIS 1.0 geomorphological and geochronological database

3.1.1 Geomorphology

The PaleoGrIS 1.0 database contains a total of 194 302 ice marginal landforms mapped over 430 500 km² of ice-free land around Greenland. The distribution of landforms is relatively sparse in the southeast (SE) and northwest (NW) regions, where the contemporary ice sheet terminates near to the coast. Consequently, the majority of ice marginal landforms mapped as part of this study are located in the southwest (SW), central west (CW), central east (CE), northeast (NE), and north (NO) regions of ice-free Greenland. Interestingly, we find ice marginal landforms are particularly abundant and well-preserved in NO Greenland, increasing our confidence here in the retreat pattern relative to other regions. This abundance and pristine nature of landforms might be due to the generally softer nature of subglacial bedrock (Pedersen et al., 2013) in this region enabling greater sediment supply for moraine-building, a drier climate that could better promote landform preservation (Niwano et al., 2021), or a potentially steadier and more monotonic retreat pattern of outlet glacier margins characterised by fewer dynamic fluctuations and readvances that cause erosion.

3.1.2 Geochronology

The PaleoGrIS 1.0 database features 1028 TCN exposure ages gathered from 45 studies and contains 423 radiocarbon ages collected from 111 studies. Following the computation of summary statistics for age groups, and the removal of outliers, this collection contains 251 TCN-derived summary event ages. This number excludes TCN exposure ages produced as past ice sheet thickness indicators (e.g. dipstick models). The radiocarbon age compilation features 370 summary event ages post age-filtering. Therefore, a total of 621 summary event ages were directly used to inform our isochrone reconstruction. Details regarding the temporal frequency distribution of summary event ages and their number per investigation are presented in Fig. 7 and Tables 2 and 3. We here acknowledge the possibility that relevant investigations and datasets may be missing from our chronological database. It is our goal that such an oversight can be corrected in future versions of PaleoGrIS.

In the PaleoGrIS 1.0 database, 90 % (n=558) of the summary event ages compiled are younger than 14 kyr BP and older than 6 kyr BP (Fig. 7). As a result, we chose to restrict our isochrone reconstruction to this time period (more justification for this decision can be found in Sect. 2.3.2). Following quality control assessment of summary event ages (Table 1), we find that high- and medium-confidence event ages each represent a quarter of the total (26 % and 25 %, respectively). Low-confidence event ages represent the largest proportion, i.e. 49 % (Fig. 7). This is in part explained by numerous event ages resulting from a single marine radiocarbon age, here attributed a low-confidence rating due to large uncertainties in past marine radiocarbon reservoir effects at high latitudes (Heaton et al., 2022) (see criteria list in Table 1). However, we stress that event ages described as low confidence may still closely estimate the true timing of margin retreat and are by no means excluded from our isochrone reconstruction. Our event age quality control assessment is merely an indicator utilised to inform challenging decisions made when tracing isochrones.

Using our geochronological database, we analyse the spatial variability in event age reliability across Greenland and its periphery. For TCN exposure ages, we find a higher relative concentration of medium- to low-confidence event ages in NO Greenland, and more specifically from Inglefield Land, Inglefield Bredning, Washington Land and around Danmark Fjord, Hagen Fjord, Independence Fjord, and the Centrumsø regions. Indeed, while studies that investigated these regions sampled both bedrock and erratic surfaces and have produced extensive datasets, their results often display old apparent ages (Ceperley et al., 2020; Larsen et al., 2018, 2020; Søndergaard et al., 2019, 2020). It has been stated that this abundance of overestimating exposure ages is most likely associated with high levels of nuclide inheritance in NO Greenland caused by insufficient subglacial erosion of bedrock and transported clasts. Future studies producing exposure ages from these regions might thus be inclined to also measure in situ cosmogenic ¹⁴C alongside other radionuclides (e.g. ¹⁰Be, ²⁶Al, and ³⁶Cl), which may enable the quantification of nuclide inheritance and more accurate estimation of Holocene deglacial exposure ages (e.g. Søndergaard et al., 2020). We note that the ice-free region to the west of the Akuliarutsip Sermia and Inuppaat Quuat glaciers (67.3–68.3° N, 50.1–54° W), in central west Greenland, also displays a high concentration of medium- to low-confidence TCN event ages (Young et al., 2020). For radiocarbon ages, however, we do not observe any clear spatial patterns in the variability in the summary event age reliability.

3.2 Regional retreat patterns and timings

In this section, we describe the PaleoGrIS 1.0 isochrone reconstruction and the pattern and timing of ice sheet margin retreat for selected regions, followed by a chronological description of the main deglacial events. This section also describes key empirical constraints that inform our reconstruction. When found, new Greenlandic place and glacier names are preferably used to describe the geography (e.g. after Bjørk et al., 2015). When not found, we refer to places, features and certain glaciers using Danish or other foreign names.

3.2.1 Ice sheet retreat in north Greenland and the timing of Nares Strait opening

Towards the onset of the studied time period (i.e. 14–13 kyr BP), our reconstruction portrays the Innuitian (ice over Ellesmere Island) and Greenland ice sheets as connected, with grounded ice from both ice masses merging along Nares Strait, as supported by empirical evidence (e.g. Jennings et al., 2011; Georgiadis et al., 2018). For that time step, we reconstruct the two connected ice sheets to feature marginal grounding lines on both the southwestern and northeastern ends of Nares strait, with respective margins terminating into Baffin Bay and the Lincoln Sea (Fig. 8). Geochronological evidence of ice sheet extent suggests that the retreat of these two grounding lines occurred somewhat simultaneously on both ends of the strait. This is supported by ages indicating onshore regions located closest to Nares Strait's openings (e.g. Inglefield Fjord and Wulff Land) were progressively deglaciated earlier than regions located towards its centre, i.e. terrestrial regions adjacent to Kane Basin, such as Washington and Inglefield lands (Fig. 8). Indeed, TCN exposure ages suggest grounded ice had retreated within Inglefield Fjord by around 11.5–11 kyr BP (Søndergaard et al., 2019). At the same time, deglacial radiocarbon ages (Funder, 1982; Kelly and Bennike, 1992) indicate ice from the Ryder, Steensby, and C. H. Ostenfeld basins had retreated towards the outermost present-day coastline of Nyeboe, Wulff, and Nares lands. After 11 kyr BP, we reconstruct grounding lines that retreated further within Nares Strait and that reached the southwestern and northeastern edges of Kane Basin by 9–8.5 kyr BP (Fig. 8). TCN ages produced from contemporary coastal regions of Inglefield Land (Søndergaard et al., 2020) and Washington Land (Ceperley et al., 2020) suggest that regions adjacent to the Humboldt Glacier lateral margins started deglaciating at around 8.5–8 kyr BP. During that time, further east, evidence suggests the Petermann Glacier front was located towards the fjord's mouth (Bennike, 2002). After ∼ 8 kyr BP, in this region, we estimate that the ice sheet had retreated inland beyond present-day coastlines, except for offshore contemporary Humboldt Glacier, where we map a grounding line that remained within ∼ 100 km outside the present-day margin. Therefore, empirical data suggest that the last ice bridge connecting the Innuitian and Greenland ice sheets over Nares Strait survived exclusively within the Kane Basin, and until 8.5–8 kyr BP (although large uncertainties remain). We thus estimate the final deglacial opening of Nares Strait to have occurred between 9 and 8 kyr BP (Fig. 8). Following this, further retreat caused the ice sheet margin to reach its present-day extent by 7.5–7 kyr BP on Inglefield Land, while shortly after (7–6.5 kyr BP) on Washington Land.

In northernmost Greenland (> 82.7° N), we reconstruct a margin retreat pattern characterised by a disconnection with glaciers and ice caps from Roosevelt Land that occurred at around 10–9.5 kyr BP. This estimation is however uncertain due to the low abundance of local geochronological constraints. Around 9–8.5 kyr BP, we estimate further retreat to generate ice sheet separation from the large Hans Tausen Ice Cap. Between 8.5 and 6.5 kyr BP, the relatively slow retreat of the more-extensive-than-present ice margins was still occurring on inter-fjord regions of Wulff, Adam Biering, and J. C. Christensen lands (Larsen et al., 2020). In the same regions, empirical evidence suggests outlet glaciers terminating in deep and wide fjords, such as Independence, Hagen, or Victoria fjords, experienced faster retreat than those located on adjacent inter-fjord regions (Larsen et al., 2020).

3.2.2 Ice sheet retreat in northeast Greenland and ice margin evolution in the Nioghalvfjerdsfjorden and Jøkelbugten regions

In this section, we describe the PaleoGrIS 1.0 isochrone reconstruction in a region characterised by a ∼ 300 km long stretch of the coast in northeast Greenland belonging to King Frederick VIII Land and which lies between Holm Land (79.8° N) and the Skaer Fjord (77.4° N). The largest contributor to ice flux in this region is the North East Greenland Ice Stream (NEGIS), which currently splits into three wide marine-terminating and fast-flowing outlet glaciers displaying surface velocities of > 1000 m yr⁻¹ in places (Joughin et al., 2018). These are, from north to south, the Nioghalvfjerdsbrae (also referred to as 79N) glacier, the Zachariæ Isstrøm glacier, and the southernmost branch which splits in two sub-outlets, the Kofoed-Hansen and Storstrømmen outlet glaciers (Fig. 9). In this region, a debate prevails regarding whether the grounded ice sheet advanced extensively on the wide continental shelf during last full glacial extent or whether it was restricted to the inner shelf (Rasmussen et al., 2022). However, recent offshore investigations (e.g. Hansen et al., 2022), including soon-to-be-published work (Ó Cofaigh et al., 2023), increasingly suggest a last full glacial ice sheet margin that was more extensive than previously drawn by Funder et al. (2011) and that likely reached the mid- to outer continental shelf. However, whether the grounded ice sheet remained confined to the prominent bathymetric Westwind (north) and Norske (south) troughs or instead also flowed eastwards over the mid-shelf bathymetric highland separating the two, the Northwind Shoal, is uncertain (Arndt et al., 2017; Pados-Dibattista et al., 2022). The large disparity between our minimum and maximum last full glacial extent scenarios reflects this debate in the literature (Fig. 9).

To remain conservative, our reconstruction in this region features a very low-confidence outermost isochrone (14–13 kyr BP) located closely on board our minimum last full glacial extent scenario, depicting a grounded ice sheet margin positioned towards the mid-shelf, while remaining confined within the Westwind and Norske troughs (Fig. 9). This 14–13 kyr BP extent is constrained in the Westwind Trough by a marine radiocarbon deglacial chronology by Hansen et al. (2022). We then tentatively reconstruct an ice sheet margin that had monotonically retreated ∼ 50 km westwards along the troughs by ∼ 12–11.5 kyr BP, as suggested by offshore data by Davies et al. (2022). However, directly south of Norske Trough, we reconstruct a less extensive ice sheet margin around the same time (12–11.5 kyr BP). This interpretation is based on TCN exposure ages by Larsen et al. (2018), indicating that by 11.5 kyr BP, the grounded ice had likely retreated west of Kap Amélie (77.5° N), into the Skærfjorden embayment and had reached the Storøen and Ambolten islands; this became a north–south-oriented archipelago acting as a topographic barrier to ice flowing eastwards from the deep Jøkelbugten basin (Fig. 9).

Based on TCN-derived event ages from Bourbon Øer (Larsen et al., 2018) and Dove Bay (Larsen et al., 2022) located ∼ 60 km further south, we reconstruct a margin that had retreated to the inner shores of Skærfjorden and that was positioned ∼ 55 km from the present-day ice front of Nioghalvfjerdsbrae by ∼ 11.5–11 kyr BP. Over the next 2000 years, our reconstruction suggests that the ice sheet margin retreated relatively quickly westwards within the Jøkelbugten basin, through Lamberts Land, and within Nioghalvfjerdsfjorden, with a former margin that was positioned within ∼ 15 km of the present-day ice sheet front in most locations of the region by ∼ 10–9.5 kyr BP. We estimate the next phase of retreat to be slower (< 20 m yr⁻¹), with an ice sheet margin remaining more extensive but near the present-day one, until ∼ 9–8.5 kyr BP. This is well supported by coeval TCN exposure ages from Bloch Nunatakker, an island located near Nioghalvfjerdsbrae's contemporary calving front, and from three sites situated ∼ 3 km from the modern lateral margins of Zachariæ Isstrøm glacier (Larsen et al., 2018). To the south of Nioghalvfjerdsbrae, and for a 300 km long stretch of the ice margin, we thus reconstruct an ice sheet margin that was as or more retreated than the present by ∼ 8.5 kyr BP. This contrasts with the mountainous region directly northwest of Nioghalvfjerdsbrae, however, where the predominantly land-terminating former ice sheet margin appears to have retreated more slowly and more steadily. This is supported by our mapping depicting a highly regular spacing of moraine ridges in this region. There, we reconstruct an ice sheet margin that remained more extensive than present for 2000 years longer than further south, until ∼ 7–6.5 kyr BP, as indicated by deglacial event ages by Bennike and Weidick (2001) and Larsen et al. (2018, 2020).

Figure 8 Time slice maps of the PaleoGrIS 1.0 isochrone reconstruction in Nares Strait region. The reconstructed ice sheet areal extent for each given time slice is displayed in each subpanel as a translucent white polygon, while the underlying opaque white layer is the present-day ice cover from BedMachine v4. The reconstructed ice sheet margins are highlighted by our isochrone polylines, which feature different colour schemes relating to our four isochrone confidence levels (see Sect. 2.3.5). We merge our isochrones with the reconstruction of former ice extent over Ellesmere Island (Innuitian Ice Sheet) from Dalton et al. (2020). While no modifications were applied to the margin extent of that dataset, the timing of Dalton et al. (2020) isochrones was in some cases modified (inferior to 1 kyr at most)to match our reconstruction. Bottom panels also feature present-day ice sheet thickness (Morlighem et al., 2017) and surface velocity data (Joughin et al., 2018). Topography and bathymetry are from the GEBCO 2022 release.

3.2.3 Ice sheet retreat in central east Greenland and ice margin evolution in the Scoresby Sund and Kangerlussuaq regions

The central east (CE) region of the Greenland Ice Sheet is characterised by two major ice drainage basins. The northernmost one is composed of ice streams and outlet glaciers flowing eastwards and into fjord systems that merge to form Scoresby Sund. The southernmost basin comprises three major ice streams flowing southeastwards, namely the Kangerlussuaq, Christian IV, and Hutchinson plateau glaciers (Rignot and Mouginot, 2012; Bjørk et al., 2015). We here describe the general pattern and timing of retreat in these two key regions (Fig. 10).

The Scoresby Sund region

At ∼ 14–13 kyr BP, we reconstruct a grounded ice sheet margin positioned towards the inner continental shelf, near the mouth of Scoresby Sund and around the location of the underwater Cape Brewster moraine complex (Dowdeswell et al., 1994; Fig. 10). The timing of this specific extent is highly uncertain and exclusively based on a TCN-derived chronology from Cape Brewster by Håkansson et al. (2007a). For this time slice, we reconstruct the ice sheet as connected with ice caps from the Liverpool Land peninsula (Fig. 10). However, our reconstruction tentatively suggests that a disconnection between the two occurred shortly after, towards ∼ 13–12 kyr BP. By ∼ 12–11.5 kyr BP, deglacial radiocarbon ages (Bennike et al., 1999) indicate that the Scoresby Sund outlet had retreated within the fjord, with its northern lateral margin resting against the southwestern slopes of Jameson Land. Following further glacier retreat and thinning, we reconstruct the disconnection of ice flowing southeastwards within Hall Bredning fjord and ice flowing northeastwards within Gaasefjord to occur between ∼ 11.5 and ∼ 10 kyr BP, based on radiocarbon deglacial chronologies by Marienfeld (1990), Ingólfsson et al. (1994), and Hansen (2001). We propose that progressive westward retreat and separation into three outlet glaciers retreating into Nordvestfjord, Harefjord, and Fonfjord, occurred between ∼ 9.5 and ∼ 8 kyr BP. After ∼ 8 kyr BP, we estimate the grounded ice extent was similar to or less than the present in this region (Funder, 1978).

The Kangerlussuaq glacier region

During the last full glacial configuration of the ice sheet, outlets of the southernmost basin in CE Greenland converged into one major lobe, commonly referred to as the Kangerlussuaq outlet glacier (Dowdeswell et al., 2010). This major outlet flowed southeastwards and is thought to have likely reached the outer continental shelf during last full glacial extent (Mienert et al., 1992), although this remains uncertain. The deglacial pattern and timing of retreat of the Kangerlussuaq outlet has been studied by a relatively large number of offshore coring (e.g. Williams, 1993, 1995; Andrews et al., 1996; Smith and Licht, 2000) and onshore (e.g. Dyke et al., 2014) investigations. Based on findings from these empirical studies, we reconstruct a grounded ice margin reaching a mid-shelf position towards ∼ 14–13 kyr BP in this area (Fig. 10). Between ∼ 14 and ∼ 10 kyr BP, our reconstruction suggests the Kangerlussuaq outlet glacier retreated in a relatively slow and monotonic manner (∼ 10 m yr⁻¹). Radiocarbon-dated marine sediment records (Smith and Licht, 2000) suggest that faster retreat between ∼ 10 and ∼ 9 kyr BP (∼ 70–40 m yr⁻¹) caused the outlet glacier to reach the Kangerlussuaq Fjord mouth by ∼ 9 kyr BP. We estimate later retreat within the fjord to have been relatively quick (∼ 60 m yr⁻¹), and a similar or more-retreated than present extent was likely reached before ∼ 8 kyr BP (Dyke et al., 2014).

Figure 9 Time slice maps of the PaleoGrIS 1.0 isochrone reconstruction in the Nioghalvfjerdsfjorden and Jøkelbugten regions. The reconstructed ice sheet areal extent for each given time slice is displayed in each subpanel as a translucent white polygon, while the underlying opaque white layer is the present-day ice cover from BedMachine v4. The reconstructed ice sheet margins are highlighted by our isochrone polylines, which feature different colour schemes (see Fig. 8 for key) relating to our four isochrone confidence levels (see Sect. 2.3.5). Bottom panels also feature present-day ice sheet thickness (Morlighem et al., 2017) and surface velocity data (Joughin et al., 2018). Topography and bathymetry are from the GEBCO 2022 release. For a more complete figure legend, the reader is advised to consult Fig. 8.

3.2.4 Ice sheet retreat in the southeast region and ice margin evolution of the Sermilik outlet glacier

To the south of the Kangerlussuaq Fjord and Kangerlussuaq Trough, a 300 km long stretch of the Greenland coastline is characterised by rugged, steep, and high-elevation mountains displaying numerous summits reaching 2000 m a.s.l., with some acting as contemporary nunataks. These high coastal mountains act as orographic barriers causing snow accumulations and resulting in the contemporary ice sheet reaching the shore along the entire coastline. This results in a lack of settlements, difficult access conditions, and thus a lack of paleo glaciological field investigations. Moreover, in this region, the offshore continental shelf does not yet feature published geochronological constraints on deglacial grounded ice margin retreat, as far as our compilation suggests. Therefore, our ice sheet margin reconstruction along this coastline displays only very low-confidence isochrones, and retreat is crudely assumed to be monotonic in nature. At a distance of 300 km further south, the next location presenting empirical constraints on past ice sheet margin evolution is the large Sermilik fjord, in which the rapidly flowing Helheim (up to 8000 m yr⁻¹; Joughin et al., 2018), Apuseerajik, and Nigertiip Apusiia tidewater glaciers terminate. The Sermilik fjord is bordered on its eastern side by the large Ammassalik Island and other peninsulas which feature steep mountains reaching > 1000 m a.s.l. (above sea level). These high topographies act as orographic barriers to moisture supply and enable sustaining numerous ice caps and mountain glaciers.

Today, the deep Sermilik fjord is thus characterised by its main tidewater glacier fronts (e.g. Helheim) terminating more than ∼ 100 km inside the fjord, while the fjord's terrestrial flanks remain heavily glaciated. We attempted to mimic this characteristic ice configuration when reconstructing the deglacial retreat of the ice sheet margin in this region. Based on extrapolating chronological constraints from the Kangerlussuaq Trough, we map a very low-confidence outermost isochrone (14–13 kyr BP) located towards the inner- to mid-continental shelf offshore the Sermilik fjord. Soon after, by 13–12 kyr BP, we locally reconstruct an ice margin that had quickly retreated towards the mouth of the Sermilik fjord and onto the shore of local coastal mountains. This relatively early retreat scenario is constrained by TCN exposure ages from the southwestern tip of Ammassalik Island (Hughes et al., 2012), located east of Sermilik fjord, and by exposure ages from the end of the Torqulertivit Imiat valley (Roberts et al., 2008), located near the western edge of the main fjord mouth. Since deglacial radiocarbon ages from the Ammassalik Fjord further east (Long et al., 2008) give similar ages to TCN exposure dates from rock surfaces located towards the mid- to inner-Sermilik fjord sides (Hughes et al., 2012), we reconstruct a Sermilik fjord outlet front that had retreated ∼ 30 km north by 11.5–11 kyr BP, while ice sheet margins remained extensive, i.e. near the coast and towards the Johan Petersen Fjord mouth, on the western side of Sermilik fjord. The Sermilik outlet thus appears to have experienced a potent state of negative mass balance during the Younger Dryas to early-Holocene rapid warming transition (GISP2 data; Alley, 2000), while ice masses occupying the fjord-side mountains either remained stable or retreated less and slower. Approximately 25 km further north, exposure ages from the Amanga island (Hughes et al., 2012), situated where the fjord splits into three, suggest the Sermilik outlet had rapidly retreated further into the inner fjord and started to split into the three main tidewater glaciers existing today (Helheim, Apuseerajik, and Nigertiip Apusiia) by approximately 11–10.5 kyr BP. Although uncertain, this rapid retreat prompts us to reconstruct front margins for these three outlets that were as or more retreated than present after 10–9.5 kyr BP, which is an early timing for this event relative to other Greenland regions, we find.

Figure 10 Time slice maps of the PaleoGrIS 1.0 isochrone reconstruction in the Scoresby Sund and Kangerlussuaq regions. The reconstructed ice sheet areal extent for each given time slice is displayed in each subpanel as a translucent white polygon, while the underlying opaque white layer is the present-day ice cover from BedMachine v4. The reconstructed ice sheet margins are highlighted by our isochrone polylines, which feature different colour schemes (see Fig. 8 for key) relating to our four isochrone confidence levels (see Sect. 2.3.5). Bottom panels also feature present-day ice sheet thickness (Morlighem et al., 2017) and surface velocity data (Joughin et al., 2018). Topography and bathymetry are from the GEBCO 2022 release. For a more complete figure legend, the reader is advised to consult Fig. 8.

3.2.5 Ice sheet retreat in southernmost Greenland

The southernmost region of Greenland is characterised by a relatively narrow (< 70 km) and shallow (typically < 200 m) continental shelf, and features high-elevation coastal mountains with several summits reaching 2000 m a.s.l. The present-day ice sheet margin is generally located within 100 km of the outer coast in this region. Compared with other Greenland regions, little ice build-up was required for the ice sheet to reach the continental shelf edge during the last full glacial extent, and this scenario is considered highly likely by previous investigations (e.g. Funder et al., 2011; Andersen et al., 2020; Levy et al., 2020).

The relatively short distance to the continental shelf combined with high local Bølling–Allerød warming and the proximity to the Irminger Current enabling warm water incursion are thought to have caused relatively early deglaciation in southernmost Greenland (Levy et al., 2020). In fact, the well-studied N14 lake record, located on the Kitsissut Islands (59.98° N, −45.18° W), features the oldest radiocarbon deglacial isolation age (∼ 13.6 kyr BP) in Greenland (Bennike et al., 2002; Puleo et al., 2022). Based on these studies, southernmost Greenland is the only region in which we map the oldest isochrone (14–13 kyr BP) onshore, towards the present-day outer coastline (Fig. 11). TCN event ages from Lindenows Fjord (Levy et al., 2020) suggest that the series of fjords located east of the Qaqortoq ice cap started deglaciating at around ∼ 12 kyr BP. Although a lack of evidence makes isochrone mapping in adjacent northern fjords highly uncertain, our reconstruction displays outlet glaciers retreating and reaching present-day configurations earlier than most Greenland regions by approximately 10 kyr BP. Therefore, the ice sheet margin was as retreated as (or more retreated) than the present for most of the Holocene in this region.

To the west of Uummannarsuaq (also known as Cape Farewell), in the region of the Nanortalik and Narsarsuaq settlements, the timing of ice sheet retreat (from Qaqortoq bay) is better constrained by empirical evidence. The majority of the local ice basin is currently drained by three fast-flowing outlet glaciers (maximum surface velocities > 1000 m yr⁻¹), namely the Eqalorutsit Kangilliit Sermiat, Eqalorutsit Killiit Sermiat, and Qooqqup Sermia glaciers (Bjørk et al., 2015; Joughin et al., 2018). TCN exposure ages from Winsor et al. (2015) indicate these outlet glacier fronts had retreated ∼ 20 km into local fjords by ∼ 13–12 kyr BP, while TCN ages from Nelson et al. (2014) suggest retreat over the next 2000 years caused the outlet glacier fronts to be within ∼ 15 km of the present day before ∼ 10 kyr BP (Fig. 11). This scenario of retreat is also consistent with a deglacial radiocarbon age produced by Weidick (1975) from the Bredefjord (the Narssaq site). To the northwest, between Qaqortoq bay and the large piedmont Sioqqap Sermia outlet glacier, our compilation features no geochronological constraints over a ∼ 200 km long stretch of ice-free coastal land. We have identified this region as a possible target for future field investigations aiming to reconstruct Holocene ice margin fluctuations. Along this region, low-confidence isochrones were thus mapped solely by linking mapped ice marginal landforms and by extrapolating retreat timings from adjacent areas. Further north, the next site to display chronological constraints on ice sheet retreat is Kuannersooq fjord (61.98° N) to the east of the Paamiut settlement, with TCN exposure ages from Winsor et al. (2015) dating an ice sheet margin retreating to the fjord mouth by 12–11.5 kyr BP. Furthermore, in this valley, Carlson et al. (2014) sampled boulder and bedrock surfaces shorter than 2 km from the present-day ice sheet margin. The resulting TCN exposure ages suggest the ice sheet had retreated towards its present-day extent by ∼ 10 kyr BP in this sector.

Figure 11 Time slice maps of the PaleoGrIS 1.0 isochrone reconstruction in southernmost Greenland. The reconstructed ice sheet areal extent for each given time slice is displayed in each subpanel as a translucent white polygon, while the underlying opaque white layer is the present-day ice cover from BedMachine v4. The reconstructed ice sheet margins are highlighted by our isochrone polylines, which feature different colour schemes (see Fig. 8 for key) relating to our four isochrone confidence levels (see Sect. 2.3.5). Bottom panels also feature present-day ice sheet thickness (Morlighem et al., 2017) and surface velocity data (Joughin et al., 2018). Topography and bathymetry are from the GEBCO 2022 release. For a more complete figure legend, the reader is advised to consult Fig. 8.

3.2.6 Ice sheet retreat in southwest Greenland and ice margin evolution in the Nuuk region

The ice-free coastal region to the east and south of Nuuk presents a complex topography with high (1500–2000 m a.s.l.) and steep-sided mountains dissected by several sinuous and deep fjords forming an archipelago. Upstream of these fjords, the majority of the ice discharge is captured by four fast-flowing outlet glaciers (maximum surface velocities 1000–2000 m yr⁻¹). From south to north, these are the Sermeq, Kangiata Nunaata Sermia, Akullersuup Sermia, and Narsap Sermia glaciers (Bjørk et al., 2015). Due to easier accessibility, several investigations studying the former ice sheet evolution have taken place in this area, and the Late Glacial and early- to mid-Holocene retreat of the ice sheet's margin is relatively well constrained. Further north, however, to the east of the Atammik settlement, a large ice-free region displaying flatter topographies with fewer overdeepenings features fewer geomorphological and geochronological constraints on ice margin retreat.

In the well-studied valleys, numerous TCN and radiocarbon dates from four parallel fjords and from three distinct investigations (Weidick, 1976; Winsor et al., 2015; Larsen et al., 2014) show consistent results. These constraints enable us to map, with good levels of confidence, an ice sheet margin positioned towards the mouths of the Nuup Kangerlua, Kangerdluarsunguak, and Sermilik fjords, while positioned towards the middle of the Ameralik fjord, at ∼ 11–10.5 kyr BP (Fig. 12). Older isochrones (14–11 kyr BP) are less well constrained and are here reconstructed following a low-confidence and monotonic retreat pattern from a mid-shelf position at 14–13 kyr BP. Several extensive investigations sampling along the fjords and towards the modern ice margin have shown that after ∼ 11–10.5 kyr BP, local outlets flowing into the Nuuk Fjord system likely experienced very rapid retreat, with their fronts reaching the inner fjords by ∼ 10 kyr BP (Weidick, 1972; Larsen et al., 2014; Young et al., 2021) (Fig. 12). According to this retreat scenario, local fjords were deglaciated in fewer than 1000 years, which represents a retreat rate of 207 ± 69 m yr⁻¹ for the Kangiata Nunaata Sermia outlet glacier during that time. Out of the 72 Greenland outlet glaciers sampled for retreat rate analysis, this rate is the sixth-highest maximum retreat rate reached over the reconstructed period. Therefore, this rapid retreat of the ice sheet margin in the Nuuk region between ∼ 11 and ∼ 10 kyr BP was likely a high-magnitude event at the ice sheet scale.

Following this, numerous multi-nuclide TCN exposure ages from near the contemporary margins of the Kangiata Nunaata Sermia, Akullersuup Sermia, and Narsap Sermia glaciers indicate local tidewater glaciers had reached an extent similar to the present by ∼ 9 kyr BP (Larsen et al., 2014; Young et al., 2021). Further north, towards a plateau (> 600 m a.s.l.) lying between the Nuup Qinngua fjord and the Saqqap Sermia outlet glacier, younger deglacial radiocarbon ages (∼ 7.5 kyr BP; Levy et al., 2017) suggest that the land-terminating ice sheet margin retreated at a slower pace. Although this hypothesis requires further evidence for validation, our reconstruction uses these data to suggest an ice margin remaining more extensive than the present until our youngest isochrone, i.e. 7–6.5 kyr BP (Fig. 12).

Figure 12 Time slice maps of the PaleoGrIS 1.0 isochrone reconstruction in the Nuuk region. The reconstructed ice sheet areal extent for each given time slice is displayed in each subpanel as a translucent white polygon, while the underlying opaque white layer is the present-day ice cover from BedMachine v4. The reconstructed ice sheet margins are highlighted by our isochrone polylines, which feature different colour schemes (see Fig. 8 for key) relating to our four isochrone confidence levels (see Sect. 2.3.5). Bottom panels also feature present-day ice sheet thickness (Morlighem et al., 2017) and surface velocity data (Joughin et al., 2018). Topography and bathymetry are from the GEBCO 2022 release. For a more complete figure legend, the reader is advised to consult Fig. 8.

3.2.7 Ice sheet retreat in central west Greenland and ice margin evolution in the Sisimiut, Disko Bay, and the Uummannaq fjord regions

Sisimiut region

The coastal region here described as belonging to central west (CW) Greenland stretches 600 km northward from the large Maniitsoq Ice Cap (66°13^′ N, 52°11^′ W) to the Uummannaq Fjord region. From the Maniitsoq Ice Cap to the southern edge of Disko Bay (Qeqertarsuup tunua), the ice-free coast is characterised by a relatively flat bedrock plateau dissected by deep fjords and extending ∼ 150 km from the outer coast to the contemporary ice margin (Fig. 13). This region, home to the towns of Sisimiut and Kangerlussuaq, has been the subject of numerous investigations reconstructing the evolution of the ice sheet margin during deglaciation (e.g. Ten Brink and Weidick, 1974; Eisner et al., 1995; Bennike, 2000; Roberts et al., 2009; Rinterknecht et al., 2009; Storms et al., 2012; Winsor et al., 2015; Kelley et al., 2013, 2015; Lesnek and Briner, 2018; Briner et al., 2020; Young et al., 2020). Moreover, ice marginal landforms are well preserved and highly concentrated in this area. As a result, the Late Glacial and early- to mid-Holocene retreat of the ice sheet margin is locally well constrained. Based on this high density of geological evidence, we reconstruct, with reasonable levels of confidence, an ice sheet margin that had retreated to within 30 km (offshore or onshore depending on the valley) of the present-day outer coast by ∼ 12–11.5 kyr BP (Fig. 13). After that time, the eastward ice sheet margin retreat appears to have been relatively monotonic in this region. Local outlet glaciers sampled, i.e. the south Russell, Isunnguata Sermia, Inuppaat Quuat, and Akuliarutsip Sermia glaciers (Bjørk et al., 2015), seem to have retreated at intermediate but steady speeds varying between 25 and 30 m yr⁻¹ once averaged. By ∼ 9 kyr BP, we reconstruct an ice sheet margin that was near (within 20 km) the inner end of local fjords (e.g. Kangerlussuaq, Nassuttooq, Qasigiarssuit). Between ∼ 9 and ∼ 7 kyr BP, empirical evidence suggests that the ice sheet started retreating slower and was located closer than 30 km from the present-day margin for ∼ 3 kyr in most locations throughout this region. This slowdown in retreat appears to have been coeval with the ice sheet margin becoming fully land-terminating in most valleys and is thus a potential consequence of the removal of calving-related ablation. The local reconstructed retreat pattern is such that the present-day ice sheet extent was not reached until ∼ 7–6.5 kyr BP.

Disko Bay

Further north, at 68.7° N, lies Disko Bay (Qeqertarsuup tunua), formerly host to a large streaming outlet glacier (also referred as the Jakobshavn Isbræ ice stream) that flowed towards Baffin Bay. This ice stream is thought to have reached the continental shelf break during the last full glacial extent and likely remained near such an extent until ∼ 14 kyr BP (Ó Cofaigh et al., 2013; Rinterknecht et al., 2014). For our outermost isochrone (14–13 kyr BP), we thus reconstruct an ice sheet margin located towards the middle-to-outer continental shelf, mostly based on data from Ó Cofaigh et al. (2013) (Fig. 13). Data from McCarthy (2011) and Rinterknecht et al. (2014) enable us to tentatively draw an ice sheet margin that had retreated towards the inner continental shelf, with ice thinning causing its lateral margin to rest against the southern coastal mountain slopes of Disko Island.

Based on offshore data by Lloyd et al. (2005) and TCN exposure ages by Kelley et al. (2013, 2015), we estimate that the further retreat into Disko Bay had caused the ice margin to be positioned 30 km offshore the mouth of the Ilulissat Icefjord, and towards the Akunnaaq and Qasigiannguit settlements further south, by 10.5–10 kyr BP. The following phase of retreat was likely characterised by local outlet glaciers (i.e. Sermeq Kujalleq, Eqip Sermia, Saqqarliup Sermia, Sermeq Avannarleq, and Akuliarutsip Sermia) retreating into fjords to the east of Disko Bay. These proglacial regions have been studied extensively and are characterised by dense mapping and numerous geochronological constraints (e.g. Donner and Jungne, 1975; Long et al., 1999, 2006; Weidick and Bennike, 2007; Briner et al., 2010; Young et al., 2013b; Carlson et al., 2014; Cronauer et al., 2016; Balter-Kennedy et al., 2021). This enables a high density of mid- to high-confidence isochrones to be drawn. To the south of Ilulissat Icefjord, empirical data suggest the ice sheet margin evolution was characterised by slow and steady retreat between ∼ 10 and ∼ 7 kyr BP. At a maximum, local outlet glacier fronts (e.g. Akuliarutsip Sermia) had retreated ∼ 20 km during that time interval, which suggest retreat rates likely below 10 m yr⁻¹. Based on these observations, we reconstruct an ice margin that reached the present-day extent by ∼ 7.5–6.5 kyr BP in this area. For the Sermeq Kujalleq outlet glacier (also known as Jakobshavn Isbræ), extensive empirical datasets from both sides of the Ilulissat Icefjord enable us to constrain, with high levels of confidence, the position of the outlet's lateral margins through time (Briner et al., 2010; Young et al., 2013). These data indicate that between ∼ 10 and 8 kyr BP, the glacier calving front was either stable, or slowly retreating, and remained near the fjord's mouth, with the glacier's right lateral margins terminating on the Ilulissat peninsula (Fig. 13). After ∼ 8 kyr BP, the calving front likely retreated to a mid-fjord position, with lateral ice margins retreating to within ∼ 5 km of contemporary ones by ∼ 7 kyr BP. We estimate the Sermeq Kujalleq had retreated to a near-present-day extent between ∼ 7 and ∼ 6.5 kyr BP.

Figure 13 Time slice maps of the PaleoGrIS 1.0 isochrone reconstruction in the Sisimiut and Disko Bay region. The reconstructed ice sheet areal extent for each given time slice is displayed in each subpanel as a translucent white polygon, while the underlying opaque white layer is the present-day ice cover from BedMachine v4. The reconstructed ice sheet margins are highlighted by our isochrone polylines, which feature different colour schemes (see Fig. 8 for key) relating to our four isochrone confidence levels (see Sect. 2.3.5). Bottom panels also feature present-day ice sheet thickness (Morlighem et al., 2017) and surface velocity data (Joughin et al., 2018). Topography and bathymetry are from the GEBCO 2022 release. For a more complete figure legend, the reader is advised to consult Fig. 8.

Uummannaq Fjord

Located ∼ 250 km north from the former Jakobshavn Isbræ ice stream, a major outlet glacier, known as the Uummannaq ice stream (Lane et al., 2014), also formerly advanced towards Baffin Bay with rapid flow, causing the formation of streamlined subglacial bedforms (Ó Cofaigh et al., 2013). Radiocarbon ages from a marine sediment core (VC45) suggest that during the last full glacial configuration, this outlet glacier likely reached the continental shelf break and remained in this position until ∼ 15–14.5 kyr BP (Ó Cofaigh et al., 2013). Given this maximum-limiting constraint, we draw our outermost isochrone (14–13 kyr BP – very low confidence) in this area towards the middle-to-outer continental shelf and map the Uummannaq ice stream margin as a protruding lobe (Fig. 14). A radiocarbon-derived event age by Bennike et al. (1994) indicates the ice sheet margin had retreated from the lower slopes of Hareoen Island prior to ∼ 12.2 kyr BP. We thus tentatively reconstruct south and north lateral margins of the Uummannaq outlet positioned towards the outer contemporary coast between ∼ 13 and ∼ 12 kyr BP, while mapping a protruding lobe along the central bathymetric trough with a margin reaching a mid-shelf position (Fig. 14).

Further east, a TCN-derived event age produced by Roberts et al. (2013) suggest the highlands of Illorsuit island were becoming ice-free by ∼ 11.5–11 kyr BP. For this isochrone, we therefore reconstruct an ice sheet margin that separates into two distinct glaciers, with a northern glacier retreating northeastwards into Karrat Fjord and a southern glacier retreating southeastwards into Uummannaq Fjord (Fig. 14). Further retreat led to divisions into 12 distinct outlets retreating in a dendritic series of narrow fjords that are now host to the calving fronts of several fast-flowing tidewater glaciers, such as Salliarutsip Sermia (72.0° N), Umiammakku Sermiat (71.7° N), Kangilliup Sermia (Rink Isbræ; 71.7° N), or Sermeq Kujalleq (Store Gletscher; 70.4° N). In agreement with TCN-derived event ages by Lane et al. (2014), our reconstruction suggests that by 10–9.5 kyr BP, the Umiammakku Sermiat and Kangilliup Sermia glacier fronts had retreated towards the fjord mouths, near the western edge of Karrat Island, ∼ 50–30 km from their present-day fronts (Fig. 14). Therefore, in the northern sector of the Uummannaq Fjord system (Karrat area), we estimate that rapid retreat of the ice sheet margin was experienced between ∼ 11.5 and ∼ 9.5 kyr BP.

Informed by data from the Qarajaq Icefjord further south (Simonarson, 1981; Roberts et al., 2013), our reconstruction also suggests a rapid retreat of the ice margin in the southern Uummannaq Fjord system, between ∼ 11 and ∼ 10 kyr BP. However, after ∼ 10–9.5 kyr BP, northern outlet glacier fronts (Karrat area) appear to have remained relatively stable, with fronts staying more extensive than present until at least 5 kyr BP and hence for significantly longer than most other Greenland regions (Lane et al., 2014). This contrasts with the ice sheet margin evolution in the southern Uummannaq sector; 60 km further south, TCN-derived event ages from Philipps et al. (2017) suggest that the lateral margins of the Perlerfiup Sermia (71.0° N) glacier were within 5 km of the present-day margin by 10.5–10 kyr BP (Fig. 14). This dataset further indicates that after ∼ 9.5 kyr BP, this outlet was either as retreated as (or more retreated than) the present-day margin. This southern sector is characterised by lower coastal mountains, more gradual slopes, wider fjords, and fewer mountain glaciers today than towards the northern Uummannaq Fjord system. This setting may have caused a smaller ice mass contribution from periphery mountain glaciers acting as tributaries than further north. At a distance of 70 km further south, our reconstruction suggests the Sermeq Kujalleq (Store) glacier was more extensive than the present until ∼ 9–8 kyr BP (Roberts et al., 2013). Overall, empirical data suggest the early- to mid-Holocene margin response of individual outlet glaciers of the Uummannaq Fjord system was complex and heterogenous, and we try to capture this variability in our isochrone reconstruction.

Figure 14 Time slice maps of the PaleoGrIS 1.0 isochrone reconstruction in the Uummannaq Fjord region. The reconstructed ice sheet areal extent for each given time slice is displayed in each subpanel as a translucent white polygon, while the underlying opaque white layer is the present-day ice cover from BedMachine v4. The reconstructed ice sheet margins are highlighted by our isochrone polylines, which feature different colour schemes (see Fig. 8 for key) relating to our four isochrone confidence levels (see Sect. 2.3.5). Bottom panels also feature present-day ice sheet thickness (Morlighem et al., 2017) and surface velocity data (Joughin et al., 2018). Topography and bathymetry are from the GEBCO 2022 release. Note that between 9–8.5 and 7–6.5 kyr BP, in this region, sections of the former ice sheet margin may have retreated beyond its present-day extent (opaque white polygon). For a more complete figure legend, the reader is advised to consult Fig. 8.

3.2.8 Ice sheet retreat in northwest Greenland and ice margin evolution in the Upernavik and Aappilattup Ikera region

To the north of the Uummannaq Fjord region, the ∼ 800 km coastline of Greenland that stretches up to Steensby Land (77° N) is commonly described as belonging to the NW Greenland region. This coastline is characterised by little ice-free land as the present-day ice sheet margin reaches the shore along most of the coastline, and this region features relatively low-lying coastal topographies. In this sector, according to our compilation, few studies constrain the timing and positions of the former Greenland Ice Sheet margin. The exception to this is the region of Upernavik, however, where two studies (Briner et al., 2013; Corbett et al., 2013) provide numerous TCN-exposure-age-derived time constraints on the deglacial evolution of the Sermeq (Upernavik Isstrøm) tidewater glacier, which drains a major catchment of the NW Greenland Ice Sheet. There, we use the deglacial TCN exposure ages from westernmost Islands located 15 km offshore the Upernavik settlement (Corbett et al., 2013) to map a retreated outermost isochrone (14–13 kyr BP) located towards the inner continental shelf, approximately 30 km west of the mouth of Aappilattup Ikera (the Upernavik Isfjord). Younger TCN exposure ages and dipstick models obtained from high topographies that lie further east along the Upernavik archipelago seem to indicate a relatively monotonic and slow (∼ 15 m yr⁻¹) ice sheet margin retreat pattern between 14–13 and 11.5–11 kyr BP. We thus reconstruct the Sermeq outlet glacier front to rest towards the outer- to mid-Aappilattup Ikera fjord by 11.5–11 kyr BP. Further inland, TCN exposure ages from an island adjacent to the historical ice margin (Briner et al., 2013) indicate faster retreat likely caused the ice sheet margin to recede towards the inner fjord by 10–9.5 kyr BP. Moreover, according to three deglacial radiocarbon ages from lake cores sampled within 5 km of the contemporary ice front (Briner et al., 2013), we reconstruct a local ice sheet margin that reached an extent that is similar to the present-day shortly after ∼ 9.5–9 kyr BP. Along the coastline located to the north of the Upernavik sector, our reconstruction mostly features very low-confidence isochrones crudely depicting a retreat pattern assumed to be monotonic, due to a lack of geochronological and land-based geomorphological evidence.

Figure 15 Ice-sheet-scale map of the PaleoGrIS 1.0 isochrone reconstruction divided in two distinct panels (north Greenland and south Greenland) for readability purposes. The individual isochrone lines mapped as part of this study, and their four respective confidence levels, are too detailed to be visualised at this scale. Instead, we show a colour map of the areas located between isochrones. Consequently, each colour-coded mask (or isochrone buffer) highlights the former location of the Greenland Ice Sheet margins between the timings of the outer and inner isochrones that it is delimited by. In order to visualise the full details of the isochrone reconstruction, along with the geomorphological and geochronological data compiled, the reader is advised to download and zoom into the PaleoGrIS 1.0 poster (A0 format) provided in the online database. Alternatively, the reader can download the PaleoGrIS 1.0 shapefile database and visualise the reconstruction at any scale in geographic information system (GIS) software.

3.3 Analysis of changes in areal extent of the Greenland Ice Sheet

3.3.1 Ice-sheet-wide areal extent evolution

Our reconstruction at 14–13 kyr BP suggests that the Greenland Ice Sheet had an approximate areal extent of 2.61 million km² (Figs. 15, 16) and lost about one-third of this areal extent (0.89 million km²) as it reduced in size to its present-day extent (1.71 million km²; Morlighem et al., 2017; Fig. 16). Between 14–13 and 9–8.5 kyr BP, our reconstruction produces a near-constant rate of ice sheet areal extent loss of 170 ± 27 km² yr. Later, between 9–8.5 and 7–6.5 kyr BP, our reconstruction suggests a rate of areal extent loss that progressively decreased through time (Fig. 16). However, as our areal extent estimates are only maximum-limiting between 10–9.5 and 7–6.5 kyr BP (see Sect. 2.4), this potential slow down remains hypothetical and could be a consequence of our limited knowledge regarding how far the ice sheet margin retreated behind the present-day position in response to the Holocene thermal maximum.

At 7–6.5 kyr BP, we estimate the ice sheet areal extent was 1.71 million km² or smaller. This suggests the Greenland Ice Sheet extent was similar to or smaller than present by this time. Although highly uncertain, our literature-based maximum and minimum full glacial extent scenarios (see Sect. 2.6) suggest that the ice sheet areal extent was likely between ∼ 3.13 and ∼ 2.94 million km² during the last full glacial configuration (Figs. 5, 16). We thus estimate that before the start of the Holocene (∼ 11.7 kyr BP), deglaciation had caused the ice sheet to lose between ∼ 26 % and ∼ 21 % of its full glacial extent. Therefore, we find that between ∼ 57 % and ∼ 51 % of the post-glacial areal extent loss occurred before the onset of the Holocene, while a significant proportion (∼ 49 %–43 %) still occurred during the early- to mid-Holocene interval.

Figure 16 Empirical estimation of the Greenland Ice Sheet areal extent evolution between ∼ 18 kyr BP and the present. The ice sheet extent during last full glacial configuration (red box) is estimated here by bracketing the maximum (at continental shelf break everywhere) and minimum extent scenarios, reconstructed after consulting the relevant literature, in an attempt to highlight the uncertainties in several debated regions (more details in Sect. 2.6). The hypothetical timing of such full glacial extent is crudely associated with maximum atmospheric and oceanic cooling over Greenland during Heinrich Stadial 1 (Buizert et al., 2018), but we stress that the timing of maximum extent is unknown and could have been reached before 18 kyr BP and until Bølling–Allerød warming started after ∼ 16 kyr BP, as highlighted by our conservative range (red box). Ice sheet areal extent estimates from the PaleoGrIS 1.0 isochrone reconstruction are shown with green circles (actual extent estimates) and orange triangles (maximum-limiting extent pointing upwards). This differentiation is due to our isochrones not covering the full ice sheet perimeter during the ∼ 10–6.5 kyr BP interval when some sections of the ice sheet margin are estimated to have been as retreated as (or more retreated than) the present margins.

3.3.2 Region-specific areal extent change

After dividing the ice sheet into seven major drainage basins (after Rignot and Mouginot, 2012; see Sect. 4), we analysed regional patterns of areal extent change through time. While all regions lost areal extent between 14–13 and 7–6.5 kyr BP, we find the timing, rate, and magnitude of ice sheet retreat differs substantially between regions (Fig. 17). By 7–6.5 kyr BP, the NO, SW, and CE regions had lost at least 40 % of their 14–13 kyr BP areal extent. The same figure is nearly twice as low (∼ 25 %) for the NW and SE regions. During that time, the CW and NE regions lost at least ∼ 30 % and ∼ 35 % of their areal extent, respectively. Furthermore, our reconstruction suggests that the rate of areal extent loss, and its evolution through time, varied significantly between regions. For instance, the SW region lost areal extent relatively quickly from the period onset (∼ 12 % kyr⁻¹–8 % kyr⁻¹) and reached an extent that is similar to the present earlier than other regions by 10–9.5 kyr BP. The NO region also lost areal extent at a relatively quick rate (∼ 8 % kyr⁻¹–6 % kyr⁻¹), but such rapid loss started later (after ∼ 12–11.5 kyr BP), and an areal extent that is similar to the present was not reached until ∼ 7–6.5 kyr BP (Fig. 17). In the NE region, on the other hand, empirical data suggest that an extent that is similar to the present was reached by approximately 8–7.5 kyr BP. Before that time, the rate of the percentage areal extent loss was relatively constant and at ∼ 7 % kyr⁻¹–5 % kyr⁻¹. In the CE region, the ice sheet areal extent loss was relatively slow (∼ 7 % kyr⁻¹-5 % kyr⁻¹) until ∼ 10–9.5 kyr BP, prior to accelerating substantially (∼ 18 % kyr⁻¹–9 % kyr⁻¹) until ∼ 8.5–8 kyr BP, by which time it had reached an extent that is similar to the present. In our reconstruction, the CE regions thus appears to have experienced a rapid and delayed collapse of a significant proportion of its total areal extent between ∼ 10 and ∼ 8 kyr BP (Fig. 17). The SE region reached an extent that is similar to the present relatively early, by 9.5–9 kyr BP, and experienced relatively slow areal extent loss prior to that (∼ 7 % kyr⁻¹–5 % kyr⁻¹). The rate of areal extent loss in the CW region decreased progressively through time, and an extent that is similar to the present was reached late in this region, after ∼ 8.5 kyr BP. The NW region follows a similar pattern but reached near-present-day extent slightly earlier, prior to ∼ 9 kyr BP (Fig. 17). Such comparisons should be treated with caution, as they are sometimes derived from low- or very low-confidence isochrones. However, the data generally suggest a latitudinal signal in the different regional retreat patterns by which southernmost regions reached a near-present-day margin position earlier than northern regions, which experienced a more delayed deglaciation.

Figure 17 The PaleoGrIS 1.0 reconstruction of percentage ice sheet areal extent loss between ∼ 14 kyr BP (0 %) and ∼ 6.5 kyr BP for each of the seven ice sheet regions spatially divided, as displayed in the inset map, after Rignot and Mouginot (2012). The colour code of the individual time series matches the inset map polygons. Importantly, we use a different symbol to highlight whether the percentage of areal extent loss, for a given region, is an absolute areal extent estimate (straight line) or a minimum-limiting areal extent loss estimate (dashed line with downward-pointing arrows). The switch from straight to dashed lines thus represents, for each region, the youngest isochrone to feature an extent that is more extensive than the present ice sheet extent along the full region-specific ice margin.

3.4 Outlet glacier retreat rates analysis

Transects were drawn across reconstructed isochrones for 72 outlet glaciers sampled around Greenland (further methods in Sect. 2.5) (Fig. 18). Measuring distances along these transects allows us to estimate that overall retreat rates of glacier fronts (or of the grounding line for marine margins) ranged between 72 and 8 m yr⁻¹ during the ∼ 14–6.5 kyr BP time period (mean and 1σ SD 37 ± 16 m yr⁻¹) (Fig. 18). Although the frequency distribution of overall retreat rates is scattered and bimodal, it is such that > 80 % of outlet glaciers display retreat rates between 56 and 16 m yr⁻¹. Although slightly different in nature, these overall retreat rates are comparable to calving front retreat rates observed by Carr et al. (2017) for marine-terminating glaciers across Greenland for the 1992–2000 period (mean and 1σ SD ∼ 42 ± 86 m yr⁻¹).

Our reconstructed ice sheet Late Glacial and early- to mid-Holocene retreat rates can be compared with deglacial retreat rates from other ice sheet reconstructions. For instance, Greenland Ice Sheet outlets retreated at paces comparable to, or perhaps slightly faster than, terrestrial margins of the British–Irish Ice Sheet, which is associated with retreat rates of < 50 m yr⁻¹ for most regions during deglaciation (Clark et al., 2022). Marine calving margins of the British–Irish Ice Sheet, however, generally retreated faster than our sampled Greenland outlets (> 50 and up to 451 m yr⁻¹). Greenland outlet retreat rates also appear to be slower than fast-retreating terrestrial margins of the Laurentide Ice Sheet, e.g. the southwest margin, which is thought to have reached retreat rates of 380–340 m yr⁻¹ following separation from the Cordilleran Ice Sheet (Norris et al., 2022). However, the Labrador sector of the Laurentide Ice Sheet is believed to have experienced retreat that is slower than and similar to the Greenland Ice Sheet Holocene retreat, with Lowell et al. (2021) reporting a near-constant mean retreat rate of 52 m yr⁻¹ between ∼ 19 and ∼ 10 kyr BP in this region. We note that our reconstructed Greenland outlet retreat rates are net rates averaged over the distances and time separating 500-year resolution isochrones. Our reconstruction may thus lead to smoothing of faster pulses of retreat (> 200 m yr⁻¹), which therefore should not be ruled out.

Our reconstruction further suggests that, overall, outlet glacier retreat rates varied significantly between regions (Fig. 18). We find that outlet glaciers from the CW region were likely the fastest retreating in Greenland between ∼ 14–13 and 7–6.5 kyr BP, with a regional mean retreat rate of ∼ 50 ± 14 (1σ SD) m yr⁻¹. This region includes the Sermeq Kujalleq (overall rate 49 ± 5 m yr⁻¹), Sermeq Avannarleq (72 ± 14 m yr⁻¹) and Store (70 ± 11 m yr⁻¹) glaciers, for instance. Outlet glaciers from the NE and NO regions also retreated relatively quickly during that period, with mean retreat rates of 44 ± 10 and 44 ± 13 m yr⁻¹, respectively. The retreat of outlet glaciers from the CE and SW regions was intermediate in speed, with mean retreat rates of 35 ± 23 and 29 ± 12 m yr⁻¹, respectively. Finally, we estimate that outlet glaciers from the SE and NW regions experienced slower retreat at that time relative to other regions, with mean retreat rates of 25 ± 5 and 24 ± 15 m yr⁻¹, respectively (Fig. 18). The latter contrasts with present-day retreat rates, as marine-terminating outlet glaciers in SE and SW Greenland currently show the third- and fourth-fastest regional mean retreat rates (136 and 117 m yr⁻¹, respectively) after NO and CW Greenland (Carr et al., 2017). Our reconstruction shows that spatial variability in outlet glacier retreat rate over the studied period is dominated by these regional patterns. Indeed, we find no statistically significant correlation between individual glacier retreat rates and their present-day front width or with the approximate width of the overdeepened troughs along which they retreated. Moreover, we argue that confidence levels in the ice-sheet-wide PaleoGrIS 1.0 isochrones are often too low to conduct quantitative analyses at the valley or outlet glacier scale.

Our results also indicate that for most ice sheet outlet glaciers sampled, the rate of retreat is highly variable throughout the reconstructed time period, i.e. between 14–13 and 7–6.5 kyr BP. Indeed, the median retreat rate of all sampled glaciers (n=72) shows intermediate values between ∼ 14 and ∼ 12 kyr BP (20–10 m yr⁻¹), increases between ∼ 11.5 and ∼ 9.5 kyr BP (30–20 m yr⁻¹), and decreases to its lowest values between ∼ 9.5 and ∼ 6.5 kyr BP (15–5 m yr⁻¹) (Fig. 19). We thus find a significant proportion of Greenland outlet glaciers experienced acceleration in their retreat pace between ∼ 11.5 and ∼ 9.5 kyr BP, when our reconstruction suggests numerous glaciers retreated at speeds > 100 m yr⁻¹ (Fig. 19). These faster retreat rates are not observed as much during other time intervals of the Late Glacial and early- to mid-Holocene. The timing of this faster retreat (∼ 11.5–9.5 kyr BP) coincides with the highest rates of atmospheric and oceanic warming reconstructed between the Younger Dryas (∼ 13–11.7 kyr BP) and the end of the pre-industrial era (Buizert et al., 2018b). Furthermore, using linear regression analyses, we find that the majority of outlet glaciers (65 %) indicate a decelerating retreat pace over the full reconstructed period (Fig. 19).

Figure 18 Map of Greenland Ice Sheet outlet glacier transects (black lines; n=72) drawn to conduct a retreat rate analysis (see Sect. 3.4) using PaleoGrIS 1.0 isochrones, the main results of which are also shown here. Retreat rates over the full reconstruction period (∼ 14–6.5 kyr BP) are denoted using pink circles characterised by sizes proportional to the overall retreat rate magnitudes. Our outermost isochrone (14–13 kyr BP), which marks the outer edge of our transect mapping, is shown with its spatially variable colour-coded confidence levels. Towards the bottom-right corner, a histogram of overall retreat rate distribution is shown to have been split into eight bins of retreat rate magnitude. Moreover, bold black numbers indicate the regional mean overall retreat rate (±1σ SD) for each of the seven ice sheet regions sampled. Towards the top-right corner, a doughnut diagram highlights the relative proportions of outlet glacier retreat rates indicating decelerating (65 %) vs. accelerating (35 %) trends over the full reconstruction period.

Figure 19 Time series denoting the evolution of the median retreat rate for all outlet glaciers sampled over the full reconstruction period (14–6.5 kyr BP). The blue polygon highlights a period of general acceleration in the Greenland Ice Sheet outlet glacier retreat speed (see Sect. 3.4). It is important to note that the youngest time slices are here associated with fewer than 72 data points, as our reconstruction features certain outlet glacier margins reaching present-day extent before others.

4 Discussion

4.1 Progress towards a robust reconstruction of the Greenland Ice Sheet retreat, knowledge gaps, and opportunities for future investigations

PaleoGrIS 1.0 is the first attempt to merge both geochronological and geomorphological markers of former grounded ice extent to produce an ice-sheet-wide isochrone reconstruction of the Late Glacial and early- to mid-Holocene evolution of the Greenland Ice Sheet margin (Fig. 15). It thus represents a significant improvement from previous ice-sheet-wide reviews and compilations (e.g. Dyke, 2004; Funder et al., 2011; Lecavalier et al., 2014; Sinclair et al., 2016). Another novelty of this work is that it aims to make the reconstruction and empirical database as accessible and usable as possible, whether it is to inform future empirical and model–data comparison investigations or to provide a format suitable to the production of future versions of the reconstruction (see online database).

This first version of PaleoGrIS, however, remains incomplete and has its limitations. As noted earlier (Sect. 2.1.1), our mapping of ice marginal landforms was conducted rapidly over a few years and deliberately at a reconnaissance level to identify only the main terrestrial landforms recording former ice margin positions. This was necessary to cover the large area and permitted us to erect the ice-sheet-wide framework we have presented. Conducting and compiling more detailed regional studies that map a wider range of ice marginal landforms, including those located offshore, is an important future task. In regions displaying fewer empirical constraints (e.g. Fig. 20), our reconstruction may act as a stimulus for such investigations and could represent a template that will undoubtedly require revision. We believe such revision would greatly benefit from becoming a wider community effort, enabling the gathering of newly acquired data and knowledge on the update scenarios of local ice sheet margin history (e.g. DATED-2, Hughes et al., 2023).

To enhance confidence in the empirical record and improve our ability to quantitatively compare paleo model simulations with observations, the glaciology community requires more field data from and around Greenland. We here used our reconstruction to identify regions that would most benefit from new constraints on the timing of former grounded ice sheet retreat. To do so, we mapped terrestrial and offshore regions displaying no or only low-confidence event ages and thus presenting either low- or very low-confidence isochrones in our reconstruction. The highlighted regions are shown in Fig. 20, and shapefiles are provided in the online database. Care should be taken, as this assessment does not include ongoing research efforts, data yet to be published, or data published after our census date (21 October 2022) (e.g. Weiser et al., 2023). Nonetheless, we find vast terrestrial and offshore areas displaying low- or very low-confidence isochrones. Offshore of present-day coastlines, we note the majority of the Greenland continental shelf remains data-scarce when compiling studies establishing the timing and pattern of grounded ice retreat from the last full glacial extent. For terrestrial regions, we also find numerous areas displaying potential for new field data collection around the full ice sheet perimeter, with a relatively greater density and surface area coverage of less studied sites in the northern half of the island's periphery (> 71° N; Fig. 20). Furthermore, a great source of uncertainty impeding a better understanding of the deglacial dynamics originates from the ice sheet being generally more extensive today than between ∼ 6 and ∼ 2 kyr BP when responding to the Holocene thermal maximum. Obtaining empirical evidence of the ice sheet minimum extent during that interval (e.g. as is attempted by the ongoing GreenDrill project; Briner et al., 2021) would thus greatly improve our capacity to reconstruct the ice sheet deglacial evolution.

Figure 20 Map of offshore and onshore regions that we interpret as data-scarce for the specific purpose of mapping and dating the deglacial evolution of the grounded Greenland Ice Sheet margin. We thus believe these regions would most benefit from new geomorphological and geochronological data. They were identified by locating regions dominated by very low- or low-confidence isochrones in our PaleoGrIS 1.0 reconstruction. Two polygon shapefiles associated with these highlighted regions can be found in the online database.

4.2 Exploring controls on the varied retreat dynamics of the Greenland Ice Sheet

Our overarching interpretation of the reconstruction is that the deglacial retreat pattern and rate of the Greenland Ice Sheet, and the magnitude of the resulting mass loss, are controlled by a series of interacting mechanisms whose relative importance differs, depending on the spatial and temporal scales analysed.

4.2.1 Ice-sheet-wide response

Linear response to climate and ocean forcing

Our reconstruction shows that at a Greenland-wide scale, the ice sheet mostly retreated in a simple radial pattern progressively stepping back from the continental shelf during the last deglaciation (Figs. 15, 16). When analysed at a 1000- to 500-year resolution, the retreat was not found to significantly diverge from a linear retreat pattern other than with local offsets and complexity arising from lagged positioning of marine versus terrestrial margins (e.g. in and adjacent to outlet glacier tongues; Fig. 12) and where the main ice sheet left residual ice caps behind. As we here focus on the ice-sheet-wide retreat signal (heterogeneous regional behaviours discussed further in Sect. 4.2.2), and because event ages are now reasonably well distributed around Greenland's periphery, we believe that this simple retreat signal would still arise even with the improved spatiotemporal distribution of geochronological markers. Along with this surprisingly straightforward retreat pattern, we observe that the ice-sheet-wide reduction in areal extent reveals a simple (negative) linear relationship with atmospheric or sea surface temperatures (Buizert et al., 2018b; Osman et al., 2021) during the Late Glacial and early- to mid-Holocene (Fig. 21). These correlations suggest climate and ocean forcing were the dominating agents of former ice sheet extent change. For instance, the overall acceleration of outlet glacier retreat that we observe between ∼ 11.5 and ∼ 9.5 kyr BP is coeval with the potent atmospheric warming that characterised the Younger Dryas–early-Holocene transition in the North Atlantic and Greenland regions (Grootes et al., 1993; Fig. 21). This faster retreat also coincides with significant ocean warming, with marine records indicating increasing sea surface temperatures offshore of most Greenland regions from ∼ 16 kyr BP and reaching maximum warmth at ∼ 10–8 kyr BP (e.g. Williams, 1993; Jennings et al., 2006, 2017; Osman et al., 2021; Fig. 21). More particularly, offshore of the SE, SW, and CW Greenland regions, warmer and saline sub-surface Atlantic water ingress from the Irminger Current (Fig. 1) is thought to have caused rapid mass loss and initial retreat of grounded ice margins from the outer continental shelves (Knutz et al., 2011; Ó Cofaigh et al., 2013; Jennings et al., 2017). Similar forcing is thought to have occurred in NE and CE Greenland, with warmer and more saline Atlantic water advected southward from Fram Strait, along the Greenland coast, and across continental shelves to marine-terminating margins, with the establishment of the west Spitsbergen, return Atlantic, and east Greenland currents during deglaciation (Hopkins, 1991; Hebbeln et al., 1994). Furthermore, the findings of former ice sheet margins likely located behind present-day ones during and following the Holocene thermal maximum (Briner et al., 2014), followed by ice margin readvances during the mid- to late-Holocene cooling that peaked with the Little Ice Age (∼ 1850; Kjeldsen et al., 2015), are more evidence of this broadly linear response to climate and ocean forcing. Overall, these observations suggest that the Greenland Ice Sheet was highly sensitive to atmospheric and oceanic warming.

Relative sea level change is also likely to have influenced the timing and the linear nature of the ice-sheet-wide retreat pattern during the Late Glacial and early- to mid-Holocene periods (Figs. 15, 16). In most coastal regions, relative sea level rose rapidly after ∼ 16 kyr BP, reaching a spatially variable high stand of 50–120 m at around 12–10 kyr BP (Funder and Hansen 1996; Gowan, 2023; Lecavalier et al., 2014; Simpson et al., 2009). The timing of this relative sea level high stand therefore also coincides with the observed peak in outlet glacier retreat rates (Fig. 19). More specifically, this rapid sea level rise is thought to have promoted the early deglaciation of interstream sectors of the marine-terminating former ice sheet, which likely presented thinner ice (< 300–400 m) grounded on shallower sections of the continental shelves and were more vulnerable to buoyant lift-off- and calving-induced retreat (Roberts et al., 2009). The importance of this mechanism has moreover been confirmed by modelling studies experimenting with different sea level forcings to nudge modelled Greenland Ice Sheet extent (e.g. Lecavalier et al., 2014; Simpson et al., 2009).

Based on the observation of such broadly linear response in ice extent, we suggest that at a whole ice sheet scale, complexities in ice dynamics are less relevant at 1000- to 500-year timescales than is perhaps perceived when analysing contemporary glaciological fluctuations and feedbacks. Alternatively, our reconstructions are not yet sufficient in density and resolution to adequately constrain retreat and advance oscillations that might have existed.

Glaciological inertia and delay in response

Although the Greenland Ice Sheet was highly sensitive to atmospheric and oceanic warming, specific glaciological responses to a forcing may vary over different timescales. Melting or calving can happen nearly instantaneously, whereas changes to ice accumulation may take thousands of years to work through the system. Dynamic changes to flow geometry and the positions of ice divides, as well as responses to glacial isostatic adjustment, may take even longer still (Rogozhina et al., 2011). Consequently, significant inertia in the system may have caused delay in the ice sheet extent response to deglacial warming. The PaleoGrIS 1.0 reconstruction suggests that the Greenland Ice Sheet margins remained more extensive than present in several regions until 7–6.5 kyr BP. The minimum Holocene extent of the entire ice sheet was thus not reached until later, most likely between 6.5 and 4 kyr BP, as was previously suggested by numerical modelling experiments (e.g. Lecavalier et al., 2014). Over Greenland, δ¹⁸O ice core records (e.g. Alley, 2000; Buizert et al., 2018b) and temperature reconstructions from paleo climate data assimilations (e.g. Buizert et al., 2018a; Erb et al., 2022) all suggest Holocene mean annual temperatures reached maximum values between ∼ 9 and ∼ 6 kyr BP, followed by gradual neoglacial cooling until the end of the pre-industrial era (∼ 1850; Kjær et al., 2022) (Fig. 21). Therefore, our ice-sheet-wide empirical reconstruction agrees with previous observations that, during the Holocene, the ice extent response of the Greenland Ice Sheet lagged the cessation of warming, potentially by thousands of years. This delay in turn caused a lag in the response of the ice sheet to positive mass balances during mid- to late-Holocene cooling. Indeed, while both summer insolation (Berger and Loutre, 1991) and temperature proxy records (e.g. Erb et al., 2022) reveal atmospheric cooling across Greenland from 8–6 kyr BP to the pre-industrial era, significant neoglacial glacier expansion did not occur until 2.5–1.7 kyr BP (Kjær et al., 2022). Finally, the Greenland Ice Sheet's response to anthropogenic warming that started around 1880–1900 was also delayed by almost a century, with the mass and ice extent loss only being of ice-sheet-scale significance from the 1980s onward (Kjeldsen et al., 2015; Yang et al., 2022). Interestingly, the latter decadal delay was shorter than the centennial- to millennial-scale lag in ice sheet's response to the earlier positive mass balances associated with mid- to late-Holocene cooling, thus possibly highlighting a positive correlation between the duration/magnitude of the climate forcing and the length of the delay in the subsequent ice extent response. However, such a lagged response observed at the ice sheet scale was not ubiquitous across the entire ice sheet margin, as certain more dynamic outlet glaciers appear to be better coupled with short-lived early-Holocene temperature changes (e.g. Jakobshavn Isbræ; Young et al., 2013b). Centennial to millennial-scale inertia of the Greenland Ice Sheet following warming has been previously suggested (e.g. Yang et al., 2022) and has important implications for ongoing climate warming because it results in committed future mass losses and thus sea level rise contributions that may last for centuries to millennia (e.g. Greve and Chambers, 2022).

Figure 21 Comparison between our reconstruction of Greenland Ice Sheet areal extent evolution and past atmospheric/sea surface temperature reconstructions from proxy and model data between 20 kyr BP and the present. (a) Blue curve shows the reconstruction of mean annual near-surface temperature at the location of the Greenland Ice Sheet Project (GISP2) ice core (72.58° N, −38.48° W) by Buizert et al. (2018a), who merged Greenland ice core temperature reconstructions (GISP2, NGRIP, and NEEM) with the Transient Climate of the Last 21 000 Years (TraCE-21ka) simulations. Red curve shows the ensemble mean reconstruction of mean annual temperature anomaly (towards GISP2 ice core location) relative to a 5–3 kyr BP reference period, by Erb et al. (2022), who use the global Holocene paleo temperature proxy database of Kaufman et al. (2020) to conduct data assimilation on HadCM3 and TraCE-21ka simulations. (b) The Greenland Summit δ¹⁸O record (Buizert et al., 2018b) that averages the GISP2 and GRIP δ¹⁸O records. (c) Reconstructions of sea surface temperature at different locations offshore Greenland, using the ensemble mean (n=500) from the Last Glacial Maximum reanalysis (LGMR) global dataset of Osman et al. (2021), who use a paleo climate data assimilation scheme to correct biases in climate simulations of the isotope-enabled Community Earth System Model (iCESM) (versions 1.2 and 1.3). Coordinates for time series are as follows: south (Cape Farewell) is 58.34° N, 43.25° W; CW (Davis Strait) is 64.72° N, 56.53° W; CE (Danmark Strait) is 66.64° N, 27.26° W; NE (Nordic Sea) is 72.25° N, 10.10° W; and NW (Baffin Bay) is 73.17° N, 62.65° W. (d) Greenland Ice Sheet areal extent reconstruction from the literature review (LGM–HS1 period), the PaleoGrIS 1.0 isochrones (this study, ∼ 14–6.5 kyr BP), and BedMachine v4 (Morlighem et al., 2017).

4.2.2 Regional responses

Moving from ice-sheet-wide to a focus on individual regions (the ice catchments of Rignot and Mouginot, 2012), we find heterogeneous retreat patterns with apparent differences in both the magnitude of regional areal extent loss and the speed of ice margin retreat (Figs. 17, 18). Spatially variable climate and ocean forcings (Fig. 21) may have caused some of these differences. For instance, the saline Irminger Current transporting warm Atlantic waters offshore of south Greenland (Fig. 1) likely caused this region to experience early, more intense oceanic and atmospheric warming (Levy et al., 2020; Fig. 21c). The inflow of relatively warmer, more saline sub-surface water originating from the North Atlantic current was likely crucial in causing sub-shelf melt near the grounding lines of marine-terminating ice sheet margins (Knutz et al., 2011; Jennings et al., 2014). It has been argued that significant leads and lags in early deglaciation across Greenland may be partly related to the complex advection patterns of this water mass, enabled by the Irminger and west Greenland currents (Fig. 1) influencing ice retreat in SE, SW, and CW Greenland. Similar mechanisms are regarded to have caused the intermediate return Atlantic current to influence ice margin retreat across the NE and CE Greenland shelves (Hebbeln et al., 1994; Nørgaard-Pedersen et al., 2003; Hansen et al., 2022). However, at regional to valley scales, we believe non-climatic and oceanic factors also help explain the variability. We discuss and hypothesise some of these mechanisms in the following sections, with a particular focus on ice-sheet- and sector-scale topographic heterogeneities and asymmetries. While numerous other non-climatic factors, such as spatial heterogeneities in ice shelf loss and sea ice cover, for instance, may also play roles in controlling such regional heterogeneities, we here only focus on what we believe to be main or first-order drivers.

Influence of heterogeneous continental shelf topographies

When the Greenland Ice Sheet reached its maximum extent (between 24 and 16 kyr BP, depending on regions; Funder et al., 2011), and when local atmospheric and, more importantly, sea surface temperatures were lowest (Grootes et al., 1993; Osman et al., 2021), the ice sheet's margins had advanced across continental shelves (e.g. CW Greenland; Ó Cofaigh et al., 2013). Such advances were facilitated by the soft and unlithified nature of shelf surface sediments promoting shear deformation of subglacial sediments, high basal sliding velocities, and ice streaming in some cases (Evans et al., 2009). Certain regions have wider continental shelves providing grounded ice with considerable distances across which to advance during build-up phases, while other sectors only feature a narrow shelf to traverse before the continental slope, deep ocean, and high calving losses are attained. For example, assuming grounded ice reached the continental shelf break in the NE, offshore of NEGIS, the ice sheet had a distance to advance of ∼ 280–300 km relative to present-day positions. In southernmost Greenland, however, this distance to attain the shelf break was instead only ∼ 80–100 km.

During full glacial conditions, precipitation and ice accumulation were very low in Greenland due to cold atmospheric and oceanic temperatures (Buizert et al., 2018). The GISP2 ice core record (Cuffey and Clow, 1997) suggests accumulations < 10 cm yr⁻¹ between 20 and 15 kyr BP, and the latest simulations of the water isotope-enabled Community Earth System Model (iCESM) suggest values of around 5 cm yr⁻¹ for that location (He et al., 2021). Therefore, in sectors experiencing larger lateral ice sheet expansion across wider continental shelves (e.g. NE and CW regions), ice accumulation was less likely to be able to compensate the required mass transfer to lower elevations during extensive margin advance. This potentially caused greater ice sheet thinning and lower surface slopes than in other sectors experiencing less expansion across the continental shelf (e.g. SW, NW, and CE regions).

Relative to a steeper ice sheet profile, a thinner sector of the ice sheet with a gentler surface slope should experience increased melting for a given rise in the equilibrium line altitude because a greater surface area is now below that altitude (Edwards et al., 2014). Warm periods such as the Bølling–Allerød and early-Holocene intervals may therefore have had varied effects on mass balance arising from regional differences in ice surface slope attained during build-up across continental shelves of different widths. Hypothetically, such mechanisms may have caused Greenland sectors with wider continental shelves and more extensive last full glacial advances (e.g. offshore CW Greenland) to have experienced faster grounding line retreat during the initial phase of deglaciation, under still predominantly marine-terminating conditions (e.g. Fig. 22b, scenario 2). Another sector where such effects may have played a role is in NE Greenland, which experienced rapid and extensive deglaciation of a wide continental shelf, while NW Greenland on the opposite side of the main ice divide instead experienced slower and more limited retreat of ice margins (Fig. 22c). Consequently, during the Late Glacial and early-Holocene periods, a significant ice divide migration may have occurred in NE Greenland. The recent findings of a possible reconfiguration of NEGIS during the early Holocene (Franke et al., 2022) could perhaps be linked to such regional adjustments in ice flow. Around Greenland, more generally, the ∼ 12 to ∼ 8 kyr BP window may have been critical for ice stream readjustments, as the ice sheet margin stepped back and thinned from mid-/inner-shelf positions into more topographically confined coastal fjords. Such confinement increased the fixing in position and channelising of ice flow, and pre-existing overdeepenings may have induced some stepped behaviour in retreat. The early-Holocene interval was thus likely associated with a stronger topographic control on ice stream position and behaviour throughout the Greenland Ice Sheet (Roberts et al., 2009).

Influence of heterogeneous onshore bed topography

The present-day onshore Greenland topography is spatially heterogenous, with the CE, SE, and SW regions characterised by prominent coastal mountain ranges, while other regions generally feature lower-lying beds, which are in some places resting below the modern sea level, such as underneath NEGIS or the Jakobshavn and Humboldt ice streams, for instance (Morlighem et al., 2017). The latter can in some cases enhance margin retreat speed and magnitude and promote outlet glacier instability (Jamieson et al., 2012). Higher and more rugged topographies near the coast present orographic obstacles to precipitation from offshore. They promote ice accumulations sustaining coastal glaciers which can contribute ice flux to the main ice sheet and locally decrease the impact of negative mass balances on ice margin retreat during warming phases. Higher coastal topographies enable steeper surface slopes near ice sheet margins, and smaller accumulation area loss for a given rise in equilibrium line altitude during warming (Fig. 22b, scenario 2). Moreover, high coastal mountains likely promote less slippery bed conditions, with a rougher bed and more bedrock areas at higher elevations potentially enabling more cold-based conditions. These mechanisms could together explain why the SE sector of the Greenland Ice Sheet is characterised by a margin located closer to the shore than other more deglaciated regions, despite significantly warmer local atmospheric and sea surface temperatures (Fig. 21; Osman et al., 2021). Indeed, while greater ocean and atmospheric warming likely caused this region to lose its former marine-terminating margins earlier and more rapidly than in other sectors, this high coastal topography feedback may have caused the ice sheet to stabilise when becoming predominantly land-terminating, with local margins remaining near the shore throughout the Holocene.

Generally, we observe a pattern of stabilisation of ice sheet margins during their marine to terrestrial transition in several regions. Indeed, an abundance of prominent moraines and other ice marginal landforms close to present-day coasts can often be observed where present-day topography reaches above the marine limit. This is common in the Ilulissat, Sisimiut, and Nuuk regions (e.g. the Kapisillit stade moraines; Young et al., 2021) (Fig. 12), for example, but also in NE Greenland, onshore of major fjords (e.g. Danmark Fjord). These landforms imply stillstands of the margin once floatation and calving fronts were lost, suggesting that calving and sub-shelf melt caused more ablation than subaerial processes in these regions during the Late Glacial to mid-Holocene retreat of the ice sheet. Here again, the role of bed topography is critical, as the width, depth, and length of fjord systems control how much margin retreat is required prior to ice becoming fully grounded. Due to the added ablation effect of calving, topographic heterogeneities between fjords and inter-fjord sectors can generate different modes of deglaciation for adjacent outlet glaciers, as observed, for example, in north Greenland (Larsen et al., 2020). Therefore, the transition from marine to land-terminating margins, which occurred at different times in different regions due to variabilities in uplift from glacial isostatic adjustment and thus in relative sea level change (Long et al., 2011) is an essential factor influencing the dynamics of former ice sheet margin retreat in Greenland. We note the spatial and temporal uncertainties in the PaleoGrIS 1.0 reconstruction are likely too large to assess such valley-scale feedback quantitatively. Moreover, while numerous regions feature robust paleo sea level records (e.g. the Disko Bugt area; Long et al., 2006), it remains a challenge to establish an accurate, time-dependent, and Greenland-wide map of paleo sea levels at the valley scale (Gowan, 2023). In most Greenland regions, our ability to test the impact of this mechanism requires thorough modelling experiments (with a coupled glacial–isostatic-adjustment model) that are calibrated against observations and run at high (< 1 km) spatial resolutions. Finally, we argue that when subglacial topography is asymmetrical on opposite sides of an ice sheet, it can cause potent ice divide migrations during ice sheet build-up and demise phases (Larsen et al., 2016). In central Greenland, for instance, significantly higher coastal topographies (up to 2500 m a.s.l.) lie underneath and around the eastern ice sheet margin relative to the western margin, which features a wide, low-lying subglacial drainage basin (Fig. 22). During strong cooling phases and ice build-up towards full glacial extent, the ice sheet in CE Greenland draining towards Scoresby Sund would have required thicker ice to overflow mountains and become unconstrained compared to CW Greenland, where ice feeding the Jakobshavn Ice Stream faces no such obstacle. This asymmetrical thickness build-up may have generated an ice sheet dome increasingly skewed towards the east during build-up to full glacial extent, resulting in a potent eastwards migration of the main ice divide in central Greenland. Conversely, during deglaciation, when eastern margins retreated towards or behind present-day positions and became increasingly constrained by coastal mountains impeding ice evacuation, a westward migration of the divide may have occurred in this region (as drawn in Fig. 22). In fact, one could argue that westward migration of the main ice divide in central Greenland might continue and possibly accelerate in the future as the ice sheet margin thins, retreats, and becomes increasingly blocked by higher topographies and reversed bed slopes to the east. We therefore speculate that ice divide migrations likely play an important role in controlling the position of ice streams and influencing ice margin retreat dynamics as noted for other paleo ice sheets (e.g. Greenwood and Clark, 2009). The influence of such mechanisms could be further tested by running high-order and high-resolution model simulations constrained by ice extent observations.

To summarise, we here raise the hypothesis that at the regional scale, the heterogeneous and sometimes asymmetric nature of Greenland's terrestrial and continental shelf topographies may be responsible for the following:

• skewed ice sheet profiles in surface elevation;

• regional variability in thickness and surface slope;

• regional variability in mass balance response to warming as a consequence of the above;

• rearrangements of certain ice streams and of surface velocity fields during retreat;

• spatial variability in the timing of marine to terrestrial transitions resulting in non-synchronous margin stabilisations; and

• major ice divide migrations during both advance and retreat phases.

Together, these rearrangements may cause inter-regional variability in the speed and magnitude of Greenland Ice Sheet margin retreat during deglaciation.

Figure 22 Schematic model describing the hypotheses raised in Sect. 4.2.2, which explore the potential roles of asymmetrical continental shelf and coastal bed topographies on the observed inter-regional variability in ice sheet margin retreat dynamics in Greenland during deglaciation. Perhaps these ideas help explain the regional variability in both rate and magnitude of retreat. The subglacial topography and bathymetry data displayed in the right-hand panels are a merged product of the BedMachine v4 and GEBCO (2022 release) DEMs.

5 Summary and conclusions

• We present PaleoGrIS 1.0, an ice-sheet-wide reconstruction of the Late Glacial and Holocene Greenland Ice Sheet margin evolution, documenting the deglacial back-stepping of ice margins to the present-day ice sheet edge. This new reconstruction is based on around 194 000 geomorphological and 1450 geochronological markers of former ice extent. It yields ice marginal isochrones at a temporal resolution of 500 years and spans ∼ 7.5 kyr from approximately 14 to 6.5 kyr BP.

• This reconstruction of the ice sheet margin evolution through time, based on our Greenland-wide mapping of ice marginal landforms, combined with a synthesis of published geochronological data, enables us to pinpoint the timing and dynamics of key glaciological events, such as the unzipping of the Innuitian and Greenland ice sheets and subsequent opening of Nares Strait, for example, which is here estimated to have occurred at ∼ 8.5–8 yr BP.

• The PaleoGrIS 1.0 isochrones enable us to quantify the Late Glacial and Holocene areal extent change in the entire Greenland Ice Sheet. We find that between 14–13 kyr BP until the present day, the ice sheet has lost approximately one-third of its areal extent (0.89 million km²). Between 14–13 and 9–8.5 kyr BP, the ice sheet experienced a near-constant rate of areal extent loss of 170 ± 27 km² yr⁻¹. We also find Greenland outlet glaciers experienced faster retreat rates between ∼ 12 and ∼ 9 kyr BP, while both rates of areal extent loss and margin retreat decreased markedly after ∼ 9 kyr BP. The timing of such faster retreat coincides with the most potent local climate warming to have occurred between the Younger Dryas (∼ 13–11.7 kyr BP) and the end of the pre-industrial era (∼ 1850). In spite of local glaciological complexities, the deglacial evolution of ice margins at the ice sheet scale is strongly correlated to, and controlled by, changes in atmospheric and oceanic temperatures. These results imply that the Greenland Ice Sheet was highly sensitive to climate warming in its recent history and will likely remain so in the future.

• We find the Greenland Ice Sheet minimum Holocene extent was reached later than ∼ 6.5 kyr BP, which is hundreds to thousands of years after the peak warming of the Holocene thermal maximum. We conclude that inertia in the ice sheet system (a mix of fast and slow glaciological responses) has caused a centennial- to millennial-scale time lag in its response to the cessation of Holocene warming. Understanding and quantifying such inertia, by better calibrating numerical ice sheet models against data, for instance, is of crucial importance to anticipating the magnitude of committed mass loss and the response of the ice sheet to ongoing global warming over the next few centuries.

• When quantifying retreat in different Greenland regions, we find heterogeneous responses in the magnitude and rate of the ice sheet areal extent loss and in outlet glacier retreat rates. While some of this heterogeneity may be related to spatial variability in climate and ocean forcing, other non-climatic factors likely play important roles in determining the regional- to valley-scale dynamics of deglacial ice margin retreat. In particular, we hypothesise that in certain Greenland regions, the asymmetrical configurations of continental shelf and onshore bed topographies on opposite sides of the ice sheet may be responsible for inter-regional variability in the magnitude and timing of margin evolution, with implications for both advance and retreat phases. Ice sheet modelling experiments could usefully explore these hypothesised controls, which may continue to be relevant regarding future ice loss and sea level rise contributions.

• Despite remarkable and accelerating efforts from the empirical community to produce an extensive library of ice extent markers around the periphery of the Greenland Ice Sheet, we find vast onshore and offshore regions with sparse data and which yield isochrones of low confidence. There is therefore much potential for improving our knowledge of deglacial Greenland Ice Sheet margin evolution by producing new robust paleo glaciological investigations in less well-constrained regions. We believe the PaleoGrIS 1.0 reconstruction and database can prove useful for identifying important study sites and motivating future onshore and offshore field activities.

• PaleoGrIS 1.0 is made available as an open-access database, which may prove useful to both the empirical and ice sheet modelling communities and for data-modelling comparison exercises. To reduce data-processing tasks and make our database usable by both communities, we provide isochrone-related datasets in both GIS (GeoTIFF/shapefile) and NetCDF formats. The database also aims to be in a format suitable to the production of improved future versions of the reconstruction.