Growth and decay of the Iceland Ice Sheet through the last glacial cycle

Constraining the dynamic evolution of past ice sheets is critical for unravelling their responses to external forcing and feedbacks over long timescales. This is particularly true in the context of marine ice sheet collapse, as this is one of the largest sources of uncertainty for future sea-level rise projections. The Iceland Ice Sheet (IIS) provides an empirically constrained case study for investigating such an instability, having retreated from a predominantly marine-based ice sheet to isolated mountain ice caps during the last deglaciation. However, previous reconstructions of the IIS have been limited by either sparse data or a restricted exploration of model parameter space, lacking a robust quantification of uncertainties. Here, we address this gap by performing a truncated history matching of the last glacial cycle of the IIS. We use the Glacial Systems Model (GSM) constrained by a curated set of geochronological data to generate an envelope of not-ruled-out-yet(NROY) ice sheet histories.

Our results indicate that numerous asynchronous ice streams effectively drain ice from the interior to the margins, resulting in an extensive yet relatively thin ice sheet. During its local Last Glacial Maximum (23.6–20.9 ka), the IIS reaches the continental shelf edge in most sectors with a total volume of 0.41 to 0.76 metres equivalent sea level (m e.s.l.). In the most extreme NROY glaciation scenarios, our model reveals an ice bridge connecting the Iceland and Greenland ice over Denmark Strait.

We find that accelerated ice discharge (at the grounding line) dominates mass loss during deglaciation. This acceleration is primarily driven by atmospheric warming through a cascade of mechanisms: surface meltwater induces hydrofracturing, leading to both ice shelf disintegration and tidewater calving, which in turn reduces buttressing and triggers rapid ice stream acceleration. The critical role of hydrofracturing in enabling model capture of deglacial data constraints is shown by explicit sensitivity experiments. This thereby supports inclusion of hydrofracturing for modelling of ongoing ice sheet response to climate change.

1 Introduction

Improving constraints on the potential of marine ice sheet collapse is crucial for projecting future sea-level rise (Robel et al., 2019). However, this large-scale instability, not yet observed in the modern era, and its associated drivers remain poorly known (Robel et al., 2019; Pattyn and Morlighem, 2020; Coulon et al., 2024). The Iceland Ice Sheet (IIS) was a predominantly marine-based ice sheet with multiple ice streams during the last glacial cycle that ended with near complete deglaciation by the mid-Holocene (Patton et al., 2017; Benediktsson et al., 2022b). It thereby offers a case study to investigate such a termination complementary to those examining terminations in Antarctica (e.g. Pollard et al., 2015) and northern Eurasia (e.g. Sejrup et al., 2022).

Reconstructing IIS evolution is challenging due to limited empirical constraints and large physical uncertainties. The latter includes: (i) oceanic and atmospheric conditions with strong dependence on the Atlantic Meridional Overturning Circulation (AMOC, as shown in Xiao et al., 2017); (ii) controls on grounding line stability and migration; (iii) calving and submarine melt processes; and (iv) ice stream drainage tied to the strong geothermal activity (Flóvenz and Saemundsson, 1993; Hjartarson, 2015). The warm/cold oscillations were controlled by the variability of the cold East Icelandic Current and the relatively warm Irminger Current derived from the North Atlantic Current (Xiao et al., 2017). The local climate was therefore contingent on AMOC strength. Based on dated marine limit shorelines, Norđdahl and Ingólfsson (2015) suggest that the marine-based part of the ice sheet collapsed during the Bølling–Allerød, speculating the key role of calving and bedrock topography in the deglaciation. The tectonic and magmatic legacy of Iceland translates to a strong geothermal heat flux (Flóvenz and Saemundsson, 1993; Hjartarson, 2015), conducive to basal melting. Geomorphologists identified streamlined subglacial bedforms, suggesting the presence of fast-flowing paleo ice streams (Bourgeois et al., 2000; Spagnolo and Clark, 2009; Principato et al., 2016; Benediktsson et al., 2022b, a). Digital elevation models reveal troughs terminating at moraines (Fig. 1, GEBCO Bathymetric Compilation Group 2023, 2023), with the troughs likely carved by ice streams during past glaciations (Spagnolo and Clark, 2009; Benediktsson et al., 2022a).

.

Figure 1 Bathymetric and topographic map of Iceland from GEBCO DEM with the present-day glacial situation (Farinotti et al., 2019).

Despite recognition of these physical processes, empirical studies can only speculate about their relative contribution. They cannot quantify the associated uncertainties, nor capture potential nonlinearities and synergistic effects of multiple forcings acting in concert. Moreover, although proxy records and glaciological imprints in mapped bedforms have enabled reconstructions of past IIS margins (Spagnolo and Clark, 2009; Benediktsson et al., 2022b, 2023a, b, c), these reconstructions lack confidence for time intervals and areas with sparse data coverage. Consequently, significant knowledge gaps persist in this empirical framework particularly for the pre-Last Glacial Maximum (LGM) configuration, the southern and eastern margins, the ice volume and sea level contribution during the LGM, as well as the drivers and controls of deglaciation (Benediktsson et al., 2022b, 2023a).

Glaciological models offer a potential tool to bridge these gaps. Previous modelling studies (Hubbard et al., 2006; Hubbard, 2006; Patton et al., 2017) have advanced our understanding of IIS dynamics, simulating a relatively fast ice sheet with a low aspect ratio and strong sensitivity to atmospheric temperature and geothermal conditions. However, these studies explored very limited parts of the parameter space, typically hand-tuning a small number of model parameters (<6) to minimize visual misfit to empirical ice margin reconstructions. As is the typical case for paleo ice sheet modelling, this approach along with very limited uncertainty assessment makes interpretation of the extent to which model results bear on the actual physical system hard to assess (given that uncertainty assessment is the specification of the relationship between data, model, and the physical system).

To address this challenge of inferential clarity, we present a truncated history matching analysis where a large ensemble of glacial cycle simulations is constrained against a curated set of paleo constraints. History matching identifies a set of simulations that are not ruled out by available data constraints given a robust uncertainty analysis. This entails full accounting for both model and data uncertainty. However, given the limited number of data constraints for the IIS, we elected to skip the deep sampling component of history matching, thus our “truncated history matching” moniker. We use the thermo-mechanically coupled glaciological 3D Glacial Systems Model (GSM, Tarasov et al., 2025) configured with 35 ensemble parameters to partially account for uncertainties in climate, basal drag, and marine ice processes.

Our results address several key research questions about the IIS. We investigate the configuration, dynamics and evolution of the IIS throughout the last glacial cycle (Sect. 3.2). We then disentangle the drivers and controls of its subsequent deglaciation (Sect. 3.2.3).

2 Methods

We describe below the glaciological model and our methodology for constraining it against paleo data.

.

Figure 2 Flow chart illustrating the main steps to calibrate paleo-ice sheets. Blue boxes indicate the model processes, orange boxes represent the constraint data, gray boxes denote the history matching loop. Rounded shapes indicate the function while squared shapes the input/output.

2.1 The ice sheet model

