An ice-sheet modelling framework to determine vulnerable regions of the Greenland Ice Sheet in the past

The contribution of the Greenland Ice Sheet (GrIS) to sea-level rise is accelerating and there is an urgent need to characterize which sectors of the ice sheet are the most vulnerable. Estimating the volume of Greenland ice that was lost during past warm periods can support efforts to constrain the ice sheet's response to future warming. Sub-ice sediment and bedrock, retrieved from deep ice core campaigns or targeted drilling efforts, yield critical and direct information about past ice-free conditions. However, it is challenging to scale the few available sub-ice point measurements to the geometry of the entire ice sheet. Here, we provide a framework for assessing sea-level potential, which we define within an ensemble of ice-sheet model simulations as the amount the GrIS has contributed to sea level when a particular location in Greenland is ice free. An assessment of dominant sources of uncertainty in our paleo ice-sheet modelling, including climate forcing, ice-sheet initialization, and solid-Earth properties, reveals spatial patterns in the sensitivity of the ice sheet to these processes and related feedbacks. We find that the sea-level potential of central Greenland is most sensitive to lithospheric feedbacks and ice-sheet initialization, whereas the ice-sheet margins are most sensitive to climate forcing parameters. We map the GrIS response to warming, in order to (1) estimate of the region(s) of GrIS that likely contributed to the first meter(s) of global sea-level change across a range of plausible deglaciation scenarios, (2) guide future sub-glacial access efforts that can provide targeted information about the response of the ice sheet to past warming, and (3) contextualize existing and future datasets within a glaciologically coherent, full-geometry framework to establish the minimum GrIS contribution to past sea level when a particular location is ice-free. Through our ensemble approach, we can assign a plausible range of GrIS contributions to global sea level for deglaciated conditions at any site. Our results identify primarily areas in southwest Greenland, and secondarily north Greenland, as best-suited for subglacial access drilling that seeks to constrain the response of the ice sheet to past and future warming.

1 Introduction

Sea-level rise (SLR) is one of the most profound economic, social and environmental issues facing humanity. Flooding associated with ongoing SLR is projected to cost up to 3 % of global gross domestic product annually by the end of this century if emissions continue unabated (Jevrejeva et al., 2018). In the United States, SLR disproportionately impacts communities of color and those in low-income areas, exacerbating issues of environmental justice (Hardy et al., 2017). Globally, the displacement of hundreds of millions of people will have cascading social, political, and environmental impacts as populations in low-lying areas, especially in the developing world, are forced inland by rising seas (Geisler and Currens, 2017). Accurately predicting the source of future SLR is critical to adaptation, because the spatial pattern of ice loss impacts that of SLR (Larour et al., 2017; Hamlington et al., 2020).

The rate of global SLR has nearly tripled since 1890 and continued to accelerate since satellite-based methods became widely used in the 1970s (Hay et al., 2015; Nerem et al., 2018). The relatively modest rates of SLR in the 19th and most of the 20th centuries were driven primarily by increased oceanic heat uptake of anthropogenic warming (Hay et al., 2015) and retreating mountain glaciers (Oerlemans, 1994). In the last half-century, SLR accelerated as glacier retreat increased globally (Hugonnet et al., 2021). In the last two decades, ice loss in Kalaallit Nunaat (Greenland) has emerged as a dominant driver of SLR (Mouginot et al., 2019; Coulson et al., 2022), very likely as the result of human perturbations to the global climate system (IPCC, 2019). The Greenland Ice Sheet (GrIS) is predicted to remain the single greatest contributor to SLR over the next half-century (Hanna et al., 2024), motivating the questions of when, at what rate, and from where will the future contributions of the GrIS to global SLR come?

Responses of the GrIS to past periods of naturally forced global warmth may give clues about its future behaviour (Briner et al., 2020). The large-scale response of the ice sheet to warm climate during the Pleistocene (the last ∼2.65 million years) has been the subject of ongoing debate, with some lines of evidence suggesting prolonged resilience and others indicating repeated retreat. Analysis of interglacial ice preserved at the base of the NEEM ice core in Northwest Greenland (Fig. 1) was interpreted to suggest that during the Eemian (Marine Isotope Stage (MIS) 5e, 125 thousand years ago (ka)) the ice sheet contributed only modestly to global sea-level despite >8 °C of warming at that site, with most of the ice-loss at the margin (NEEM community members, 2013); this is consistent with coupled climate-ice-sheet modelling under MIS 5e boundary conditions (Helsen et al., 2013). Reyes et al. (2014) argued that during MIS 11, while a portion of the southern GrIS diminished, the ice sheet lost no more than 30 %–40 % of its volume. Argon isotopes in the basal ice of the southern GrIS (Dye-3; Fig. 1) reinforce this interpretation and demonstrate that some basal ice there is as old as 400±170 ka (Yau et al., 2016). The same approach yielded a basal ice age of 970±140 ka for GRIP ice beneath the summit of the ice sheet (Yau et al., 2016). These findings are in accord with an analysis of the basal material from the GISP2 ice core (Fig. 1), which found that site had been continuously covered by ice for the last 2.7 Ma (Bierman et al., 2014). These results, which argue the GrIS did not completely melt during the Pleistocene, are unsurprising in the context of previous modelling results. DeConto et al. (2008) demonstrated the sensitivity of GrIS glaciation to atmospheric CO₂ levels, which have been near or below 280 ppm deglaciation threshold for at least the past 3 million years (Bereiter et al., 2015; Martínez-Botí et al., 2015); more recent studies have also demonstrated that a perturbation of short enough magnitude and/or duration (e.g. <1 ka) causes changes that are broadly reversible, and do not lead to Greenland's large-scale deglaciation (Bochow et al., 2023).

In contrast, an emerging line of direct observations from the analysis of subglacial materials recovered via ice-core and sub-ice drilling has revolutionized the study of past ice-sheet stability. In Antarctica, studies have focused on the reconstruction of ice-sheet thickness from interior sites using a so-called “dipstick” approach (Halberstadt et al., 2023; Jones et al., 2015), by determining the surface exposure history of nunataks that stick out of the ice-sheet today or at some time in the past. In Greenland, long-archived basal rock and sediment from GISP2, at the centre of the ice sheet, and Camp Century, closer to the margins, are yielding new insight into ice-sheet stability through the Pleistocene. Cosmogenic-nuclide dating (Schaefer et al., 2016), multi-proxy analysis (Christ et al., 2021), and optically stimulated luminescence dating (Christ et al., 2023) in these archives provide in situ evidence that the ice sheet was absent from the site within a particular window of time; although these methods can still be ambiguous about the exact timing of past ice-free conditions, as these methods continue to develop and more samples become available, these windows are narrowing. Balter-Kennedy et al. (2021) added yet another dimension to the possibilities afforded by subglacial bedrock samples, showing that down-core cosmogenic nuclide measurements (e.g. into the bedrock) can yield insights to the total amount of erosion at a particular location, and thus a time-integrated history of ice dynamics (e.g. basal thermal state and basal velocity). Such breakthroughs, combined with a range of other data sources including bedrock elevation, ice dynamics, ice-surface conditions, subglacial lithology, and subglacial thermal state, are currently being used to plan the acquisition of subglacial samples that promise new insights into the ice-sheet history (Briner et al., 2022). But what does this knowledge of a particular location being ice-free tell us about the larger GrIS vulnerability, and in turn, global sea level?

Global to regional sea-level reconstructions have placed constraints on past Pleistocene sea-level highstands (e.g. Dutton et al., 2015, 2022), although disentangling the relative contributions of specific ice sheets is difficult when relying on these indicators (Hay et al., 2014; Vyverberg et al., 2018; Barnett et al., 2023), particularly for magnitudes of sea-level rise on the order of a few meters (Khan et al., 2017; Dyer et al., 2021). Paleo ice-sheet modelling can provide a complementary method for predicting where the ice sheet is most vulnerable, but differences between modelling approaches and in particular model forcing makes it difficult to assess with high confidence which parts of the GrIS margin are the most responsive to warming and at which time (Plach et al., 2018). Previous studies have used ice-sheet models to quantify the response of the GrIS to specific periods of past warmth, and come to different conclusions about the resilience (e.g. Helsen et al., 2013) and geometry (e.g. Helsen et al., 2013; Stone and Lunt, 2013; Robinson et al., 2011) of the ice sheet, even for the same Pleistocene interglacial (e.g. MIS 5e, Plach et al., 2018). Differences in the modelled footprint of the ice sheet result primarily from differences in the experimental design, and in particular, the approach taken to climate forcing (Plach et al., 2018). However, the relative role of different processes and forcings, including uncertainty in surface mass balance, solid-Earth feedbacks, and ice-sheet initialization has not been fully assessed (Edwards et al., 2014), and is especially lacking at regional scale.

Here, we outline a novel approach using an ensemble of ice-sheet model simulations, not tied to a particular interglacial but rather encompassing a range of past warming scenarios, to define sea-level potential throughout Greenland. We apply this methodology across Greenland to reveal regional patterns in the sensitivity of the ice sheet to different processes and feedbacks, and demonstrate how this information can be used to guide sub-ice drilling efforts to gain the greatest insights into past sea level. In documenting the utility of this methodology, we also suggest ways forward for incorporating yet additional sources of uncertainty and providing data-based constraints that will further strengthen the inference of sea-level potential in the future.

Our methodology provides a framework for answering questions motivated by the recovery of subglacial materials from specific regions (e.g. Christ et al., 2023). Given that we know a location was ice-free at a certain time, this framework constrains how much mass the ice sheet must have lost for that to happen, providing both a robust lower-limit and a range of plausible ice volume loss numbers. The approach allows us to pose, and answer, critical questions about the future of the ice sheet with a novel perspective, such as “what are the most likely source regions for future contributions to sea-level rise from Greenland?”

2 Methods

Here we describe our modelling approach and the decisions made with respect to each of the ensemble parameters (Table 1).

Table 1 Summary of parameters used to define the ensemble. Detailed information on parameter choices in provided in Sect. 2.2.1 through 2.2.5.