The Glacial Systems Model (GSM) is a 3-D thermo-mechanically coupled hybrid Shallow Ice Approximation (SIA)/Shallow Shelf Approximation (SSA) ice-sheet model that resolves ice streams, ice shelves, and grounding line migration. The GSM includes a permafrost resolving bed-thermal model. It incorporates fully coupled visco-elastic glacio-isostatic adjustment that accounts for the changing load of the adjacent Greenland Ice Sheet. The GSM uses a combined positive degree day (PDD) and novel positive temperature short-wave parameterization to account for the direct impact of changing orbital forcing on surface melt. The GSM also incorporates the marine ice cliff instability for tidewater ice margins as well as hydro-fracturing (Weertman, 1973; Van der Veen, 1998; Buck, 2023) following the approach of Pollard et al. (2015). As for any subgrid process, the large scale difference between the GSM grid resolution and actual terminal crevasse size reduces the need for accurate process representation as long as critical dependencies are captured within associated ensemble parameters. As such, hydro-fracturing is parametrized as a quadratic contribution from surface runoff to surface crevasse crack propagation. In previous modelling of Antarctica (Pollard et al., 2015), the choice of quadratic dependence was made to enable present-day preservation of ice shelves while facilitating marine ice collapse of the West Antarctic ice sheet during the Eemian interglacial to better capture the inferred sea level highstand. As this choice imposes a potential structural uncertainty on the model, the GSM includes a flag for linear dependence on surface runoff. Calving and hydro-fracturing process uncertainties are further addressed by three relevant ensemble parameters (coefficients for surface runoff contribution, calving rate, and summer temperature threshold for calving). A complete description of the GSM is available in Tarasov et al. (2025).

We use a lat-lon grid with resolution of 0.0625° × 0.125° ($\sim \mathrm{7}\times \mathrm{6}$ km). Our choice of grid resolution represents a trade-off between the need to resolve key physical processes and computational efficiency. The chosen grid resolution explicitly resolves the major ice streams, while keeping the computational cost manageable for running a large ensemble of simulations.

For climate forcing, the GSM uses a combination of an asynchronously coupled 2D energy balance climate model (EBM) and both yearly and monthly resolved glacial indices to interpolate between LGM monthly temperature, precipitation and evaporation climatologies (Paleo-Modelling Intercomparison Project 3, Braconnot et al., 2012) and present-day climatologies (ECMWF ERA5 reanalysis, van Wessem et al., 2018). We use the North Greenland Ice Core Project (NGRIP) ice core record (North Greenland Ice Core Project members, 2004), transient simulation output from a climate model (Geng et al., 2025, Fig. 3), and monthly mean EBM air temperature for the 60 to 40° southern latitude band for the glacial indices.

Figure 3 Glacial indices derived from smoothed temperature anomalies (50 years running mean, interpolated to 50 year time steps) from the NGRIP ice core (North Greenland Ice Core Project members, 2004) and seasonal variations from the LCice model (February, April, and July, Geng et al., 2025). Raw signals are shown in translucent.

The GSM is forced by the ocean temperature field of the TraCE-21 ka deglacial simulation (Liu et al., 2009) performed with version 3 of the Community Climate System Model (CCSM3). Prior to 22 ka, the ocean forcing is extrapolated from the TRACE chronology using a glacial index, similarly to the atmospheric forcing. To partly address uncertainties in the ocean forcing, an ensemble parameter adjusts interglacial to glacial phasing of the index.

The bedrock topography is derived by subtracting the observed ice thickness field (Farinotti et al., 2019) from the GEBCO surface elevation DEM (GEBCO Bathymetric Compilation Group 2023, 2023).

Hank and Tarasov (2024) have demonstrated that the surge pattern over the Hudson Strait is very sensitive to the geothermal heat flux. However the 4 km deep geothermal heat flux under Iceland is not well-constrained. For this work, two geothermal heat flux fields are subject to weighted average according to an ensemble parameter (Fig. S1 in the Supplement). The two fields are generated by interpolating boreholes from Flóvenz and Saemundsson (1993) and Hjartarson (2015).

Our sediment distribution is the result of two merged fields: the NOAA offshore dataset (Straume et al., 2019) and a parameterized terrestrial field with 200 m of sediment at slope divergences and 0 m at convergences in high-resolution bedrock topography. We upscale the latter to the GSM resolution to preserve the high-resolution information. Seismic surveys indicate thick sediment layers in the valleys, likely deposited over the course of previous glaciations (Black et al., 2004).

2.2 Paleo-constraints

To evaluate model plausibility, we compiled a comprehensive IcelandICE database of geological constraints that capture the spatial and temporal evolution of the IIS throughout the last glacial cycle (orange boxes in Fig. 2). Our database is a compilation and refinement of published geochronological data (Norđdahl and Pétursson, 2005; Licciardi et al., 2007; Benediktsson et al., 2022b, 2023a, b, c), supplemented by geomorphological evidence (Spagnolo and Clark, 2009; Benediktsson et al., 2022b, a) (Fig. 4). The database is available on the online repository https://doi.org/10.5281/zenodo.19453940 (Goffin et al., 2026).

The geochronological constraints are classified into two categories. One category is observations of past ice extent (paleoEXT) which consists of marine radiocarbon ages. The other is past ice thickness (paleoH) which consists of terrestrial cosmogenic and radiocarbon ages (Fig. 4). The radiocarbon ages are calibrated to calendar years. We introduce a multi-tier data quality assessment. Only data with rigorous age calibration, clear interpretation, and providing clear constraint value are assigned to tier 1. Most of the past ice thickness data was assigned to tier 1. However, outliers (i.e., ages that are anomalously anti-phasing with other surrounding data) were systematically ranged into tier 4 (e.g. Sandfell cosmogenic age which is too old).

Data with robust calibration ages and high constraining value but subject to widely different, yet justified, interpretations were classified as tier 1A. These data were explicitly tested for impact of constraint on the simulations within the context of the interpretative debate. Tier-2 data supplement Tier 1 data, providing finer details of the deglacial ice-sheet-thinning history. Tier-3 data comprise lower-quality constraints, marked by low-confidence in uncertainty specifications and ambiguous interpretation.

.

Figure 4 Our Iceland constraint database summary plot. The topographic map was generated using GEBCO DEM and the geomorphological imprints of the IIS on the Iceland shelf are based on Spagnolo and Clark (2009), Benediktsson et al. (2022b, a)

Past ice extent (paleoEXT) constraints provide ages for the deglaciation of marine ice. Such data are from marine sedimentary cores and require a strong understanding of the stratigraphic units. For instance, some age may reflect the mixing with older sediments introducing ambiguity in the interpretation. Therefore, only data with a clear stratigraphic history are classified as tier 1. Since most of the IIS was marine-based at the LGM, tier-1 marine data play the most critical role in constraining the key deglacial changes. We compile the main radiocarbon ages from Benediktsson et al. (2022b, 2023a, b, c) while ensuring that the dates from the source studies have acceptable uncertainties. In the Reykjafjarđaráll-Húnaflóadjúp trough on the Northern shelf (∼20.2° W, 66.9° N), a date of 15.0±0.6 ka from post-glacial muds overlying a diamicton constrains the timing of ice retreat from the trough (Andrews and Helgadóttir, 2003). On the western shelf, a deglaciation age of 15.3±0.9 ka in the Jökuldjúp trough [24.21° W, 64.29° N] is based on an extrapolation down to ice-contact sediments from three radiocarbon dates in the overlying marine mud (Jennings et al., 2000; Norđdahl and Pétursson, 2005; Benediktsson et al., 2023a). Both of these data are classified as tier 1 and provide the most powerful constraint for the ice sheet retreat over the continental shelf. Tier 1A paleoExt data introduce important complexities and highlight areas under debate. For instance, in the Djupall trough on the Northwest shelf (∼24.0° W, 66.6° N), radiocarbon dates (calibrated) of 22.3±1.0 ka and 36.5±1.5 ka suggest that glacier ice was absent during the LGM, and that the area was potentially occupied by shorefast sea ice instead (Geirsdóttir et al., 2002; Andrews et al., 2018). This is in contrast with geomorphological evidence suggesting more extensive ice cover. Furthermore, older Tier 1A dates from the Reykjafjarđaráll-Húnaflóadjúp trough indicate that ice had culminated there before 28 ka and persisted until the onset of deglaciation around 15 ka (Andrews et al., 2000).