We have framed this study in terms of five model parameters that have been established as having a major impact on deglaciation processes across a range of timescales, and thus represent the major uncertainties in paleo ice-sheet reconstruction: ice-sheet initialization (Solgaard and Langen, 2012; Goelzer et al., 2018), lithospheric relaxation time (e.g. Austermann et al., 2015; Abe-Ouchi et al., 2007), precipitation lapse rate, spatial climatology pattern (Helsen et al., 2013), and rate of interglacial warming (Bochow et al., 2023). We aim to establish end-members for each parameter that illustrate the main sources of uncertainty in defining sea-level potential while leaving room for future efforts to explore the parameter space in finer detail (Table 1). Two major factors drive temperature change over Greenland on multi-centennial (i.e. deglacial) timescales: CO₂ (Buizert et al., 2018) and insolation (Helsen et al., 2013), resulting in different patterns and magnitudes of warmth; yet our knowledge of the distribution of past temperature change is limited, particularly for past warm periods. To encompass a range of spatial temperature patterns, we use two different starting climatologies representative of warmth driven by greenhouse gases (pre-industrial) versus high boreal summer insolation (Holocene Thermal Maximum; HTM), which are well-constrained from a synthesis of ice-core, observational, and model-derived outputs (Buizert et al., 2018; Fig. 1). Temperature change can be rapid or slow; we use four rates of temperature change to account for this. Deglaciation may begin from an ice-sheet configuration that looks similar to modern or it may have continued straight out of a glacial maximum with the ice margin advanced to the continental shelf (Funder et al., 2011); we use three different starting geometries to account for this. We also examine the impact on deglaciation patterns of increasing precipitation as a function of warming and different lithospheric response times, which determine the rate of glacial isostatic adjustment. Through our ensemble approach, we examine the response of the ice sheet to many different scenarios, and integrate the responses to constrain the sea-level potential of any particular site (defined as the amount the GrIS has contributed to sea level when a particular location in Greenland is ice-free). For the error propagation, we encapsulate multiple sources of uncertainty by capturing different end-members in the climate forcing, initialization, and solid-Earth model, thereby placing uncertainties on estimates of sea-level potential that stem from these unknowns. Each unique set of parameters is subject to four different rates of atmospheric warming, allowing us to capture how uncertainty in these parameters affects the way the ice sheet retreats under diverse warming scenarios (Table 1). We map the GrIS response to warming, in order to (1) estimate of the region(s) of GrIS that are most likely to contribute to the first meter of global sea-level change, (2) guide future sub-glacial access efforts that can provide targeted information about the response of the ice sheet to past warming, and (3) contextualize existing and future datasets within a glaciologically coherent, full-geometry framework to establish the minimum GrIS contribution to past sea level when a particular location is ice-free (e.g. Christ et al., 2023). From our resulting map, we can infer the range of plausible sea-level potential of any part of the GrIS, regardless of when the most recent deglaciation occurred. Stated differently, the map illustrates which segments of the GrIS are most vulnerable under diverse climatic warming scenarios. We first provide details on the model and simulation set-up, and then explain our calculation of sea-level potential.

Figure 1 Ice-sheet model forcing and initialization. (A) Two climatologies are used to initialize the climate forcing (Mean Annual Surface Temperature is shown). The first is from the Holocene Thermal Maximum, and the second is modern (preindustrial). The difference between the two climatologies shows that the HTM climate is warmer in North and West Greenland by up to 2 °C. (B) Three starting ice-sheet configurations are used in the ensemble: a “paleo-start” initialized by running the model through a glacial cycle, a modern or “present-day start” initialized based on modern thickness observations (Morlighem et al., 2017), and a modelled Last Glacial Maximum ice sheet. Thickness differences are shown relative to BedMachine (Morlighem et al., 2017). Coloured dots indicate deep ice-core sites discussed in the text: Camp Century (Deep Red), NEEM (Pink), NGRIP (light green), EGRIP (orange), GRIP/GISP2 (dark green), Dye-3 (light blue).

2.1 Ice-Sheet Model

We used a three-dimensional thermomechanical ice-sheet model (Penn State Ice-Sheet Model) that uses a hybrid ice-flow law that efficiently bridges between fast-flowing areas of streaming ice (Shallow Shelf Approximation) and inland areas of low velocity and high driving stress (Shallow Ice Approximation) (Pollard and DeConto, 2012a). We solve for the thermal state of the ice-sheet by considering vertical diffusion of geothermal heat from below (constant geothermal flux of 42 mW m⁻²) and above (surface climate), heat generated by internal friction and friction at the basal boundary, and advection from new accumulation and ice-flow. Isostatic adjustment is calculated using an elastic lithosphere-relaxing asthenosphere model, where the response time of the bedrock to a changing ice and ocean load is a free parameter that we vary in our ensemble. The model has been validated against other ice-sheet models under a range of conditions (e.g. Cornford et al., 2020; Pollard and DeConto, 2020) and has been used extensively for paleo and future ice-sheet simulations in Antarctica (DeConto et al., 2021) and the Northern Hemisphere (Han et al., 2021). The model uses a Weertman-type sliding law for basal ice motion with basal sliding coefficients calculated through an inverse scheme that iteratively adjusts sliding to reduce the mismatch between the modelled and observed ice-sheet geometry for the present-day ice sheet (Pollard and DeConto, 2012b). For our experiments, we calculate basal sliding coefficients through an inversion with modern climate forcing everywhere beneath the ice-sheet; ice-free areas are then assigned a sliding coefficient that reflects one of two end-members based on their elevation, with bedrock below sea-level assigned a weak bed and subaerial bedrock assigned a strong bed. For an initial state at LGM, the ice margin advances onto the submarine continental shelf with a relatively weak bed (Fig. 1b). The model has been applied to understand paleoclimate scenarios where both the boundary conditions and model forcing differ substantially from modern-day. It is therefore well suited for the relatively long integrations run here, and even longer (100 ka–1 Ma) integrations that offer complementary ways to study the evolution of the ice sheet on glacial-interglacial timescales in future work. We use a linear temperature lapse-rate correction of 5 °C km⁻¹ (Abe-Ouchi et al., 2007) to downscale the 40 km² climate forcing to the 10 km² ice-sheet grid, and to dynamically adjust the ice-sheet surface temperature as the ice geometry evolves (e.g. height-mass balance feedback, Weertman, 1961).

2.2 Ensemble Design

We ran ninety-six simulations varying four key parameters: starting climatology, lapse rate for precipitation, aesthenosphere relaxation time, and starting geometry. Each combination of parameters was subject to four different rates of glacial-to-interglacial warming applied for 10 000 years (see Sect. 2.2.3). In analysing the ensemble, we focus on one set of experiments initialized with a modern geometry and one set of experiments initialized with an LGM geometry so as not to weight the results toward modern, resulting in 64 simulations evaluating sea-level potential for each 10 km by 10 km grid box. This approach and its implications are discussed in greater detail in Sect. 2.2.5.

2.2.1 Initial Climate Forcing

A primary control on the spatial pattern of GrIS deglaciation is the pattern of surface mass balance (SMB). Climate reconstructions exist continuously for the last 21 kyr at a range of spatial and temporal resolutions, because during this period ice cores, climate models and modern observational data overlap (Buizert et al., 2018; Badgeley et al., 2020; Osman et al., 2021). These reconstructions can be used to infer past patterns of surface mass balance to force ice-sheet models. Conversely, for most past warm periods before 21 ka, we have less confidence in high-resolution climate reconstructions. Thus, we select two representative time periods from the Holocene to represent end-members in the SMB forcing (Fig. 1). First, we select a time slice in the Holocene Thermal Maximum at 8.5 ka when warm (especially summer) temperatures were driven by Earth's orbital configuration, resulting in a more developed ablation zone in northern Greenland, and a reduced ablation zone in western Greenland (Fig. 1). The second time-slice chosen is pre-industrial (1850 CE), when warm annual temperatures were driven by increased atmospheric CO₂ relative to glacial periods, resulting in minimal melting in northern Greenland and a well-developed ablation zone in western Greenland. These timeslices thus represent two end-members: lower CO₂/high insolation (HTM) and higher CO₂/lower insolation (preindustrial). In this way, they are representative of two known modes of interglacial warmth, and capture both the spatial and seasonal patterns associated with them. However, they do not capture different spatial/seasonal patterns that might be associated with climates warmer than modern, and we assume both spatial and seasonal patterns stay fixed as we conduct the interglacial warming experiments. Both forcings come from a blended model-data reconstruction that includes seasonally resolved spatial and temporal variability (Buizert et al., 2018) downscaled to 40 km resolution (Fig. 1, Sect. 2.2.5).

2.2.2 Lapse rate applied to precipitation

Ice-sheet margin retreat inland of its present-day position during past warm periods can be volumetrically offset by increased accumulation in precipitation-limited areas whose temperature remains far below the melting point. As climate warms, the atmosphere's capacity to hold water is enhanced, which can lead to increasing precipitation rates (and thus ice sheet thickening) for inland ice-sheet regions (Payne et al., 2021). We account for this by considering a precipitation correction that increases precipitation by 2 % per degree of temperature increase in each grid cell. This “precipitation-lapse rate correction” enables us to consider the impact of the feedback between a warming atmosphere and its moisture content/capacity in calculating ice-volume changes. We consider this value to be a plausible upper-bound, as it has been found to accurately reproduce glacial-interglacial changes in precipitation rate (Ritz et al., 2001; Abe-Ouchi et al., 2007). Simulations where the precipitation-lapse rate correction is ignored do not have a mechanism for increased precipitation to offset ice loss as temperature rises and melting increases, so these ensemble members are expected to predict greater total ice-sheet mass loss prior to deglaciation.

2.2.3 Rate of interglacial warming

The GrIS volume decreased in response to past variations in natural forcing, including orbital changes, changes in ocean circulation, and atmospheric greenhouse gasses. Mechanisms driving past climate warming vary among interglacial periods (Past Interglacials Working Group of PAGES, 2016) and the timing of ice-sheet retreat during pre-LGM interglaciations is only coarsely constrained (Schaefer et al., 2016), so we leverage the best-studied periods of past ice-sheet retreat to understand possible rates of interglacial climate warming. During the last deglaciation, Greenland's mean annual temperature increased by ∼18 °C between 18 and 12 ka, an average rate of 3 °C per millennium (Buizert et al., 2014). However, the total temperature change during the summer season, which is largely responsible for controlling ice-sheet melt, was closer to 12 °C (Buizert et al., 2018). During the early Holocene, more muted warming (∼3 °C over 3 kyr) drove GrIS retreat to behind its present-day margin in many sectors (e.g. Bennike and Weidick, 2001; Larsen et al., 2016; Young et al., 2021). Thus, to capture a reasonable range of warming rates based on paleoclimate evidence in our ensemble, we subject the ice sheet to an idealized interglacial warming rate ranging from 1.0 to 2.0 °C kyr⁻¹ in increments of 0.33 °C. We choose to apply this rate of warming for 10 000 years as a representative Pleistocene interglacial length, and note that the sea-level potential for a given site may be different for an interglacial of a different duration. For a shorter-duration interglacial, the subset of ensemble members that deglaciate most quickly could be analysed. For a longer-duration interglacial, the simulations could be run further forward with continued warming.

2.2.4 Solid-Earth relaxation time

Solid-Earth dynamics influence ice-sheet stability (Austermann et al., 2015) and have changed beneath Greenland as a function of time (Rogozhina et al., 2016) and potentially in response to fluctuations of the GrIS itself (Stevens et al., 2016). A great deal of uncertainty remains in the details of the solid-Earth beneath the ice sheet (Kappelsberger et al., 2021). Here, we choose to use two lithospheric relaxation times (500 and 3000 years) to represent end-member scenarios, with the former representing hot, low-viscosity (fast-responding) mantle like that underlying northeast Greenland today (Fahnestock et al., 2001), and the latter is a standard value for relaxation time that has been calibrated against measurements of glacial-isostatic adjustment in Antarctica (Le Meur and Huybrechts, 1996; Coulon et al., 2021). The strength of this approach within our modelling framework is that it allows us to quantify where solid-Earth processes are likely critical for understanding sea-level potential without having full knowledge of the true Earth structure beneath Greenland.