Past ice thickness is constrained by cosmogenic nuclide exposure ages from mountain summits, nunataks, as well as radiocarbon ages from sites that constrain local deglaciation timing (paleoH). We compiled and recalculated the cosmogenic ³He ages of table mountain summits from Licciardi et al. (2007) using version 3 of the online exposure age calculator formerly known as the CRONUS-Earth calculator (Balco et al., 2008). We use the same locally calibrated Iceland ³He production rate from Licciardi et al. (2006), time-varying geomagnetic scaling LSDn (Lifton et al., 2014), and the same isostatic adjustments applied in Licciardi et al. (2006). Exposure ages are calculated with no erosion, in line with observed preservation of original lava flow surface morphologies at the sample sites. In a scenario wherein a subaerial lava flow erupts to form a table mountain summit shield that is continuously exposed following emplacement, the exposure age should be equivalent to the post-glacial time of eruption, and inherited cosmogenic isotopes are a non-issue (Licciardi et al., 2007). Geomorphic evidence suggesting glaciation of various Icelandic table mountain summits has been described in some previous studies (e.g. van Bemmelen and Rutten, 1955; Eason et al., 2015). These infer a more complex exposure scenario involving initial emplacement of the table mountain summit shield, quickly followed by ice coverage when the glacier surface healed over the edifice after eruptive activity ceased, then subsequent re-exposure when the ice sheet ultimately thinned to reveal the summit. If this latter case applies, exposure ages on such table mountain summits would be dating the timing of deglaciation rather than post-glacial lava emplacement age. We do not have the means of distinguishing between these two exposure history scenarios with existing isotopic data, but we note that in both scenarios the ³He exposure ages provide a deglacial summit age and therefore provide a pertinent empirical constraint for ice thickness history. No adjustments are made for snow cover, which if applied would result in older exposure ages for the table mountain summits, but the time-integrated local snow history at the sample sites is poorly known. The recalculated table mountain ³He ages are within ∼1 % of the ages originally reported in Licciardi et al. (2007). Key Tier 1 paleoH constraints include cosmogenic ³He ages from several table mountains such as Herđubreiđ (10.4±1.0 ka, 1522 m), Búrfell (10.8±0.6 ka, 949 m), and Bláfjall (14.3±0.8 ka, 1087 m, Licciardi et al., 2007). These sites provide crucial lower limits on ice thickness in the interior of Iceland. Additionally, radiocarbon dates from marine shells, such as those from Vopnafjörđur (11.0±0.3 ka, 50 m) and Hvalvík (14.7±0.4 ka, 40 m), constrain the timing of local deglaciation.

While not directly dated, extensive geomorphological evidence provides further context for understanding the IIS configuration and dynamics. Streamlined subglacial ridges and cross-shelf bathymetric troughs indicate the presence of fast-flowing ice streams that drained the ice sheet interior (Spagnolo and Clark, 2009; Principato et al., 2016; Benediktsson et al., 2022a). Terminal moraines identified at or near the continental shelf edge in most sectors support an extensive LGM ice sheet (Patton et al., 2017). Additionally, the radial pattern of these landforms suggests a primary ice divide near the center of Iceland (Spagnolo and Clark, 2009).

The spatial distribution of the geological constraints is uneven. As shown in Fig. 4, data are scarce over the southern and eastern regions of the Icelandic continental shelf. Providing a best-estimate age and age uncertainties for the LGM extent of the IIS is thus challenging and somewhat subject to interpretation of chronological and stratigraphic data (Benediktsson et al., 2022b, 2023a). However, based on the existing empirical data from different sites on the western, northwestern, and northern continental shelf, the LGM peak extent can be securely bracketed between maximum ages of 45.1 and 39.8 ka and minimum ages of 22.3, 16.2, and 15.3 ka (median probability ages based on calibrated age ranges at the 95.4 % ±2σ level), with a best-estimate age of 23–28 ka (Benediktsson et al., 2022b, 2023a).

The peak extent during the Younger Dryas is comparably well constrained both spatially and chronologically. The empirical data suggest that it was asynchronous across the ice sheet, most probably owing to internal glacier dynamics, topographic control, and local or regional climate conditions. However, where the age of the Younger Dryas ice margin is well constrained, it broadly falls between 12.8 and 12.1 ka (Benediktsson et al., 2023b). In other areas, the relative age of the Younger Dryas ice margin can be established from cross-cutting relationships with Preboreal deposits that are typically dated between 11.1 and 11.6 ka, yielding a wider range of 12.5±3.0 ka.

Our database is designed to constrain last glacial cycle model simulations. For a compilation with comprehensive Holocene coverage, we refer the reader to Harning et al. (2026).

2.3 Experimental design

We produced a large ensemble of simulations based on Latin Hypercube sampling to efficiently explore the model parameter space, ensuring broad coverage of parameter interactions (blue boxes in Fig. 2). Each of these simulations is characterized by a parameter vector that specifies the value of all ensemble parameters. Sensitivity tests were used to identify which potential GSM ensemble parameters had significant impact on various model run statistics. As a result of these tests, the GSM was configured with 35 ensemble parameters (9 for ice dynamical/basal drag, 15 for climate forcing and surface mass balance, 5 for ice calving/subshelf melt, 1 for deep geothermal heat flux, and 2 for basal hydrology). Each simulation is initialized with the present-day ice thickness field (Farinotti et al., 2019) and runs from 122 to 0 ka (2000 CE). Each simulation was then evaluated against the data (see History Matching in the next Sect. 2.4).

In addition, we conducted sensitivity experiments to isolate the role of hydrofracturing and absolute sea level change in deglaciation, as well as to partially assess the impact of uncertainties in basal hydrology, the parametrized dependence of hydrofracturing on grid cell runoff, and those due to lateral variations in earth rheology. Detail on these experiments and results thereof are in Sect. 3.2.4 below.

Overall, in this study, we conducted a total of ∼7500 simulations: 6000 for ensembles, and 1500 for sensitivity experiments.

2.4 History Matching

History matching identifies a set of model chronologies that are not ruled out given available data constraints and robust uncertainty analysis (gray boxes in Fig. 2). As such, it aims to “bracket reality” as opposed to the much more difficult task of determining a meaningful most likely chronology (Tarasov and Goldstein, 2021). Each simulation is scored against key metrics and subsequently deemed implausible or not via multiple sievings. To be meaningful, this must account for both model and data uncertainties. The resulting products are sets of NROY (not-ruled-out-yet) chronologies.