2.2.5 Ice-sheet initialization

Following its expansion to the continental shelf edge at the end of the Last Glacial Maximum (LGM; 21 ka), rising global CO₂ drove the GrIS to recede and approach its present-day margin (Cuzzone et al., 2019). A warm summer orbit during the HTM (∼8 ka) led to continued ice-sheet margin retreat inland from its current position (e.g. Young et al., 2021). These two phases of retreat illuminate two distinct ways that Greenland may have deglaciated more fully in the past: either quickly following a glacial period, when the GrIS was defined by extensive marine-terminating margins, or after reaching a modern-like “interglacial” state, when much of the ice-sheet is land-terminating and ice-ocean interactions are mostly confined to narrow fjords. To capture both these possibilities, we start our simulations from either an ice sheet that has been run to equilibrium with LGM climate conditions, or a modern ice-sheet. For the latter, we ran a set of simulations present-day ice sheet initialized to match ice thickness observations and a set of paleo-start simulations that has been spun-up to modern through a glacial cycle (Buizert et al., 2018). To initialize the present-day ice sheet, we used an observational data set of ice extent and thickness (Morlighem et al., 2017). For the paleo-start simulations, we equilibrated the sheet with an LGM climate forcing (Buizert et al., 2018) for 80 kyr followed by an evolving climatology for 21 kyr (Buizert et al., 2018). The reconstruction combines a reanalysis product that uses meteorological station records and regional climate model output (Box et al., 2009), three ice-core-based temperature reconstructions (Buizert et al., 2014), and a transient simulation of the last deglaciation (Liu et al., 2009) to produce a high-resolution, transient reconstruction of the last 21 kyr; a complete description of the methods used can be found in Buizert et al. (2018). The two modern initial ice sheets differ in their geometry near the margin but have a similar ice-volume (Fig. 1b).

Figure 2 Ensemble design. An example of our results shown for one location in West Greenland. (A) Ice thickness at the green dot in panel (C) plotted for all ensemble members. Each simulation is represented by one thin line. Simulations that reach thickness =0 at some point during the deglaciation are used to calculate sea-level potential for this site. Purple and red lines correspond to purple and red histograms in panel (B). (B) Histogram of outcomes for the location shown with the green dot in panel (C). The contribution of Greenland to global sea level when this site becomes ice-free ranges from 2.0 to 3.2 m. The full ensemble is shown in grey. The ensemble members which all have the precipitation lapse rate turned off are superimposed on the histogram in purple. The ensemble members with a HTM climatology are superimposed in orange. In our ensemble, the sea-level potential for this site is most sensitive to spatial climatology, because knowing that parameter with certainty would reduce the spread of the ensemble by the greatest amount. (C) Greenland footprint associated with ice-free conditions for the location in West Greenland identified with a green dot. Blue is land. Black regions indicate that every simulation is ice-free at the same time that this location deglaciates, whereas white regions are still ice-covered in every simulation when this location becomes ice-free. Other colors indicate the number of simulations predicting deglaciation at that location.

2.3 Sea-level potential calculation

Figure 2 shows the method that we apply to calculate sea-level potential, and the sensitivity of the sea-level potential to each of our ensemble parameters. For each 10 km model grid cell, we analysed the ensemble to find the first timestep the site becomes ice-free in each simulation (Fig. 2a). For the first ice-free timestep in each simulation, we save the ice-sheet volume and extent, and collate these values across the whole ensemble (Fig. 2b, c). We define the median of the ice-sheet volume histogram for each location as the sea-level potential for that location. For each site, we consider the uncertainty to be the ensemble spread. In our method, the sea-level potential is specific to each site and to our ensemble. Repeating this process for each grid cell, we also look at our results along every dimension of our ensemble, enabling us to calculate the importance of each parameter for each site and rank which of the considered parameters dominate the ensemble spread at each location. We define the sensitivity to each parameter as the subensemble spread for each parameter separately divided by the full ensemble spread. When these spreads are equal, the sensitivity score is one and the ensemble spread does not change if we consider only a particular value of that parameter. When a parameter spans a smaller range than the full ensemble does, this number is less than one, indicating that knowledge of this parameter reduces the spread of sea-level potential. The parameter that each location is most sensitive to is determined by identifying the smallest sensitivity score, and the relative importance of each parameter is determined by the ranking the scores from lowest to highest.

Figure 3 Greenland's sea-level potential. (a) Colors indicate sea-level potential, defined as the median amount that Greenland has contributed to global sea level when that grid cell has become ice-free. Size of each dot indicates the uncertainty (quantitatively, the ensemble spread as in Fig. 2b; smaller dots represent a greater uncertainty). Black outline highlights regions where ice-free conditions are associated with median sea-level potential less than 2 m, and the spread is less than 1.5 m. Dots are plotted for areas where the ice thickness is greater than 600 m and sea-level potential is greater than zero within a 20 by 20 km area; areas within the modern limit of the ice sheet (underlying white area) that do not meet these criteria are not plotted. (b) Sea-level potential only (meters sea level equivalent). (c) Confidence: Histogram width only (meters sea level equivalent).

3 Results