The scoring uses all tier 1, 1A and 2 constraint data (past ice thickness score and past ice extent score). This implausibility metric is a conceptual distance between the model output and the associated data constraint normalized by model and data uncertainties (Eq. 1). The past ice thickness score is calculated by comparing the paleoH observation to the corresponding model grid cell (± spatial uncertainty ≈5 km) and time slice (± time uncertainty =500 years). We used a similar approach for the tier 1A past ice marine extent constraint, where the score is based on whether the model grid cell was covered by grounded ice, floating ice, or was ice-free, depending on the sample type. The samples can either infer the presence of open marine conditions (OMC), sub-ice shelf (SIS), or proximal to grounding line (PGL). As per the two tier 1 past marine extent constraints, the past deglaciation score is calculated by evaluating the capability of the model to expand over the data points during Marine Isotope Stage 2 (MIS 2) and subsequently deglaciate within the calibrated 2σ range data ages. The implausibility values I_i for each datum and each parameter vector is calculated as follows:

$$\begin{array}{}\text{(1)}& I_i^{\mathrm{2}}(\mathit{x}\mathbf{,}\mathit{c}\mathit{m})={\displaystyle \frac{{\left(M_i(\mathit{x}\mathbf{,}\mathit{c}\mathit{m})-d_i\left(x\right)-{\mathit{ϵ}}_{\text{total}}\right)}^{\mathrm{2}}}{{\mathit{\sigma }}_{\text{struct}}^{\mathrm{2}}+{\mathit{\sigma }}_{\text{obs}}^{\mathrm{2}}}}\end{array}$$

where M_i(x,cm) is the model output at location and time x with parameter vector “cm”, d_i(x) the datum value, ϵ_total the mean bias error term, and σ_struct, and σ_obs the structural and observational standard deviations respectively.

In this study, we applied three different sievings to produce 3 distinct subsets.

The first and most critical sieving consists of ruling out simulations if they are deemed implausible according to the top tier-1 data. In particular, we ruled out simulations if either of the following thresholds were met:

• the implausibility with respect to the two top tier marine constraints is above 3, a rejection threshold based on the 3σ rule (Pukelsheim, 1994). The two marine constraints are the datum at [24.21° W, 64.29° N] (Andrews et al., 2000) and the datum at [20.23° W, 68.86° N] (Jennings et al., 2000; Benediktsson et al., 2022b), with inferred 2σ deglaciation time intervals of [13.52–16.51] ka and [13.41–17.16] ka respectively.

• the implausibility with respect to the top tier 1 past ice thickness data is above 3.

The resulting NROY_tier1 subset is the main product of the present study that is used to bracket the last glacial cycle evolution and serves as the main basis for our analysis (see Sect. 3.2).

The second sieving provides an IC_{data_tier1} subset based on the same sieving as NROY_tier1 but without accounting for the structural uncertainty. This gives a narrower set of chronologies that are nominally consistent with the data. The third sieving consists of ruling out simulations deemed implausible by considering both the tier-1 and the contentious tier-1A data, within model and data uncertainty. The same rejection threshold of 3 was applied. Here, the resulting product is a set of NROY chronologies: NROY_tier1-1A.

2.5 Ice discharge assessment

We used two different approaches to estimate the mass of ice discharged into the ocean. Our first approach calculates the ice discharge at the grounding line using the component method: we subtract the total net mass balance from the surface mass balance of grounded ice, D= SMB−MB.

The second approach calculates the ice discharge over gates placed within bathymetric troughs and valleys. That is, where ice streams are likely to appear. To estimate the relative role of mass loss from ice streaming, we also computed the ice stream discharge normalized by the ice sheet volume. A similar approach was used by Stokes et al. (2016) for the Laurentide Ice Sheet. Throughout this paper, we use the term “ice streams” to encompass both ice streams sensu stricto (fast-flowing regions bounded by slower-moving ice) and outlet glaciers (fast-flowing ice confined within bedrock valleys and fjords).

3 Results

This results section first evaluates model performance against geochronological constraints (Sect. 3.1), then presents the bracketed chronology of IIS evolution through the last glacial cycle based on the main NROY subset (Sect. 3.2).

3.1 Model fits to data constraints and scoring

Here the 2000 final large ensemble members are compared to the tier-1 and tier-1A data using 3σ implausibility thresholds.

The NROY_tier1 (tier 1) subset has 107 members. The IC_data-tier1 subset (only constraint data uncertainties) has 44 members. And the NROY_tier1-1A subset (including the contradictory 1A constraint data) has only 23 members.

While there is a relatively low misfit score with respect to paleoH data (most NROY_tier1 scores <2σ), marine past ice extent data scoring is the limiting factor (2σ< NROY_tier1 score <3.5σ). This discrepancy primarily stems from the poor score related to the Djúpáll data-point (see SITID 2004 in paleoExt dataset – see Sect. 3.1.1).

3.1.1 Past marine ice extent

The marine ice extent data provide critical constraints for evaluating model performance across Iceland's continental shelf sectors. Here we focus on the data-model comparison for the two NROY subsets NROY_tier1 and NROY_tier1-1A. Each of the NROY members are compared against the tier 1 and tier 1A past ice extent data which can infer proximal-to-GL (PGL) conditions, sub-ice shelf (SIS) conditions, or open marine conditions (OMC). Overall, while the tier 1 observations are bracketed by the NROY sub-ensembles (site 2003, 2006), tier-1A data fitting presents a greater challenge (2002, 2004, Fig. 5a).

.

Figure 5 Past ice extent (a) and past ice thickness (b) misfit scores. The data ID locations are shown in Fig. S2.

In the Reykjafjarđaráll-Húnaflóadjúp trough located on the Northern shelf (∼20.2° W, 66.9° N), the NROY sub-ensembles perform well with the tier-1 observation (2003) but most of the NROY members struggle to bracket the oldest datum (2002). The latter PGL constraint (2002) suggests ice free conditions at 28.6 ka, whereas most NROY simulations indicate sustained glaciation around 30 ka or earlier.

In the Northwest Djúpáll trough (∼24.0° W, 66.6° N), none of the NROY members bracket the data point (2004), suggesting that the trough was ice free at ∼22.3 ka. However, the Djúpáll data are subject to methodological limitations including core depth, the number and location of cores, datable material, reservoir corrections, and interpretations of different stratigraphic units, all of which introduce ambiguity. Andrews et al. (2002) claim that the exact location of the LGM margin is uncertain due to insufficient data, but that the Djúpáll trough was ice-free from 36 ka based on dated foraminifera in diamict and ice-rafted debris (IRD) in overlying sediment. Subsequent studies (Geirsdóttir, 2004; Chesley, 2005; Principato et al., 2006; Quillmann et al., 2012) propose revising this chronology, indicating that Djúpáll trough was most likely ice-filled during the LGM. The older dates of ∼36 ka may reflect reworking of older sediments rather than indicating a restricted LGM extent in Djúpáll. Consequently, more recent studies depict the ice sheet margin at the continental shelf edge at the LGM (e.g. Geirsdóttir et al., 2009; Geirsdóttir, 2011; Andrews et al., 2018).

Moreover, this conflicts with multiple nearby observations (e.g., SITID 2003, 2006, and 2007), which indicate persistent grounded ice across the continental shelf from ∼25–15 ka. The scenario implied by this data (2004), with the Djúpáll trough ice-free by ∼22.3 ka while grounded ice persisting further south in the Jökuldjúp trough [24.21° W, 64.29° N] until $\sim \mathrm{15.3}\pm \mathrm{0.9}$ ka (Jennings et al., 2000) [SITID = 2006], is glaciologically dubious.