Sea-level potential is generally lower near the ice margins (the first regions to deglaciate in our ensemble) and higher in both Central Greenland and in areas of high topography along the southeast coast (the last regions to deglaciate in our ensemble) (Fig. 3b). In southern Greenland, sea-level potential increases gradually from the west coast toward the ice divide. In southeast Greenland, areas of high sea-level potential along the margin are resilient to deglaciation, though there is a region between the southern dome and the main dome of the ice-sheet with low sea-level potential. In the north, sea-level potential is generally low near the ice margin and increases towards the ice divide, with a more gradual increase in the northeast and a sharper increase in the northwest. An area of higher sea-level potential extends from the northwest corner of the GrIS inward toward central east Greenland, forming a core of high sea-level potential there.

The uncertainty in sea-level potential is greatest in North and South Greenland, indicating that ice-free conditions here are associated with a range of ice-sheet geometries (Fig. 3c). The lowest uncertainties are well-correlated with the regions of highest sea-level potential, near where the ice sheet is thickest today and along the southeast coast, where ice caps covering the alpine peaks are the last vestiges of the GrIS to melt (Fig. 3c). However, there are also relatively low uncertainties throughout West Greenland, which generally decrease towards the ice divide, indicating that the initial stages of deglaciation show a greater variety of ice geometries but as deglaciation progresses these geometries tend to converge. We combine the sea-level potential (median of histogram) and uncertainty (width of histogram) to produce a map that highlights areas that (a) deglaciate when the ice sheet has contributed less than 2 m to SLR and (b) have an ensemble spread (histogram width) of less than 1.5 m (Fig. 3a). These regions indicate likely regions for the first 1–2 m of ice loss in the past. This map reveals regions in Greenland where ice-free conditions are associated with a narrow band of contributions to sea level (Fig. 3a), primarily in west Greenland, but also at select coastal sites in northwest and northeast Greenland (Fig. 3a).

Figure 4 Parameter sensitivity test. (A) Shows which ensemble parameter exerts the strongest control on the ensemble spread of sea-level potential. (B) Sensitivity to starting the simulation from Last Glacial Maximum conditions. (C) Sensitivity to a reduced lithospheric relaxation time. (D) Sensitivity to neglecting a precipitation lapse rate correction. (E) Sensitivity to starting from a climatology from the Holocene Thermal Maximum.

We find that initial geometry, aesthenosphere relaxation time, lapse rate for precipitation, and starting climatology play a dominant role in determining sea-level potential for different parts of the ice sheet in our ensemble. The starting climatology and precipitation-lapse rate generally play a greater role near the ice-sheet margin, with lithospheric response time and starting ice-sheet geometry playing the dominant role in inland regions (Fig. 4a). To identify the most important drivers for each region, we generated an estimate of the parameter sensitivity for each of our four ensemble parameters (Fig. 4b–e). We find that North and West Greenland are most sensitive to the choice of starting climatology, driven by the differences in SMB patterns between the chosen reconstructions for those sectors (Fig. 4a). However, broad regions in Northwest and South Greenland are most sensitive to the precipitation-lapse rate scaling. In Central and East Greenland, both the initial configuration (LGM versus modern) and a more responsive solid-Earth are the main factors that drive variance in the ensemble. Across our ensemble, these regions are consistently the last to deglaciate. Thus, we find that areas near the margin are most sensitive to differences in surface mass balance, whereas inland regions are more sensitive to slower processes like glacial isostatic adjustment that impact ice-sheet evolution over thousands of years. In general, regions that are sensitive to the starting geometry are also sensitive to lithospheric relaxation, and regions sensitive to climatology are also sensitive to precipitation lapse rate. However, there are some exceptions to this trend; for example, in North Greenland climatology and starting geometry are the two parameters that have the greatest impact on sea-level potential.

Figure 5 Ensemble illustration. For any location beneath the ice sheet (here a location in Northeast Greenland, marked with a red box), our approach allows us to determine which ensemble parameters are responsible for the greatest amount of uncertainty in sea-level potential. Here, for each of the four ensemble parameters tested, the minimum and maximum configurations of the ice sheet when that location is ice-free are illustrated. We consider sea-level potential at this location relatively insensitive to the fast-responding lithosphere and HTM climate, because ice-free conditions at this site are associated with a wide range of ice-sheet geometries for a given value of this parameter. In contrast, for a given value of precipitation lapse rate or initial ice-sheet configuration, the ice-free geometries are more similar. Thus, improved knowledge of this parameter when this location is ice-free results in a more tightly-defined sea-level potential.

4 Discussion

4.1 Parameters underlying variability in sea-level potential