On the western shelf, in the Jökuldjúp trough (2006), the NROY models, particularly the NROY_tier1-1A subset, perform relatively well. The larger model-data mismatch for some NROY_tier1 simulations reflects the large spread in modelled grounding-line retreat within the Jökuldjúp trough (see details in Sect. 3.2). This spread is a direct result of parameter variations within the NROY_tier1 ensemble. In particular, the highly non-linear grounding-line dynamics in this deep marine embayment are very sensitive to parameters controlling basal sliding and calving rates. Consequently, while the NROY_tier1-1A sub-ensemble captures the overall deglacial pattern, individual simulations can diverge from the timing of the constraint.

3.1.2 Past ice thickness

Cosmogenic exposure ages, and to a lesser extent radiocarbon ages, constrain past ice thickness enabling robust comparison between modelled and empirical ice sheet thinning (Fig. 5b). We evaluate NROY_tier1 and NROY_tier1-1A ensembles against each tier-1 past ice thickness observation using two complementary scoring components. First, “presence of ice” scoring gives zero misfit if modelled grid-cell ice thickness exceeds the empirical value prior to the data age. Second, ”absence of ice” scoring gives zero misfit if modelled ice thickness is below the inferred ice thickness at the data age (Fig. S3). Both misfit scores increase proportionally with temporal offset between modelled and empirical values.

While the NROY subsets successfully bracket all paleo ice thickness observations, certain NROY simulations yield higher misfit scores for specific data points. For example, at Herđubreiđ (SITID = 1005), the empirical age of deglaciation ranges between 11.4 and 9.4 ka, but in one of the NROY simulations (run identification number nn555), deglaciation occurs earlier at 15 ka, resulting in a high misfit score for the “presence of ice” (score = 19). However, given the early deglaciation, the “absence of ice” condition after 9.4 ka is satisfied and the associated score is zero. Conversely, further north at Bláfjall (SITID =1009), the empirical deglaciation window is 15.1–13.5 ka, whereas in one of the NROY simulations (nn599), deglaciation occurs later, at 10 ka. This leads to a low misfit score for the “presence of ice” but a high misfit (score = 36) for the “absence of ice”, as the model grid cell should not have ice after 13.5 ka. Persistent data-model mismatches are likely in good part attributable to inadequate grid resolution for resolving the complex topography of Iceland.

3.2 Last glacial cycle IIS chronology

The evolution of the IIS volume and area through the last glacial cycle for the full initial ensemble and more-data constrained NROY_tier1 sub-ensemble are presented in Figs. 6a and S4. Out of 2000 ensemble members, 107 are NROY. Ensemble variance is not well correlated with ice volume (Fig. S5). However, once the ice sheet expands beyond the present-day coastline (i.e., once the ice area threshold of $\sim \mathrm{100}\times {\mathrm{10}}^{\mathrm{3}}$ km² is crossed), the ensemble variance grows significantly, reflecting regionally divergent rates of grounded ice advance (Fig. 6b). Depending on the ensemble parameter vector choice, the grounding line fluctuates considerably in light of the gentle shelf slope. Furthermore, the relative variance, i.e., the variance normalized by the mean ice volume, peaks during the deglaciation (Fig. S5).

Figure 6 (a) Time series of IIS volume during the last glacial cycle. Initial full ensemble (gray shaded area) and sieved not-ruled-out-yet (NROY_tier1) sub-ensemble (blue distribution) with median NROY_tier1 (black curve). (b) Density distribution of grounded ice within the NROY_tier1 sub-ensemble at different time slices. White contours on the 25, and 12 ka panels show empirical ice-margin reconstructions for 25.5±2.5 ka (Benediktsson et al., 2022b, 2023a), and 12.5±3.0 (Benediktsson et al., 2023b), respectively.

3.2.1 Pre-LGM

For the Last Interglacial spin-up state, we initialize the model with the present-day ice thickness field (Farinotti et al., 2019). At the onset of the last glacial cycle, ice in our NROY_tier1 ensemble expands from two main centers: the eastern part of the central highlands from the existing ice caps, particularly Vatnajokull and, a smaller glaciation over the Northwest peninsula (Fig. 6b). By ∼90 ka the two glaciated ice bodies merge into an ice sheet covering the entire terrestrial landmass. The ice sheet is characterized by a double dome configuration, established on the foundations of the initial two ice bodies. Subsequently, at ∼83 ka, a rapid retreat separates the two domes.

Between MIS 5a and the end of MIS 3 (80–30 ka), the grounding line advances further and fluctuates across the continental shelf, constrained at its maximum by the shelf edge. Within the NROY_tier1 sub-ensemble, the main variations of extent occur over the Northern and Southeastern continental shelf. For similar ice volumes, the direction of ice expansion can vary significantly, contingent on the parameter choice. For some parameter vectors, the absence of ice over the northern shelf (70 ka in Fig. 6b), can be attributed to the depth of bathymetric troughs (Fig. 1). Conversely, simulations show that grounded ice can extend across this area when anchored on shallow parts of the continental shelf such as the Kolbeinsey Ridge (∼ 19° W, 67° N), part of the Mid-Atlantic Ridge. Furthermore, the Kolbeinsey Ridge can also act as a pinning point that buttresses and stabilizes the ice shelf, thereby affecting grounding line dynamics.

Owing to the underlying soft till, numerous troughs, and high geothermal heat flux, the ice sheet develops multiple ice streams across the continental shelf. Most of these ice streams activate and deactivate independently (Fig. S7), with basal velocities periodically dropping to zero, indicating complete shutdowns (Fig. S6a, d–f). Given this independence and the relatively stable, cold glacial climate (with insufficient oceanic or atmospheric warmth to induce substantial surface or sub-shelf melt (Fig. S8e), and thus little change in ice shelf buttressing) we infer that ice-stream activation is primarily driven by internal thermodynamic “binge-purge” cycling rather than external forcing. Preferential drainage pathways from the ice-sheet interior to the margin vary between simulations and through time within individual simulations, often favouring one trough over adjacent ones. In contrast, some ice streams remain continuously active in the Southwestern and northern shelf sectors (Fig. S6b, c).

The main limitation predating the LGM is lack of constraints. Consequently, no scenarios of ice expansion pathway can be ruled out and the ice margin position remains highly uncertain.

3.2.2 The Last Glacial Maximum

At MIS 2, the IIS advances further onto the shelf, specifically along the western margin. The ice sheet maintains its double-dome structure, with the primary dome centred over the present-day Vatnajokull region and a secondary dome over the Northwest peninsula (Fig. 7). The simulated NROY_tier1 ice volume spans from 0.76 to 0.41 m eustatic sea level equivalent (m e.s.l.), with a median of 0.6 m e.s.l. at the local LGM (23.6 to 20.9 ka). The large variance in LGM volume is primarily due to differences in the extent of the grounding line across the Southwestern shelf, the potential link with Greenland as well as variations in the overall ice thickness.

During the global LGM, summer air temperatures vary within our NROY_tier1 sub-ensemble, but remain below freezing point, with little to no surface melt (Fig. S8e). The IIS expansion is therefore mainly constrained by the continental shelf edge and reaches it in most sectors. This is supported by the presence of numerous moraines all along the shelf edge (Fig. 4). The exception is in the southwest where the grounding line position is more variable across NROY models (Fig. 6b). Our NROY_tier1 simulations at 25 ka bracket the empirical ice-margin reconstructions for 25.5±2.5 ka (Benediktsson et al., 2022b, 2023a).

.