Incorporating multiple sources of uncertainty in our ensemble design allows us to look in detail at how each parameter impacts the sea-level potential at each site (Fig. 5). We note that sea-level potential, and its associated uncertainty (our ensemble spread) is ensemble- and model-dependent; here we have endeavoured to design an ensemble that captures major modes of interglacial climate variability and used an ice-sheet model that is well-suited to studying periods with different boundary conditions, and future work to explore a larger ensembles and compare multiple models will be beneficial. Nevertheless, in this work, we find that areas with a wide spread of sea-level potential (e.g. Camp Century; Christ et al., 2023) can be thought of as places where ice-free conditions at that site are associated with a wide range of potential ice-sheet geometries; the ice sheet can grow and shrink in different regions while ice remains at that site. For example, we find sea-level potential for the Camp Century site is 3.2 m with a spread of 1.4–5.6 m SLE. An analysis of subglacial material from Camp Century found that site was ice-free during MIS 11; our method thus provides a plausible lower bound on the GrIS contribution to global sea level at that time (Christ et al., 2023). Uncertainty associated with the sea-level potential could be reduced by adding constraints on simultaneously ice-free conditions at more than one location or by considering a subset of our parameter space. For example, at Camp Century, simulations that start with an LGM ice sheet require that the ice sheet loses +2.2 m of sea level equivalent prior to Camp Century deglaciating, and simulations that use a HTM-like climatology require +2.7 m before Camp Century deglaciates (Christ et al., 2023). Thus, knowledge of when a site deglaciated (i.e. during a period of insolation-driven enhanced summer warming relative to present) can inform our interpretation. Our work shows that parameter sensitivity varies geographically; for example, locations nearby to Camp Century are more sensitive to initialization and less sensitive to spatial climatology (Fig. 5). For many sites, the shape of the full histogram suggests that future work to more fully sample the parameter space will be beneficial; nevertheless our end-member approach here shows that the uncertainty in sea-level potential can be greatly reduced in some regions through knowing one of our ensemble parameters more precisely.

Southwest and North Greenland, where there are broad areas of low sea-level potential that could be accessed by subglacial sampling, are most sensitive to the spatial climatology pattern (Fig. 4a). The inclusion of a precipitation-lapse rate correction and initializing the simulations with a LGM ice sheet geometry are the dominant parameters in some sub-regions, for instance in Northwest Greenland and Southwest Greenland. Considering the sensitivities of each individual parameter, the spatial climatology, precipitation-lapse rate, and LGM initialization all play some role in controlling the ensemble spread in the regions of Greenland where the first few meters of SLR are likely to be sourced. In contrast, accounting for a shorter lithospheric relaxation time only impacts the ensemble spread around the most resilient portions of the ice sheet; by the time the ice margin has reached these areas, Greenland has most likely contributed >4 m to SLR (Fig. 3b). This may reflect a critical role for solid-Earth processes in dictating the location of the ice-sheet margin in Central Greenland and aligns with a region that has previously been argued to have a higher geothermal flux and a more viscous mantle (Fahnestock et al., 2001; Rogozhina et al., 2016; Stevens et al., 2016). While lithospheric response exerts a dominant control on sea-level potential in Central Greenland, this source of uncertainty is not likely to impact the regions where the first meter of SLR will come from.

Initializing the model with a LGM ice-sheet geometry has the greatest impact in Central Greenland, where the LGM ice sheet was thinner than the modern ice sheet due to an arid LGM climate. However, this parameter is of secondary importance for North and South Greenland and has the least impact in West Greenland. Neglecting a precipitation-lapse rate has the strongest control on the ensemble in Northwest and South Greenland, where separate ice domes exert a strong control on ice dynamics (Fig. 4d). The dominance of the precipitation-lapse rate illustrates the importance of accounting for changes in precipitation as temperature and surface elevation changes over peripheral ice-domes during periods of deglaciation (such as Northwest and South Greenland), as this can increase resilience to the elevation-surface mass balance feedback (Weertman, 1961; Edwards et al., 2014).

The use of a HTM climatology influences deglaciation in North, West, and South Greenland. Holocene melt records are available in North Greenland, including at NEEM (NEEM community members, 2013), GISP2 (Alley and Anandakrishnan, 1995), and Agassiz ice cap (Koerner and Fisher, 1990). The climate record from NEEM also indicates a greater sensitivity to HTM conditions, showing an early Holocene mean annual warming of 6 °C, relative to 2 °C at Summit (Lecavalier et al., 2017; Dahl-Jensen et al., 1998). Our ensemble does not include variations in ocean or indirect sea ice forcing, which likely played a role in past deglaciation scenarios (Irvalı et al., 2020). In particular, although the precipitation lapse rate we apply reveals regions where changes in precipitation strongly impact deglaciation, changes in precipitation that are not associated with temperature changes (e.g. changes in moisture pathways associated with changing sea-ice cover; Koenig et al., 2014) are not captured by our ensemble set-up. At the fjord scale, ocean warming can have a distinctive impact on ice-sheet dynamics and thus should be considered in future work (e.g. Straneo and Heimbach, 2013; Wood et al., 2021). However, because the modern GrIS is mostly terrestrial and the influence of ocean forcing is often limited to within ∼10 ice thicknesses of even marine-terminating glaciers (Felikson et al., 2017), ocean forcing is expected to have less of an impact on deglaciation than changes in surface climate – although future work to directly test this hypothesis is warranted, in particular for looking at time periods when the GrIS was larger than today but smaller than glacial maximum conditions. In addition to playing an important role in the modern-day ablation zone of West Greenland, central northwest Greenland is particularly sensitive to this parameter. This area corresponds to a low-lying part of Greenland's topography and is on the ice divide connecting northwest Greenland with the central dome of the ice sheet. The dominance of the climatology here reflects the important role of HTM-like conditions (enhanced warming in the North and West) for driving deglaciation further once the northwest GrIS has disintegrated.

A key control on patterns and rates of deglaciation in regions of low sea-level potential is the applied SMB forcing (e.g. Plach et al., 2018). In our ensemble, the starting SMB fields, particularly the extent of the ablation zone, play an important role in ice-sheet geometry during deglaciation across all scenarios and play a dominant role in sea-level potential for regions that consistently deglaciate before Greenland has contributed its first meter to SLR. Surface mass balance is difficult to accurately reconstruct for past interglacial periods (e.g. Helsen et al., 2013). Our approach avoids direct reconstruction of SMB for a particular interglacial by considering a range of plausible forcings and identifying the range of sea-level potential associated with the uncertainty in the climate forcing. By including both a pre-industrial and HTM climate forcing, we capture two well-documented modes of interglacial climate in Greenland (e.g. Buizert et al., 2018). However, other modes of surface climate are possible, and may become dominant in the future as boundary conditions and forcings evolve (e.g. Koenig et al., 2014; Sellevold et al., 2021). Our climate forcings come from one reconstruction, and future work to include other work that reconstructs climate using different methods would complement and expand our analysis (e.g. Badgeley et al., 2020). Nevertheless, our results confirm the primacy of correctly predicting the spatial patterns of climate over Greenland (Edwards et al., 2014) for inferring past sea-level change, and suggest that selecting sites that have lower uncertainty in their sea-level potential will increase the impact of subglacial observations.

4.2 Comparison with other modelling studies

Ice-sheet modelling experiments investigating GrIS response to past warmth have resulted in divergent conclusions about ice-sheet stability. Many previous studies found that West Greenland responded most strongly to past interglacial warm periods (e.g. Greve, 2005; Robinson et al., 2011; Born and Nisancioglu, 2012; Helsen et al., 2013; Sommers et al., 2021). At the same time, other studies have found that North Greenland is also highly sensitive to past interglacial warmth (e.g. Stone and Lunt, 2013). Some studies show both West and North Greenland responding to past warmth simultaneously (Robinson et al., 2011; Born and Nisancioglu, 2012; Aschwanden et al., 2019). Our sensitivity-mapping approach allows us to consider how and why these results may differ from other studies that modelled Greenland deglaciation patterns. For example, we find that whether Northern Greenland is an early contributor to SLR is dependent on the choice of a HTM-like climate forcing (Fig. 4a). Our approach is distinct because rather than considering one particular warm period, our ensemble encapsulates a range of deglaciation scenarios and treats them as all equally likely. This allows us to overcome the challenges associated with perfectly simulating a particular time period in favor of identifying the patterns that are common among different styles of deglaciation.

4.3 Implications for interpreting sea-level records

Our approach complements far-field sea-level records by providing a method to quantify a minimum sea-level contribution for the GrIS during periods when directly dated subglacial samples are available, regardless of the final geometry of the ice sheet at the time of maximum retreat during an interglacial (e.g. Dyer et al., 2021; Barnett et al., 2023). Terrestrial records from West Greenland have revealed this area was particularly sensitive to warming during the HTM (e.g. Larsen et al., 2016; Young et al., 2021) and our results confirm this as a persistent feature of GrIS response to warming. Although most far-field sea level records of past warm periods have focused on the maximum sea level attained, our work may provide a framework for thinking about where we may look for evidence of retreat early in a warm period. Our ensemble identifies SW Greenland as earliest to deglaciate regardless of uncertainty in the climate forcing and other parameters examined here. In part of the parameter space, namely those simulations with a HTM climate forcing and no precipitation lapse rate, north Greenland is also among the first regions to deglaciate.

At present, the ensemble has been designed specifically to demonstrate how subglacial material documenting past ice-free conditions, in combination with numerical ice-sheet modelling, can provide estimates of continent-wide sea level contribution that account for dominant sources of uncertainty. Due to hysteresis effects (Robinson et al., 2012), sea-level potential may be different for a growing versus shrinking ice sheet; here we focus on the application of sea-level potential to presently ice-covered sites as they become deglaciated, but we highlight that future work to examine the hysteresis of sea-level potential is warranted and may be useful for understanding what feedbacks and processes are most important for ice-regrowth, which is a critical process for understanding the full response of the ice sheet to climate change during the Pleistocene (Pico et al., 2018; Dalton et al., 2019). In future work, we plan to incorporate other kinds of constraints, for example the total amount of time that a site was ice free or the need for ice to re-advance over the site as opposed to re-nucleating as a separate ice cap before coalescing with the rest of the ice sheet.

5 Conclusions

We present the calculation of sea-level potential as a novel method to constrain past ice-sheet geometry and its associated uncertainty within an ensemble-based ice-sheet modelling framework. Our results reveal regions of the GrIS which, when ice-free, are associated with a narrow range of ice-sheet geometries. Future programs to collect samples from beneath the ice-sheet margins and interior, including the U.S. National Science Foundation-funded GreenDrill (Briner et al., 2022), and Green2Ice, an ERC Synergy Grant funded by the European Union, in combination with our results, may provide novel constraints on paleo sea level contributions from the GrIS. Moreover, this modelling approach can be used to inform future drilling, as regions we identify represent locations where information about past ice-free conditions can be most directly translated into information that can inform efforts to adapt to sea-level change. Future work to expand the ensemble of past ice-sheet geometries is required to further develop the usefulness of this method. In particular, we expect that improved knowledge of the past spatial mass balance patterns and relationships between temperature and precipitation change will have the greatest impact on our results. Our results reveal distinct spatial patterns of ice-sheet sensitivity to different physical processes, which can provide input to scientific communities working on understanding these processes in space and time; for highly vulnerable regions of GrIS such as southwest Greenland, the spatial climatology pattern, treatment of precipitation-lapse rate, and ice-sheet starting geometry influence how the ice-sheet deglaciates. Reducing uncertainty in sea-level projections for both paleo and future scenarios will require efforts to better constrain these parameters. More precise knowledge of e.g. past climate forcing can be readily incorporated into our experimental design so that sea-level fingerprinting and local sea-level impact predictions are informed by the most relevant sources of paleo-data as Earth's climate continues to warm.