Figure 7 Ice sheet surface elevation (gray shadings) and bedrock elevation (colour shadings) of the mean NROY_tier1 at 20 ka. Ice surface-height contours are delineated at 450 m intervals, revealing the double-dome configuration. Indicated ice shelf extent (white filled area) was determined from the ensemble-mean ice shelf geometry, excluding grid cells where mean ice thickness was <25 m. The minimum (dashed line) and maximum (solid line) NROY_tier1 ice (grounded + floating) extent are shown in orange.

The ice sheet is bounded by ice shelves, particularly along the northern margin. Some NROY_tier1 simulations exhibit an ice bridge over Denmark Strait connecting the Northeastern Iceland ice shelf to the Greenland Ice Sheet (Fig. 7). The latter reached the shelf edge over its Southeastern margin during the LGM (Andrews et al., 2000; Dowdeswell et al., 2010; Funder et al., 2011; Vasskog et al., 2015). Even though not located within our model domain, core MD99-2323 at a depth of 1062 m, situated over the southern Denmark Strait, shows evidence of bioturbation during the LGM (Andrews et al., 2017), which has been interpreted as indicating open water conditions at that location. However, Riddle et al. (2007) documented benthic life under the Amery Ice Shelf, demonstrating that bioturbation does not preclude ice shelf cover. Andrews et al. (2017) further infer that the Denmark Strait served as the primary gateway for iceberg passage, as the IRD sediments in core MD99-2323 are sourced from northern and northeastern Greenland, the Arctic Basin, and the European Arctic, though a contribution from the clockwise Irminger Current cannot be ruled out. Moreover, paleoceanographic data of the northern North Atlantic at the LGM reflect relatively warm conditions with perennial sea-ice cover confined to the northeast Greenland margin area (Kucera et al., 2005). Conversely, there is evidence for a reduction in the strength of the Western Boundary Undercurrent (WBUC) around Southern Greenland during glacial stages (Müller-Michaelis and Uenzelmann-Neben, 2014). This is likely due to reduced Denmark Strait Overflow (Andrews et al., 2018), potentially attributed to the presence of an ice shelf. In some simulations, grounded ice extends fully across Denmark Strait (Fig. 6b).

The presence of a thick, soft basal sediment layer in marine sectors and valley bottoms as well as high geothermal heat flux facilitate fast basal sliding. This partly explains why all simulations, regardless of the parameter choice, depict a highly dynamic ice sheet with numerous fast-flowing ice streams that effectively drain the ice from the interior to the margin (Fig. 8a). Another factor is the presence of numerous topographic troughs that over time guide ice stream drainage of the ice sheet (Fig. 1 and 8b). Most if not all of these troughs are likely the result of the two-way feedback between topographic excavation by subglacial erosion and topographic steering and confinement of ice streaming. Over time and location there are varying degrees of topographic steering and confinement with some ice streams clearly topographically confined, and others crossing over topographic ridges (Fig. 8b).

.

Figure 8 Basal velocity for one physically self-consistent NROY_tier1 run at 20 ka (run identification number ic10593). (a) The basal velocity field is shown in colour. The background shows the bed topography in greyscale. The streamline ridges are also shown (orange lines) (Benediktsson et al., 2022a). (b) Bedrock topography is shown in colour. Overlain are basal velocity contours at 200 m yr⁻¹ intervals (black lines), the grounding line (dotted white line), and the ice margin (solid white line).

Given the temporal variation in ice stream configuration for a given simulation as well as the variation between ensemble members (Fig. S7), no single simulation at a given time is likely to be fully consistent with relevant geological data such as streamline ridges (Fig. 8a). A further challenge in such comparison is the poor age control on such features, with only indirect inferences that they were formed during the local LGM (Hoppe, 1968; Bourgeois et al., 2000; Spagnolo and Clark, 2009; Principato et al., 2016; Benediktsson et al., 2022a).

Such effective drainage results in an extensive ice sheet with a small mean ice thickness on the order of 811±71 m as per the NROY sub-ensemble. This is lower than crude estimates based on volume-area scaling (using a scaling exponent of 1.23 from Paterson, 1994, and an area of 250 000 km² giving a mean ice thickness =1026 m), and even lower than past models for Iceland (939 m for Hubbard et al., 2006, and 1172 m for Patton et al., 2017).

3.2.3 The last deglaciation

The IIS loses ∼95 % of its mass during the last glacial-interglacial transition (∼20.9 to ∼9 ka). The deglacial mass loss is dominated by calving until ∼11.5 ka. The IIS deglaciation history is characterized by a two-step evolution with five distinct phases: three of moderate ice retreat or even re-advance (I, III, and V in Fig. 9), and two of intervals of rapid ice loss (II and IV in Fig. 9).

.

Figure 9 Climate forcing (input) and mass balance components (output) of the IIS over the last deglaciation (19–9 ka) for the NROY_tier1 subset. (a) Global mean eustatic sea level reconstruction. (b) Ocean temperature at 191 m depth from TraCE-21ka deglacial simulation (Liu et al., 2009) for southwestern (red) and northern (blue) grid cells. Gray line represents the mean across all grid cells with ±1σ uncertainty (shading). (c) Glacial indices derived from smoothed temperature anomalies (50 years running mean, interpolated to 50 year time steps) from the NGRIP ice core (North Greenland Ice Core Project members, 2004) and seasonal variations from the LCice model (February, April, and July, Geng et al., 2025). (d) Mass balance components of grounded ice including net mass balance (blue), surface mass balance (orange dashed), and ice discharge at the grounding line (yellow dashed), with shading indicating uncertainty (±1σ). (e) Total ice mass (blue, left axis) and normalized mass fluxes (right axis) including surface melt, snow accumulation, calving, and sub-shelf melt. Vertical dashed lines represent the deglaciation phase limits.

First, between the end of the local LGM and 14.6 ka (Phase I), the grounding line retreats moderately across the South-western shelf (Fig. 6b). The mass loss is tied to a slight offset between ice discharge and surface mass balance (Fig. 9d). This is driven by sea level rise alone (Fig. 9a) as the climate is still relatively cold (Fig. 9b, c), consistent with TraCE-21000 simulations showing negative surface temperature anomalies relative to 22 ka until 15 ka over the North Atlantic region (Liu et al., 2009; He, 2011). However, we cannot rule out that the lack of an apparent submarine melt role in initiating deglaciation may be due to limitations of the relatively coarse resolution GCM used to generate the TRACE ocean temperature field. Between the end of the local LGM and ∼14.6 ka, the atmospheric temperatures surrounding Iceland for the NROY_tier1 simulations remain below the threshold necessary for substantial surface melt (Fig. 9c, e).

Between 14.6 and 14 ka (II), ice discharge significantly exceeds the surface mass balance (Fig. 9d) and the marine-based ice sheet collapses. The ice sheet transitions from a marine to a terrestrial setting (Fig. 6b). The mass loss is primarily driven by atmospheric warming during the Bølling–Allerød. We hypothesize that atmospheric warming drives increased ice discharge through hydrofracturing. Accordingly, we conducted a sensitivity experiment without hydrofracturing (see Sect. 3.2.4 below). An increase in surface melting, particularly at lower altitudes (i.e., on the ice shelves over the northern margin and tidewater glaciers over the southern margin; Fig. S9), intensifies surface crevasse growth through hydrofracturing, thereby increasing calving flux (Fig. 9e). The resulting ice shelf disintegration (Fig. S10) together with the tidewater front retreat reduces back stress, which would enhance ice flux across the grounding line and thereby at the calving front. The disintegration of the ice shelves leads to less total sub-shelf melt (Fig. 9e) in spite of an increase in ocean temperature forcing. The marine ice sheet instability (MISI) potentially further accelerates the grounding-line retreat (Ritz et al., 2015). Based on several NROY runs, it appears that the Kolbeinsey ridge (∼19° W, 67° N in Fig. 1) may act as a pinning point that buttresses the ice, maintaining continuous and sustained drainage of the ice streams until ∼14.6 ka. Once this anchoring stops, the ice sheet retreats rapidly across the mid-outer shelf. Such a combination of feedbacks results in the steepest rate of mass loss throughout the deglaciation (Fig. 9e). This aligns with previous geological inferences of marine ice sheet collapse during the Bolling-Allerød (Norđdahl and Ingólfsson, 2015).

Following the phase of rapid mass loss, most of the NROY_tier1 runs have ice regrowth from 14 to 11.8 ka (III). This is driven by both an increase in accumulation and a substantial reduction in surface melt, which simultaneously increases the surface mass balance (Fig. 9d) and reduces calving by inhibiting hydrofracturing (Fig. 9e). Simulations from the whole ensemble exhibit a range of behaviours over this time interval: there is either a strong decrease in deglaciation rate, stabilization, or a re-advance of the IIS (Fig. S4). Simulations with smaller ice volumes generally exhibit proper re-advance, while those with larger ice volumes experience a decreased deglaciation. At 12 ka, our minimum NROY_tier1 extent closely matches the empirical ice-margin reconstructions for 12.5±3.0 ka (Fig. 6b, Benediktsson et al., 2023b).

By 11.8 ka (IV), the ice sheet shrinks rapidly. The ice loss is primarily driven by atmospheric warming at the end of the Younger Dryas. This results in a large increase in surface runoff which gives rise to a sharp decrease in surface mass balance as well as a brief increase in calving (Fig. 9d, e). The increased surface melt increases hydrofracturing and therefore calving until it is offset by decreasing the marine fraction of the ice margin.

At 11.5 ka, the IIS margin generally coincides with the present-day coastline, aligning with past empirical inferences based on shorelines and moraines (Norđdahl and Pétursson, 2005; Sigfúsdóttir et al., 2018; Benediktsson et al., 2023b). From 11.5 to 9 ka (V), the mass change becomes progressively dominated by the surface mass balance. By 10.7 ka, the magnitude of surface ablation offsets the snow accumulation, resulting in a negative total surface mass balance. The thinning accelerates as surface melt extends to higher elevations (Fig. S11), strengthening the melt-elevation feedback. Once the main ice dome is gone, the residual ice margin retreats gradually to the higher elevations until reaching a configuration comparable to present-day conditions at ∼9 ka. This is consistent with empirical evidence indicating that ice extent was similar to present prior to 9 ka (Benediktsson et al., 2024).

Our above analysis demonstrates that the IIS deglaciation is dominated by ice discharge during the Phase II rapid retreat episode. Large volumes of ice discharge are associated with the acceleration and expansion of ice streams which drain interior ice to the margin (Fig. 10). This acceleration is initiated by atmospheric warming, which triggers both ice shelf disintegration and tidewater glacier retreat via hydrofracturing (see Sect. 3.2.4 on hydrofracturing below). With the elimination of ice shelf buttressing and reduced backstress at tidewater margins, ice streams accelerate. Here, ice stream acceleration is therefore triggered by ice shelf removal and tidewater glacier retreat (see Sect. 3.2.4 on hydrofracturing sensitivity), in correspondence with recent observations of ice stream and outlet glacier acceleration in Antarctica and Greenland (Rignot et al., 2004; Wuite et al., 2015; Bell and Seroussi, 2020; King et al., 2020), rather than internal processes as simulated for the Laurentide Ice Sheet during last glacial cycle (Hank and Tarasov, 2024) or IIS throughout MIS 2 (this study). Therefore, contrary to inferences for the Laurentide ice sheet (Stokes et al., 2016), we conclude that ice streams can play a critical role in changing ice sheet mass balance. However, once the IIS fully retreats to non-marine margins (phase V), ice streams largely shut down (Fig. S7), consistent with geomorphological signatures (Aradóttir et al., 2023, 2024).

Figure 10 Ice stream discharge (a) and ice stream discharge normalized by the ice sheet volume (b) of the NROY_tier1 subset. The lines represent the means and the shaded areas the 2σ range. Vertical dashed lines represent the deglaciation phase limits as per Fig. 9. Ice stream discharge was calculated over gates placed within bathymetric troughs and valleys.

On the basis of this data-model comparison, we conclude that the IIS deglaciation is most likely caused by a non-linear series of mechanisms. That is, once a temperature threshold is crossed, surface runoff on the ice shelf increases calving via hydrofracturing, leading to both ice shelf disintegration and tidewater glacier front retreat. The resulting elimination of ice shelf buttressing and reduced backstress at tidewater margins trigger ice stream acceleration and thereby increased ice discharge into the ocean. Hydrofracturing thereby provides a coupling between atmospheric forcing, specifically surface warming and rain, and enhanced grounding line discharge of ice.

3.2.4 Sensitivity experiments

Hydrofracturing is a key component of both ice shelf and tidewater calving in the GSM (c.f. Eqs. 31 and 32 in Tarasov et al., 2025). To decipher the relative role of hydrofracturing in the IIS deglaciation, we conducted three deglacial sensitivity experiments on the full NROY set of ensemble parameter vectors. In turn, these experiments were: without hydrofracturing, with hydrofracturing only for tidewater margins, and with hydrofracturing only for ice shelves.

The NROY ensemble without hydrofracturing has a much reduced rate of net mass loss during the 14.6 to 14.0 ka phase II interval over which the NROY ensemble with hydro-fracturing loses more than 60 % of its mass (Fig. 11). The experiment therefore shows that hydrofracturing plays a key role in the non-linear response of simulated ice sheet retreat to atmospheric warming for at least our NROY set of parameter vectors. With hydrofracturing, the NROY ensemble has an approximate doubling in ice discharge across the grounding line at ∼14.5 ka, which is also the start of the interval of strongest mass loss from ice calving (Fig. S8). Without hydrofracturing, there is no increase in ice discharge across the grounding line during the whole phase II interval (Fig. S12). From 15 to 11.6 ka, the no hydrofracturing ensemble has moderate grounding line retreat. Subsequently, at ∼11.6 ka, once a certain threshold of temperature is crossed, the rate of mass loss increases significantly. Here, contrary to the reference set-up, surface mass balance completely dominates the mass loss, with subshelf melt and, to an even lesser extent, calving playing secondary roles (Fig. S12). The limited calving and sub-shelf melt are not enough to trigger major ice shelf disintegration or tidewater calving (Figs. S12, S13b). Therefore, ice stream acceleration is substantially reduced compared to the reference experiment with hydrofracturing (Fig. S14d).

.

Figure 11 Sea level contribution relative to 15 ka from experiments with hydrofracturing (black), without hydrofracturing (red), without hydrofracturing only for tidewater margins (orange), and without hydrofracturing only for ice shelf (blue) based on the 107 NROY_tier1 parameter vectors. The lines represent the medians and the shaded areas show the 2σ ranges.

The experiments with hydrofracturing only for tidewater margins, and with hydrofracturing only for ice shelves indicate that both mechanisms contribute comparably to total hydrofracturing (Fig. 11). Thus, regardless of whether margins terminate in ice shelves or grounded tidewater glaciers, hydrofracturing remains critical in deglaciation. The normalized ice stream discharge in the no tidewater hydrofracturing experiment is slightly higher than in the experiment without any hydrofracturing (Fig. S14f), while the no ice shelf hydrofracturing experiment is comparable to the experiment without any hydrofracturing (Fig. S14h).

Given the model's high sensitivity to hydrofracturing and its potential reliance on this process to capture critical tier-1 marine constraints, there is the possibility that our analysis is skewed by this reliance. To address this, we reran the whole 2000 member initial ensemble with hydrofracturing turned off over the whole glacial cycle. We then sieved it against the tier-1 data using the same 3σ implausibility threshold applied to NROY_tier1. Only 2 members meet the NROY criterion, and neither are nominally consistent with the data. In these instances, the model fails to reproduce the last deglaciation dynamics of the IIS, in particular the collapse of the marine ice sheet suggested by the data. As such, for at least the GSM, hydrofracturing is required in order for the IIS simulation to be consistent with our tier I constraints. However, other potential IIS deglaciation drivers not considered here, such as episodic volcanic intervals (Lin et al., 2022), cannot be ruled out.

These results highlight the potentially critical role that hydrofracturing can play in marine ice sheet deglaciation. Pollard et al. (2015) draw similar conclusions in their reconstruction of the Antarctic penultimate deglaciation. We therefore suggest that recent projections (e.g. Payne et al., 2021; Klose et al., 2024) may underestimate the relative contribution of atmosphere in driving future marine ice sheet mass loss unless they incorporate hydro-fracturing (DeConto and Pollard, 2016; DeConto et al., 2021; Coulon et al., 2024). This is particularly likely as 60±10 % of Antarctic ice shelves that buttress upstream ice has been inferred to be prone to hydrofracturing if surface crevasses are filled with water (Lai et al., 2020).

To partly address structural uncertainties in the representation of hydrofracturing (i.e. uncertainties not covered by ensemble parameters), the NROY_tier1 parameter vector set was re-run with linear (as opposed to quadratic) dependence of hydro-fracturing on surface runoff for the whole glacial cycle. The linear coefficient was scale matched to that of the default quadratic dependence for a 0.5 m yr⁻¹ surface runoff. This scale choice approximately matches mean ensemble deglacial runoff over the ice shelves during the 15 to 14 ka interval. The impact of this structural change is nearly unbiased on the ensemble (Fig. S17), with less than 5 % impact on 15 ka grounded ice area for most simulations. The impact for three other metrics (20 ka ice volume, ice thickness misfit score, and marine ice extent misfit score) is larger over the ensemble, but in each case it is still less than 5 % for the majority of simulations.

Another potential driver of marine ice sheet collapse is sea level rise due to increasing ocean volume during deglaciation as well as associated changes in the gravitational field. To isolate the relative impact of this, we carried out a sensitivity experiment repeating the 15 to 0 ka deglacial interval of the simulations with the geoid held constant at its 15 ka value ($\approx -\mathrm{100}$ m around Iceland relative to the Earth's center of mass). This was applied to the whole NROY_tier1 set of simulations. This change in forcing had minimal impact on the ensemble. A few simulations had slightly reduced mass loss rates (maximum difference at any time of ∼0.02 m e.s.l.), while the rest had almost no visually discernible response (Fig. S15).

The GSM visco-elastic GIA solver is restricted to a spherically symmetric earth rheology. However, there are likely significant lateral variations in the visco-elastic earth structure under Iceland. To bound potential uncertainties due to such variations, we reran the NROY_tier1 set with nominal stiff and soft bound earth rheology values. These rheologies were respectively assigned upper mantle viscosities of 5×10¹⁹ and 5×10¹⁸ Pa s along with respective lithospheric thicknesses of 96 and 46 km in accord with the previously inferred range of values for a deglacial timescale (Le Breton et al., 2010; Auriac et al., 2013). Simulations with the stiffer earth rheology tended to have a bit more 20 ka ice volume (ranging to 14 % larger, Fig. S16). In contrast, the stiff bound earth rheology gave less 15 ka grounded ice area (ranging to 10 % smaller) as well as larger “iceTot” ice thickness misfits against cosmogenic age constraints. The latter is also indicative of insufficient regional ice. The difference in half-space timescale for this rheology range (approximately 1050 years) is not much larger than the nominal ice sheet dynamical timescale during early to mid deglaciation (average grounded ice thickness/average accumulation ≈850 years). As such, a non-negligible impact is not unexpected.

Finally, to partly address uncertainties due to the impact of a simplified “leaky bucket” local sub-glacial hydrology model, the NROY_tier1 set was rerun with a 10 m (as opposed to the default 2 m) limit on subglacial water thickness. In comparison to the reference ensemble, differences in 15 ka grounded ice area, 20 ka ice volume, and IceTot score were all less than 5 % for all but a handful of ensemble members (Fig. S18). The paleoExt marine misfit score again showed the largest impact, but differences in individual scores were still less than 10 % for all but two of the NROY_tier1 parameter vectors. A much more detailed assessment of the impact of basal hydrology model complexity in the context of Hudson Strait scale ice stream surge cycling (Drew and Tarasov, 2023) found that most of the differences were within the range of what was already covered by parametric uncertainties for the simplistic basal hydrology model employed herein.

4 Conclusions

We have presented a truncated history matching analysis of the IIS evolution during the last glacial cycle. The GSM was configured with 35 ensemble parameters, covering climatological (including ocean), glaciological, and solid-Earth components. Out of 2000 simulations, 107 were not ruled out by the constraint data. Additional structural uncertainties related to basal hydrology, water collection into surface crevasses, and lateral variations in earth rheology were experimentally assessed, and only lateral variations in earth rheology were found to be non-negligible.

Throughout the last glacial period, the IIS fluctuated extensively across the continental shelf with numerous asynchronous ice streams. The ice sheet exhibited an extensive, shallow profile with a double-dome configuration. At the local LGM (23.6–20.9 ka), ice volumes range from 0.41 to 0.76 m e.s.l. NROY ensemble members include those with an ice bridge between the Greenland and Iceland ice shelves. For a few ensemble members, the ice is grounded across Denmark Strait. We know of no geological data that can clearly either support or refute such an ice bridge. High-resolution sea floor mapping and sediment coring across the shallower portions of the strait could potentially shed light on this question and should be explored in future studies.

The IIS remained a relatively stable ice sheet until the Bølling–Allerød, with only minor retreats driven by sea level rise. The rapid deglaciation (14.6–14 ka) was associated with the collapse of the marine-based ice sheet, and started only once discharge significantly outweighed surface mass balance. The mass loss was primarily driven by the atmosphere through disintegration of ice shelves and tidewater calving via hydrofracturing, which in turn induces enhanced ice streaming and thereby increased ice flux to the calving margins.

Contrary to past inferences for the Laurentide ice sheet (Stokes et al., 2016) which had extensive terrestrial ice streams as well as more confined marine-based ice streams, our data-constrained numerical experiments indicate that enhanced marine-based ice streaming can drive net mass loss when there is enough surface melt to induce hydro-fracturing.

On the basis of the demonstrated pivotal role of hydrofracturing in enabling the model to capture deglacial constraints, we conclude that hydrofracturing is likely a critical mechanism in marine ice sheet retreat (as previously concluded by Pollard et al., 2015, for Antarctica).