Modeled Greenland Ice Sheet evolution constrained by ice-core-derived Holocene elevation histories

The Cryosphere, 19, 3599–3622, 2025 https://doi.org/10.5194/tc-19-3599-2025 © Author(s) 2025. This work is distributed under the Creative Commons Attribution 4.0 License. Modeled Greenland Ice Sheet evolution constrained by ice-core-derived Holocene elevation histories Mikkel Langgaard Lauritzen1, Anne Solgaard3, Nicholas Mossor Rathmann1, Bo Møllesøe Vinther1, Aslak Grindsted1, Brice Noël4, Guðfinna Aðalgeirsdóttir2, and Christine Schøtt Hvidberg1 1Niels Bohr Institute, University of Copenhagen, Copenhagen, Denmark 2Institute of Earth Sciences, University of Iceland, Reykjavík, Iceland 3Geological Survey of Denmark and Greenland, Copenhagen, Denmark 4Laboratoire de Climatologie et Topoclimatologie, SPHERES, University of Liège, Liège, Belgium Correspondence: Mikkel Langgaard Lauritzen (mikk[email protected]) Received: 16 July 2024 – Discussion started: 24 July 2024 Revised: 17 June 2025 – Accepted: 17 June 2025 – Published: 10 September 2025 Abstract. During the Holocene, the Greenland Ice Sheet (GrIS) experienced substantial thinning, with some regions losing up to 600m of ice. Ice sheet reconstructions, paleoclimatic records, and geological evidence indicate that, during the Last Glacial Maximum, the GrIS extended far beyond its current boundaries and was connected with the Innuitian Ice Sheet (IIS) in the northwest. We investigate these long-term geometry changes and explore several possible factors driving those changes by using the Parallel Ice Sheet Model (PISM) to simulate the GrIS thinning throughout the Holocene period, from 11.7ka ago to the present. We perform an ensemble study of 841 model simulations in which key model parameters are systematically varied to determine the parameter values that, with quantified uncertainties, best reproduce the 11.7ka of surface-elevation records derived from ice cores, providing confidence in the modeled GrIS paleo evolution. We find that since the Holocene onset, 11.7 ka ago, the GrIS mass loss has contributed 5.3 ±0.3 m to the mean global sea-level rise, which is consistent with the ice-core-derived thinning curves spanning the time when the GrIS and the Innuitian Ice Sheet were bridged. Our results suggest that the GrIS is still responding to these past changes, having raised the sea level by 23±26 mmSLEka−1 in the last 500 years. Our results have implications for future ice sheet evolution, which should account for this long-term transient trend. 1 Introduction During the Last Glacial Maximum (LGM), approximately 20 ka ago, Earth was covered by large ice sheets, including the Laurentide, Fennoscandian, Innuitian, and Greenlandic ice sheets, and the global mean sea level was 125–134m lower than today (Lambeck et al., 2014; Yokoyama et al., 2018). Geological evidence suggests that the Greenland Ice Sheet (GrIS) extended to the continental shelf and was connected to the Innuitian Ice Sheet (IIS) at the Nares Strait (England et al., 2006). Toward the end of the last glacial period, the Bølling– Allerød interstadial brought abrupt warming to the Northern Hemisphere 14.7 ka ago, followed by cooling in the Younger Dryas stadial 12.9ka ago (Rasmussen et al., 2006). The Holocene interglacial began 11.7ka ago, bringing temperatures that were locally up to 15°C warmer in Greenland (Andersen et al., 2004). However, temperature reconstructions vary by several degrees, which is crucial for simulating the GrIS Holocene evolution (Nielsen et al., 2018). Following the Holocene Thermal Maximum, 6–9 ka ago, Greenland temperatures have shown a long-term decreasing trend (Vinther et al., 2009) but anthropogenic forcing has since reversed the course of natural temperature change, resulting in a global increase in temperatures since pre-industrial times (Eyring et al., 2021). Accurately modeling the historical evolution of the GrIS is essential for evaluating and calibrating ice sheet models. Ice sheet models respond to climate change over Published by Copernicus Publications on behalf of the European Geosciences Union. 3600 M. L. Lauritzen et al.: Modeled Greenland Ice Sheet evolution constrained by Holocene elevation histories a range of timescales and are rarely in a steady state (e.g., Lauritzen et al., 2023). However, several ice sheet model studies have overlooked a calibration of their temporal evolution and only focused on the evolution of ice temperature, neglecting other delayed responses, such as bedrock dynamics. For example, the ISMIP6 protocol does not require calibration (Nowicki et al., 2020), while most of the ISMIP6 ensemble simulations underestimate the observed IMBIE consensus mass loss from the GrIS (The IMBIE Team, 2020; Aschwanden et al., 2021). Recent advances have addressed this by calibrating an ice sheet model to satellitebased gravimetry-derived mass-loss data of the GrIS (Aschwanden and Brinkerhoff, 2022) but the satellite-based calibration data period only covers 22 years at the time of writing, and there is no guarantee that it gives a sensible longterm response. Calibrating the model to align with present-day observations of ice thickness and velocities risks capturing only the present-day state while being on a wrong state trajectory; that is, neglecting the long-term memory of the ice sheet and the response of the bedrock to past changes in ice load. These differences in past trajectories affect the projected future mass loss in this century, as demonstrated by A ¯ dalgeirsdóttir et al. (2014). To simulate time periods before the satellite era, ice sheet modeling must rely on proxy data from paleo-climatic records and ice extent markers for constraining and validating the long-term transient response of the ice sheet (state trajectory) over these considerably longer timescales. Past temperatures can be inferred from oxygen isotope measurements. When water evaporates from the oceans and precipitates over the GrIS, a temperature-dependent fractionation process alters the ratio of oxygen isotopes in the water – a relationship first used by Dansgaard et al. (1969) to infer past temperatures from oxygen isotope measurements at Camp Century (CC). Vinther et al. (2009) used this temperature dependence to derive a GrIS-wide oxygen isotope signal by assuming that the Renland and Agassiz (see Fig. 1) icecore sites are located within restricted ice domes where ice thickness remains constant. This GrIS-wide oxygen isotope signal was then subtracted from the oxygen isotope signals at CC, NGRIP, GRIP, and Dye 3 (see Fig. 1) to derive local surface-elevation histories, after correcting for upstream effects. These surface-elevation histories provide constraints for modeling the GrIS throughout the Holocene, offering valuable insights into the ice sheet’s response to past climate changes and helping to improve the robustness of model predictions. Previous studies have attempted to model elevation changes derived from ice cores. Notably, Lecavalier et al. (2017) modeled Holocene surface-elevation changes at CC using temperature anomalies from the Agassiz ice cores, suggesting that early Holocene temperatures were 7°C higher than today. However, they did not account for the buttressing effect of the IIS, a key driver of thinning (MacGregor et al., Figure 1. Model domain showing the present-day bedrock topography from Morlighem (2022), Jakobsson et al. (2020), and the GEBCO Bathymetric Compilation Group (2023) with the presentday ice cover from Morlighem (2022) and RGI Consortium (2023) shown in white. The ice-core sites discussed in the text (CC, NGRIP, GRIP, Dye 3, Renland, and Agassiz) are shown together with the glacier catchment basins (NW, CW, SW, SE, CE, NE, NO) from Mouginot and Rignot (2019) and extended out to the exclusive economic zone of Greenland (Flanders Marine Institute (VLIZ), Belgium, 2023) and constitute our extended continental shelf (ECS) domain, see text. 2016), and focused only on relative elevation changes, without reconstructing absolute elevation history. More recently, Tabone et al. (2024) successfully modeled elevation changes at the GRIP site, attributing them to the onset of the Northeast Greenland Ice Stream (NEGIS). In this study, we use the Parallel Ice Sheet Model (PISM) to model the long-term transient response of the GrIS to past climatic changes and the collapse of the IIS bridge during the Holocene. By varying 20 influential model parameters in an ensemble of 841 members, we show that it is possible to model the ice-core-derived elevation histories rather than just the thinning if the grounding line can advance to the continental shelf and the GrIS can connect to the IIS. We use this setup to constrain the model parameters for the ensemble and to estimate the GrIS long-term evolution with quantified uncertainties. Using the calibrated model, we investigate the Holocene ice sheet mass loss and assess the ongoing longterm response of the modeled GrIS and bedrock dynamics. The Cryosphere, 19, 3599–3622, 2025 https://doi.org/10.5194/tc-19-3599-2025 M. L. Lauritzen et al.: Modeled Greenland Ice Sheet evolution constrained by Holocene elevation histories 3601 Table 1. The 20 parameters that are varied in our ensemble of simulations. The temperature reconstruction is sampled discretely. The estimated parameter values are given as the mean plus or minus the standard deviation of the posterior PDFs, except for the temperature reconstruction, which is given as the mode of the posterior PDFs. Parameter Description Spin-up Range Estimate Combined CC NGRIP GRIP Dye 3 Atmosphere 1T Temperature reconstruction 3 1–5 1 1 4 2 3 fsPDD param. snow (mmK−1d−1) 5.04 5.7–8.9 6.5 ±0.7 7.1 ±0.8 7 ±1 7.5 ±0.7 7.8 ±1.0 fiPDD param. ice (mmK−1d−1) 12.5 7–10 7.7 ±0.6 8.5 ±0.8 8.6 ±0.8 8.7 ±0.9 8.4 ±0.7 0Atmospheric lapse rate (Kkm−1) 5 4–9 5.4 ±0.7 7 ±1 6 ±1 6 ±1 6 ±1 ω↓Southern precip. scaling (%K−1) 5 0–4.5 2 ±1 2 ±1 2 ±1 2 ±1 2 ±1 ω↑Northern precip. scaling (%K−1) 7 0–9 2 ±1 3 ±2 2 ±1 4 ±2 3 ±2 Ocean Hcr Threshold for thickness calving (m) 50 50–150 96 ±17 87 ±25 97 ±30 112 ±30 106 ±30 σmax Characteristic stress (MPa) 1 0.8–1.2 0.92 ±0.09 1.03 ±0.10 1.0 ±0.1 1.0 ±0.1 1.0 ±0.1 ˙mo ↓Melt rate south of 71° N (ma−1) 300–500 394 ±29 409 ±65 400 ±56 391 ±56 391 ±50 ˙mo ↑Melt rate north of 80° N (ma−1) 10–30 20 ±6 19 ±6 21 ±6 19 ±6 19 ±5 τOcean melt onset (ka) 4–8 5.6 ±0.6 5.5 ±0.7 6 ±1 6 ±1 6 ±1 1τ Ocean melt ramp-up time (ka) 0–2 1.1 ±0.5 1.1 ±0.6 0.9 ±0.6 1.0 ±0.6 0.7 ±0.6 Dynamics nSSA Creep exponent for the SSA (1) 3.3 3.2–3.4 3.28±0.04 3.35 ±0.04 3.33 ±0.04 3.25 ±0.04 3.23 ±0.03 ESIA Enhancement factor for the SIA (1) 3 2.5–3.3 3.0 ±0.2 2.7 ±0.2 2.9 ±0.2 3.1 ±0.2 3.1 ±0.2 qBasal sliding power coefficient (1) 0.8 0.7–0.9 0.82 ±0.05 0.79 ±0.05 0.81 ±0.06 0.80 ±0.06 0.83 ±0.06 δEffective pressure parameter (%) 2 1.5–2.5 2.1 ±0.2 2.0 ±0.2 2.0 ±0.3 2.0 ±0.3 2.1 ±0.3 ϕmin Minimal till friction angle (°) 10 5–10 8 ±1 8 ±1 7 ±1 7 ±1 8 ±1 ϕmax Maximal till friction angle (°) 42 40–45 43±1 43 ±2 42 ±1 43 ±1 43 ±1 zmin Lower elevation cutoff (m) −700 −600 to −300 −421 ±60 −481 ±80 −443 ±84 −433 ±80 −447 ±77 zmax Upper elevation cutoff (m) 700 0–500 271 ±172 238 ±118 231 ±134 252 ±154 271 ±146 2 Model setup To model the Holocene evolution of the GrIS, we use the open-source Parallel Ice Sheet Model (PISM) version 2.1 (Bueler and Brown, 2009; Winkelmann et al., 2011) at 20 km resolution for the spin-up, refined to 10km for the last 20 ka. PISM is a three-dimensional thermomechanically coupled model that solves both the Shallow Ice Approximation (SIA) and the Shallow Shelf Approximation (SSA) in a hybrid scheme, capturing both slow-moving interior flow and fast flow in ice streams and outlet glaciers. At the ice–ocean boundary, PISM includes sub-grid parameterizations to model grounding-line advance and retreat (Gladstone et al., 2010). The model parameters, listed in Table 1, are varied in our ensemble simulations unless otherwise specified. 2.1 Model domain The model domain, shown in Fig. 1, spans 6.7 ×106km2, from the continental shelf in the east to the Canadian Arctic Archipelago in the west. A north polar stereographic projection with a standard parallel at 70° N and central longitude of −45° W (ESPG 3413) is used. The projection introduces distortions of up to +5% in the north and −11 % in the southwest relative to the central latitude and longitude. PISM uses a flat-Earth approximation, so volume is not conserved when transforming thickness between projections. Volumes and mass-loss rates are reported using the actual grid area. To partition mass between Greenland and Canada, we introduce an extended continental shelf (ECS) mask, corresponding to Greenland’s exclusive economic zone (Flanders Marine Institute (VLIZ), Belgium, 2023). This divides the two regions at the Nares Strait and Baffin Bay, extending to the continental shelf in the north, east, and south. For massloss partitioning when the grounding line advances, we extend the basins from Mouginot and Rignot (2019) to the ECS using nearest-neighbor extrapolation. The present-day bedrock topography over Greenland is from BedMachine v5 (Morlighem, 2022) and extended using data from IBCAO v4.2 (Jakobsson et al., 2020) and the GEBCO Bathymetric Compilation Group (2023) to cover the larger domain, in that order of preference, to get the best bedrock available. At the lateral boundary, a Dirichlet boundary condition of zero ice thickness is used; the influence of the majority of the Laurentide Ice Sheet is thereby neglected. The north and south are bounded by open ocean, while Iceland and Svalbard are just visible toward the east. At the base of a 2 km deep bedrock thermal layer, the thermal heat flux from Shapiro (2004) is applied constantly in time. https://doi.org/10.5194/tc-19-3599-2025 The Cryosphere, 19, 3599–3622, 2025 3602 M. L. Lauritzen et al.: Modeled Greenland Ice Sheet evolution constrained by Holocene elevation histories 2.2 Atmospheric forcing To model the surface mass balance (SMB), we apply a positive degree day (PDD) scheme to calculate the surface melting. This approach bases the SMB solely on temperature, T, and precipitation, P. In the PDD scheme, the surface melt rate is proportional to the extent to which the temperature exceeds the freezing point (e.g., Braithwaite, 1985), ˙ms∝max(0°C,T ), (1) where we use two constants of proportionality: one for snow, fs, and another for ice, fi. To force the PDD model, we use a 12-month reference climatology based on the multi-year monthly averages of temperature and precipitation for the period 1960–1989 from RACMO (Noël et al., 2015, 2018, 2019). Since our model domain is not covered by a single RACMO simulation, we combine different simulations (see Fig. A3). We merged three areas with precipitation data: Greenland, the northern Canadian Arctic Archipelago, and the southern Canadian Arctic Archipelago (Noël et al., 2018), treating areas outside these regions as having no precipitation. For temperature, we used RACMO2.3p2 at 5.5km for Greenland (Noël et al., 2019), combined with a broader 11 km simulation (Noël et al., 2015) for the rest of the area. The mean precipitation and summer temperatures for the resulting climatology are shown in Fig. A2. Following Nielsen et al. (2018), we account for paleo temperature changes by applying a spatially uniform timevarying temperature anomaly, 1T , along with a lapse rate adjustment, 0, which modifies the surface temperature based on deviations from the RACMO surface topography. The temperature reconstructions are shown in Fig. A1. Insolation changes are not included in this approach. Reconstructions 1 and 2 use the GRIP ice core with linear and quadratic transfer functions from Huybrechts (2002) and Johnsen et al. (1995), respectively. Reconstruction 3 uses the NGRIP core with the same transfer function from Huybrechts (2002), while Reconstruction 4 is the GrIS-wide reconstruction from Vinther et al. (2009) – the only one that accounts for elevation change. Reconstruction 5 is based on the NGRIP core with an isotope diffusion inversion scheme (Gkinis et al., 2014). The Holocene Thermal Maximum is only captured by Reconstructions 3, 4, and 5, while Reconstructions 1 and 2 suggest a more constant Holocene climate. Reconstruction 1 was used by the SeaRISE project (Bindschadler et al., 2013). Since the vapor pressure scales approximately exponentially with temperature in the Clausius–Clapeyron relation, we account for paleo precipitation changes by scaling the reference precipitation field with a time-dependent scaling factor, exp(ω(φ)1T (t)). Here, ωhas the latitude dependence ω(φ(x,y)) =       ω↓φ≤φp ↓ ω↓+φ−φp ↓ φp ↑−φp ↓ω↑−ω↓φp ↓≤φ≤φp ↑ ω↑φp ↑≤φ ,(2) where φis the latitude and φp ↓=60° N and φp ↑=75°N are chosen to cover most of Greenland; ω↓and ω↑are the southern and northern precipitation scaling parameters, respectively, which are varied in our ensemble. This approach allows for different precipitation histories in northern and southern Greenland, in contrast to the uniform scaling used in many previous modeling attempts (e.g., Nielsen et al., 2018). 2.3 Ocean forcing Following Aschwanden et al. (2019), we take the sub-shelf ocean melt to be separable in space and time: ˙mo(x,y,t) = ˙mo x(φ(x,y)) ˙mo t(t), (3) with the spatial dependence controlling the present-day melt rate given by ˙mo x(φ(x,y)) =     ˙mo ↓φ≤φo ↓ ˙mo ↓+φ−φo ↓ φo ↑−φo ↓˙mo ↑− ˙mo ↓φo ↓≤φ≤φo ↑ ˙mo ↑φo ↑≤φ ,(4) where φo ↓=71° N and φo ↑=80° N, following Aschwanden et al. (2019), while ˙mo ↓and ˙mo ↑are the upper and lower melt values, which we vary. To allow the formation of an ice bridge to Canada, the sub-shelf melt rate is scaled by ˙mo t(t) =     0t≤τ t−τ 1τ τ≤t≤τ+1τ 1τ+1τ ≤t ,(5) such that there is no ocean melt for times earlier than τ, while it increases to present-day values in the time 1τ, inspired by the rapid change in ocean temperatures found by Clark et al. (2020). In addition to sub-surface melt, ice is calved off at the ocean front at a rate that is proportional to the tensile von Mises stress and inversely proportional to a characteristic parameter, σmax (Morlighem et al., 2016). Additionally, all ice thinner than Hcr is calved off, and a eustatic sea-level forcing from Imbrie and McIntyre (2006) is applied, changing the ocean level by 130m in the last 19ka. 2.4 Ice dynamics The constitutive relation that relates the strain rate, ˙ij , to the stress, τij , in the ice sheet is ˙ij =EAτn−1 eτij ,(6) where Eis the enhancement factor, τeis the effective deviatoric stress, nis the creep exponent, and Ais the ice The Cryosphere, 19, 3599–3622, 2025 https://doi.org/10.5194/tc-19-3599-2025 M. L. Lauritzen et al.: Modeled Greenland Ice Sheet evolution constrained by Holocene elevation histories 3603 softness, which depends on temperature, pressure, and water content of the ice. The enhancement factor and the creep exponents are taken to be different for the SIA and the SSA. We use ESSA =1.3, nSIA =3 while varying ESIA and nSSA following Aschwanden and Brinkerhoff (2022). The numerical value of Ain Eq. (6) is not changed; only the units are adjusted. The basal sliding velocity ubin the SSA is related to the basal shear stress τbthrough the pseudo-plastic power law: τb= −tan(ϕ)Ntill ub uq th|ub|1−q,(7) where qis the sliding exponent and uth =100 ma−1is a characteristic speed. The till friction angle, ϕ, is parameterized as a continuous function of bedrock topography that increases linearly from ϕmin to ϕmax between zmin and zmax. The effective pressure on the till, Ntill, decreases exponentially with the water level in the till, Wtill (Bueler and van Pelt, 2015): Ntill =min(P0,˜ N0δP0 ˜ N0Wtill/W max till ).(8) It decreases to a fraction δof the overburden pressure, P0, when the water level reaches its maximum value, Wmax till =2 m. The parameter ˜ N0=5.6 ×108Pa represents the effective pressure that the till would have at zero water content if it were not capped by the overburden pressure. 2.5 Earth deformation and initialization The bedrock response to changes in ice load is given by the visco-elastic bed deformation model of Lingle and Clark (1985) and Bueler et al. (2007), with flexural rigidity D=5×1024 Nm and upper mantle viscosity η=1×1021 Pas. We initialize the GrIS by running the model from −100 to −20 ka (all times are relative to 2ka CE) at 20 km grid resolution using the parameters listed in Table 1 and with the initial bedrock topography taken to be the same as at the present day (Morlighem, 2022; Jakobsson et al., 2020; GEBCO Bathymetric Compilation Group, 2023). Since the modeled ice extent, and consequently the surface elevation, is sensitive to ocean melt and sea-level forcing, we apply an artificial correction to the bed topography to ensure that the present-day ocean mask closely aligns with observations. To achieve this, we iteratively adjust the bedrock at −20 ka (as illustrated in Fig. 2) so that the modeled bedrock topography at the end of the simulation better matches the observed present-day topography. A simulation with 20 km resolution is run from −20 ka to the present, and the deviation between the modeled and observed bedrock topography is used to update the initial bedrock according to b0 i+1=b0 i+Kbobs −b1 i,(9) Figure 2. Model ensemble experiment. The ice sheet is initialized at −100 ka using present-day geometry and run through the last glacial period at 20 km resolution. The bedrock is then iteratively updated at −20 ka to reduce the modeled present-day bedrock topography deviation. After finding a suitable bedrock topography, the ice sheet model is branched off at −20 ka, and an ensemble of simulations is run at 10 km resolution. While the xaxis depicts time, the yaxis is only used to reflect that the states differ. where bobs is the observed present-day topography (Morlighem, 2022; Jakobsson et al., 2020; GEBCO Bathymetric Compilation Group, 2023), b0 iis the modeled bedrock topography at −20 ka, and b1 iis the modeled bedrock topography at the present day. This iterative correction method is similar to the approach used by van Calcar et al. (2023), who employed a comparable scheme to improve present-day topography. The relaxation parameter K=0.7 is introduced to prevent overcompensation from any potential positive feedback associated with the updated bedrock, although such feedback may not be significant, given that the bedrock–mass balance feedback is likely to be negative. After 20 iterations, the root mean square error (RMSE) of the bedrock decreased from 77.4 to 3.3 m, as shown in Fig. 3. The impact of these iterations on the surface elevation is illustrated in Fig. A6. At −20 ka, the simulation is branched using the adjusted bedrock (b0) from the last step (shown in Fig. A8), and an ensemble of simulations is run at a 10km grid resolution until the present day, varying the 20 parameters listed in Table 1. 3 Bayesian inference To account for model uncertainty and assess the importance of model parameters, we run an ensemble of simulations from −20 ka to the present, varying the 20 parameters listed in Table 1. The dynamic parameters are based on those varied by Aschwanden and Brinkerhoff (2022), while the atmospheric and oceanic parameters are introduced in Sects. 2.2 and 2.3. To effectively sample the parameter space, 841 parameters are drawn using the second-order orthogonal Latin Hypercube Sampling (LHS) design (Tang, 1993). This ensures that all pairs of parameters are sampled uniformly and reduces the risk of clustering. For each of the four ice-core sites, we calculate the likelihood of observing the ice-core-derived elevation hishttps://doi.org/10.5194/tc-19-3599-2025 The Cryosphere, 19, 3599–3622, 2025 3604 M. L. Lauritzen et al.: Modeled Greenland Ice Sheet evolution constrained by Holocene elevation histories Figure 3. Iterative bedrock adjustment. (a) Modeled present-day bedrock elevation deviation, compared with Morlighem (2022), Jakobsson et al. (2020), and GEBCO Bathymetric Compilation Group (2023). (b)–(e) Zeroth, first, second, and 19th iterations of the modeled presentday bedrock elevation deviation from observed topography. tory, given the modeled elevation change, with model parameters m: ρhobs i|m∝Y j exp  −hi,j (m)−hobs i,j 2 2σ2 iβ  ,(10) where idenotes the ice-core site and jis the time step of the ice-core samples to which modeled elevations are interpolated; σiare the uncertainties from Vinther et al. (2009), derived from the spread of δ18O values in two parallel records from the Agassiz Ice Cap and the uncertainty in bedrock uplift at the Agassiz and Renland sites. The lapse rate uncertainty used to derive elevation changes is not included. Present-day and past elevations are weighted equally to avoid biasing the likelihoods toward the present configuration. Following Aschwanden and Brinkerhoff (2022), we introduce β=100 to account for autocorrelation in the uncertainties, reducing the number of degrees of freedom by a factor of 100, which corresponds to a decorrelation time of 2000 years. Additionally, we calculate the combined likelihood of the ice-core-derived elevation changes for all sites, which is proportional to the product of the site-specific likelihoods, assuming no spatial correlation between the drill sites: ρ(hobs|m)∝ 4 Y i=1 ρhobs i|m.(11) The parameters are sampled uniformly over the ranges specified in Table 1, which focuses on the volume of parameter space that shows the highest likelihood in an initial ensemble of simulations. The posterior joint probability density functions (PDFs) are then given by Bayes’s theorem: ρm|hobs i=ρhobs i|m ρhobs iρ(m), (12) where the prior distribution ρ(m)is taken to be uniform within the intervals listed in Table 1. From the five posteriors, we get five PDFs of ice sheet evolution through the Holocene, from which we estimate relevant observables, listed in Table 2. Unless stated otherwise, all model results are based on the combined posterior PDF. To evaluate the effectiveness of the sampling, we compute the effective sampling size for each of our five normalized posteriors: neff =1 Pkρmk|hobs i2,(13) where kdenotes the sample member. The effective sample size is the number of equally weighted samples that would yield the same variance of the mean as the weighted set of samples. If only one ensemble member has a non-zero likelihood, the effective sample size is 1. Conversely, if all members have the same likelihood, the effective sample size equals the actual sample size, namely 841. Although constraining the simulations to present-day observations would increase confidence in our modeled present-day state, we avoid doing so because this study focuses on the transient evolution of the GrIS and how the icecore-derived surface-elevation histories from Vinther et al. (2009) can be used to constrain the ice sheet’s evolution. 4 Results 4.1 Surface-elevation evolution Figure 4 shows the modeled and ice-core-derived surface elevations during the Holocene at the four ice-core sites. The site-specific modeled elevations closely match the ice-corederived reconstructions, reproducing the substantial thinning observed at CC and Dye 3, as well as the more moderate thinning at the interior sites GRIP and NGRIP, with RMSEs ranging from 12 to 53.6m. The combined elevation estimate is lower than the site-specific estimates at CC and Dye 3 at the onset of the Holocene, while it is too high at GRIP. The corresponding histories of bedrock elevation and ice thickness associated with these surface changes are shown in Fig. A7. The Cryosphere, 19, 3599–3622, 2025 https://doi.org/10.5194/tc-19-3599-2025 M. L. Lauritzen et al.: Modeled Greenland Ice Sheet evolution constrained by Holocene elevation histories 3605 Table 2. Estimates of key observables for the past and present of the GrIS, as well as observables for the simulations restricted to the present-day Greenland mask (Grl) and the ECS. All observables are calculated within the ECS mask at model resolution and do not include Canada. Observable Estimate Restricted Combined CC NGRIP GRIP Dye 3 Prior Grl ECS Elevation history RMSE CC (m) 88.1 53.6 116.8 154.3 170.7 119 223.4 165.9 NGRIP (m) 26.8 44.9 12 69.8 69.4 34.9 67.8 49.9 GRIP (m) 58.6 123 86 27.2 29.5 62 64.5 59.5 Dye 3 (m) 99.1 153.4 116.4 87.2 55.4 106.1 123.5 88.5 Present-day configuration Grounded ice volume (mSLE) 9.0 ±0.1 9.5 ±0.3 9.1 ±0.3 8.7 ±0.3 8.7 ±0.2 9.0 ±0.4 8.55 9.08 Grounded area (106km2) 1.99±0.02 2.00 ±0.03 1.94 ±0.05 1.95 ±0.05 1.93 ±0.04 1.95 ±0.06 1.82 1.98 Falsely grounded (106km2) 0.19±0.01 0.21 ±0.02 0.17 ±0.03 0.18 ±0.03 0.17 ±0.03 0.18 ±0.04 0.07 0.19 Missing grounded (106km2) 0.08±0.01 0.08 ±0.01 0.11 ±0.02 0.11 ±0.02 0.12 ±0.02 0.10 ±0.02 0.13 0.08 Ice thickness RMSE∗(m) 420.8 495.2 435.9 396.4 394.9 428 313.4 418.7 Bed topography RMSE (m) 27 42.2 19 18.5 17.5 20.2 58.3 58.3 Surface speed RMSE∗(ma−1) 84.4 80.7 79.5 82.7 82 79.1 95.3 96.3 Configuration at −12 ka Grounded ice volume (mSLE) 15.7 ±0.3 16.0 ±0.7 15.8 ±0.5 15.1 ±0.6 15.8 ±0.7 15.3 ±0.9 9.03 14.31 Grounded area (106km2) 2.96±0.03 2.93 ±0.03 2.93 ±0.04 2.96 ±0.05 3.01 ±0.05 2.93 ±0.05 2.08 2.87 dvdtlast 500 a (mmSLEka−1)−23 ±26 −31 ±27 −74 ±133 −52 ±106 −75 ±143 −60 ±112 −18.19 −70.81 Time of collapse (kab2k) 4.9 ±0.5 4.9 ±0.7 6 ±1 6 ±1 5 ±1 6 ±1 neff 9.28 41.36 151.87 99.67 28.56 841.00 ∗RMSEs are calculated within the present-day observed grounded mask. To illustrate the effect of allowing the ice sheet to advance beyond its present-day boundaries, we ran two additional simulations: one restricted from advancing beyond the present-day GrIS coast and another restricted from advancing beyond the ECS mask. Both simulations started from the unrestricted branch-off point at −20 ka and used the same parameters as the ensemble member with the highest combined likelihood, although these parameters may not necessarily be optimal for either of the restricted runs. The simulation restricted to the present-day GrIS coast could not reproduce the observed surface-elevation history at CC, NGRIP, and Dye 3, showing the importance of a dynamic grounding line. The simulation restricted to not advancing beyond the ECS performed better but also failed to reproduce the thinning at CC, showing the effect of including the Canadian Arctic Archipelago when modeling the GrIS Holocene history. The RMSEs associated with the four ice-core-derived elevation histories are listed in Table 2 for the five estimates. 4.2 Inferred parameters The ice-core-derived surface-elevation histories provide constraints on the model parameters, with the marginal PDFs for each of the five posteriors shown in Fig. 5. The degree to which individual parameters are constrained varies, and the estimated values are summarized in Table 1. Notably, the site-specific and combined estimates of the SIA enhancement factor, ESIA, differ substantially. The sitespecific estimate based on CC is 2.7 ±0.2, while that based on GRIP is 3.1 ±0.2. Similarly, the site-specific estimates of the SSA creep exponent, nSSA, differ, with higher values for CC and NGRIP than for GRIP and Dye 3. Among the five temperature reconstructions, Reconstruction 1, being the coldest throughout the Holocene, has the highest combined probability (61 %) and the highest sitespecific probability for CC (40%), where the largest surface thinning is observed. In contrast, Reconstruction 2 is more than 1 °C warmer during the early Holocene and performs worse than Reconstruction 1. Reconstructions 4 and 5, the warmest of the set, have near-zero probability at CC but perform better at other sites; notably, Reconstruction 4 has the highest site-specific probability for both NGRIP and Dye 3. The northern precipitation parameter, ω↑, is more constrained by the northern sites CC and NGRIP, where it has the most influence. Likewise, the southern precipitation parameter, ω↓, is most constrained by Dye 3. Both parameters are estimated to be 2 ±1 %K−1, which is substantially lower than the default of 7.3%K−1introduced by Huybrechts (2002), resulting in less accumulation in the warm periods of the Holocene and more accumulation in the cold glacial periods, where the ice sheet builds up. The onset of sub-shelf ocean melt is well constrained by CC, occurring at 5.6 ±0.6 ka before the present. https://doi.org/10.5194/tc-19-3599-2025 The Cryosphere, 19, 3599–3622, 2025 3606 M. L. Lauritzen et al.: Modeled Greenland Ice Sheet evolution constrained by Holocene elevation histories Figure 4. Observed and modeled surface elevation over the last 11.7ka for the ice-core sites Camp Century (CC), NGRIP, GRIP, and Dye 3. The blue lines are the ice-core-derived surface elevations from Vinther et al. (2009), and the blue envelopes denote 1 standard deviation. The orange solid lines show the combined PDF means, while the green solid lines indicate the site-specific PDF means, for each site. Shaded orange and green envelopes represent the corresponding 16th–84th percentile ranges. The dashed lines are the ensemble members with the highest combined likelihood (orange) and the highest likelihood for each site (green). The orange dash-dotted and dotted lines are simulations with the same parameters as the best ensemble member but restricted to the ECS (dash-dotted) and the present-day land margin of the GrIS (dotted). 4.3 Modeled Holocene evolution From the branch-off point at −20 ka until the onset of the Holocene (11.7 ka ago), the modeled ice sheet bridges the gap between Canada and Greenland across the Baffin Bay and the Nares Strait. Figure 6 shows the ice sheet configuration at −12, −9 ka, and the present day, while Fig. 7 presents the volume and area evolution from the branch-off point to the present. The model clearly responds to the change in resolution at the branch-off point, showing a positive drift in volume, though this shock appears to have stabilized before the start of the Holocene. By −12 ka, the ice sheet reaches its glacial maximum extent, grounding on the continental shelf and through the Nares Strait. At −12 ka, the GrIS has a modeled grounded area of 2.96 ±0.03 ×106km2within the ECS. This is 49% or 0.98 ±0.05 ×106km2larger than the present-day modeled area and it is 0.9% larger than the minimum LGM extent and 5.6 % smaller than the maximum LGM extent from Leger et al. (2024). Compared with the modeled present-day GrIS, the modeled grounded volume is 6.6±0.4 mSLE larger at −12 ka. Additionally, the grounded volume above flotation at −12 ka was 5.3 ±0.3 mSLE greater than at present. The modeled times of the last glaciation are shown in Fig. 8. Outside the ECS, the IIS and Laurentide Ice Sheet are cut off at the domain boundary with a Dirichlet boundary condition of zero thickness. This moves the ice divide at Baffin Island farther to the east than if it had been connected to a complete Laurentide Ice Sheet. Together, they have a grounded area of 1.20 ±0.03 ×106km2and a grounded volume of 5.0 ±0.2 mSLE. During the Holocene collapse of the IIS, the ice divide at the GrIS moves toward the west and the ice streams reorganize in northern Greenland, as shown in Fig. 6. This divide migration could explain the onset of NEGIS, as found by Franke et al. (2022), and the shutdown of the older, more northern ice stream, as observed by Jansen et al. (2024). Figure 9 shows the rate of change of grounded ice for the ensemble member with the highest combined likelihood for the seven basins of the GrIS. The GrIS rate of change becomes negative at −10.7ka and exhibits two distinct peaks: one at −7.8 ka, with a mass-loss rate of 548 Gta−1, and another at −4.95 ka, following the onset of sub-shelf melting, with a mass-loss rate of 511Gta−1. It continues to be negative for the rest of the Holocene, except for a few times during the last 2 ka, where the average mass-loss rate is 23.7Gta−1. The mass-loss rates stated here are averaged over 50 years. The Cryosphere, 19, 3599–3622, 2025 https://doi.org/10.5194/tc-19-3599-2025 M. L. Lauritzen et al.: Modeled Greenland Ice Sheet evolution constrained by Holocene elevation histories 3607 Figure 5. Kernel density estimates of the inferred marginal PDFs for the 20 model parameters that we varied. The units on the yaxes are the inverse of those on the xaxes. Figure 6. Time slices showing modeled surface speed, streamlines, bed topography, and ice shelf extent for the ensemble member with the highest combined likelihood at −12, −9ka, and the present day (PD). The present-day locations of the ice-core sites Camp Century (CC), NGRIP (NG), GRIP (GR), and Dye 3 (D3) are shown. https://doi.org/10.5194/tc-19-3599-2025 The Cryosphere, 19, 3599–3622, 2025 3614 M. L. Lauritzen et al.: Modeled Greenland Ice Sheet evolution constrained by Holocene elevation histories Appendix A Figure A1. Paleoclimatic temperature anomalies derived from 18O measurements at GRIP and NGRIP using linear transfer function (Huybrechts, 2002) and quadratic transfer function from Johnsen et al. (1995) and 18O measurements at Renland and Agassiz (Vinther et al., 2009) and 18O measurements at NGRIP using an inversion scheme (Gkinis et al., 2014). The temperature anomaly from Lecavalier (2017) is also shown as a reference. Figure A2. (a) Annual mean precipitation and (b) summer (June, July, and August) mean 2 m temperatures for our 30-year reference climatology (1960–1989). The Cryosphere, 19, 3599–3622, 2025 https://doi.org/10.5194/tc-19-3599-2025 M. L. Lauritzen et al.: Modeled Greenland Ice Sheet evolution constrained by Holocene elevation histories 3615 Figure A3. Overview of different domain boundaries used to patch together the 30-year reference climatology (ZGRN11, FGRN055, SCAA, NCAA) and the bedrock topography (BedMachine, IBCAO). Figure A4. Scatter plots of the last 500 years of mass-loss rates vs each parameter varied in our ensemble; ρis the Pearson correlation. https://doi.org/10.5194/tc-19-3599-2025 The Cryosphere, 19, 3599–3622, 2025 3616 M. L. Lauritzen et al.: Modeled Greenland Ice Sheet evolution constrained by Holocene elevation histories Figure A5. Modeled surface elevation at Camp Century, color coded for each parameter. The Cryosphere, 19, 3599–3622, 2025 https://doi.org/10.5194/tc-19-3599-2025 M. L. Lauritzen et al.: Modeled Greenland Ice Sheet evolution constrained by Holocene elevation histories 3617 Figure A6. Observed and modeled surface elevation over the past 11.7ka at the ice-core sites Camp Century, NGRIP, GRIP, and Dye 3. The blue lines represent ice-core-derived surface elevations from Vinther et al. (2009), with the blue envelopes indicating 1 standard deviation. The orange, red, green, and black lines correspond to the modeled surface elevations for the 0th, 1st, 2nd, and 19th iterations of the bedrock adjustment, respectively. Figure A7. Modeled bedrock elevation and thinning over the past 11.7 ka at the ice-core sites Camp Century, NGRIP, GRIP, and Dye 3. The red envelopes represent the ensemble-estimated mean and standard deviation, while the green envelopes show the site-specific estimate. The black dots indicate the observed present-day bedrock elevation from Morlighem et al. (2017). https://doi.org/10.5194/tc-19-3599-2025 The Cryosphere, 19, 3599–3622, 2025 3618 M. L. Lauritzen et al.: Modeled Greenland Ice Sheet evolution constrained by Holocene elevation histories Figure A8. Model state at the branch-off point at −20 ka. (a) Modeled surface velocity with streamlines and ice shelf extent. The present-day locations of the ice-core sites Camp Century (CC), NGRIP (NG), GRIP (GR), and Dye 3 (D3) are overlaid. (b) Bedrock topography at −20 ka relative to sea level. (c) Difference in bedrock topography at −20ka compared with the present-day observed topography (bPD −b20ka). (d) Modeled bedrock uplift rates at −20 ka. Code and data availability. PISM is open-source software that can be downloaded from https://github.com/pism/pism (last access: 4 September 2025; https://doi.org/10.5281/zenodo.10202029, Khrulev et al., 2023) (Bueler and Brown, 2009; Winkelmann et al., 2011). Surface-elevation data from the four ice-core locations are available upon request. The presented RACMO data are available upon request and without conditions from Brice Noël ([email protected]). Oxygen isotope records from GRIP and NGRIP are accessible at https://doi.org/10.1594/PANGAEA.55091 (Johnsen, 1999) and https://doi.org/10.1594/PANGAEA.586886 (North Greenland Ice Core Project Members, 2007), respectively. Model output and glacier catchment basin data are available at https://doi.org/10.5281/zenodo.15681862 (Lauritzen, 2025). Video supplement. A video showing GrIS evolution through the Holocene for the most likely ensemble member can be found at https://doi.org/10.5446/68337 (Lauritzen et al., 2024). Author contributions. MLL, CSH, and AS designed the study. MLL prepared the data, performed the model runs, and carried out the subsequent analysis. All authors discussed and improved the paper. Competing interests. At least one of the (co-)authors is a member of the editorial board of The Cryosphere. The peer-review process was guided by an independent editor, and the authors also have no other competing interests to declare. Disclaimer. Views and opinions expressed are those of the authors only and do not necessarily reflect the views or opinions of the European Union or the European Research Council Executive Agency. Neither the European Union nor the granting authority can be held responsible for them. Publisher’s note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors. Acknowledgements. We would like to thank Benoit Lecavalier for providing the temperature anomalies derived from the Agassiz ice core in Lecavalier et al. (2017). We would also like to thank the two anonymous reviewers and the editor, Alexander Robinson, for their help in improving the manuscript. Financial support. Mikkel Langgaard Lauritzen was funded by the Independent Research Fund Denmark through the project GreenPlanning (grant no. 0217-00244B). Christine Schøtt Hvidberg, Nicholas Mossor Rathmann, and Aslak Grindsted received funding from the Novo Nordisk Foundation (grant no. NNF23OC0081251), the Independent Research Fund Denmark (DFF) (grant no. 203200364B), and the Villum Foundation (grant no. 23261). Anne Solgaard was funded by the European Union (ERC, Green2Ice, 101072180). Brice Noel was funded by Fonds de la Recherche Scientifique de Belgique – F.R.S. (FNRS) (grant no. 34805166). Review statement. This paper was edited by Alexander Robinson and reviewed by two anonymous referees. The Cryosphere, 19, 3599–3622, 2025 https://doi.org/10.5194/tc-19-3599-2025 M. L. Lauritzen et al.: Modeled Greenland Ice Sheet evolution constrained by Holocene elevation histories 3619 References A¯ dalgeirsdóttir, G., Aschwanden, A., Khroulev, C., Boberg, F., Mottram, R., Lucas-Picher, P., and Christensen, J.: Role of Model Initialization for Projections of 21st-Century Greenland Ice Sheet Mass Loss, J. Glaciol., 60, 782–794, https://doi.org/10.3189/2014JoG13J202, 2014. Albrecht, T., Winkelmann, R., and Levermann, A.: Glacial-cycle simulations of the Antarctic Ice Sheet with the Parallel Ice Sheet Model (PISM) – Part 2: Parameter ensemble analysis, The Cryosphere, 14, 633–656, https://doi.org/10.5194/tc-14633-2020, 2020. Andersen, K. K., Azuma, N., Barnola, J.-M., Bigler, M., Biscaye, P., Caillon, N., Chappellaz, J., Clausen, H. B., Dahl-Jensen, D., Fischer, H., Flückiger, J., Fritzsche, D., Fujii, Y., Goto-Azuma, K., Grønvold, K., Gundestrup, N. S., Hansson, M., Huber, C., Hvidberg, C. S., Johnsen, S. J., Jonsell, U., Jouzel, J., Kipfstuhl, S., Landais, A., Leuenberger, M., Lorrain, R., Masson-Delmotte, V., Miller, H., Motoyama, H., Narita, H., Popp, T., Rasmussen, S. O., Raynaud, D., Rothlisberger, R., Ruth, U., Samyn, D., Schwander, J., Shoji, H., Siggard-Andersen, M.-L., Steffensen, J. P., Stocker, T., Sveinbjörnsdóttir, A. E., Svensson, A., Takata, M., Tison, J.-L., Thorsteinsson, Th., Watanabe, O., Wilhelms, F., White, J. W. C., and North Greenland Ice Core Project members: High-Resolution Record of Northern Hemisphere Climate Extending into the Last Interglacial Period, Nature, 431, 147–151, https://doi.org/10.1038/nature02805, 2004. Aschwanden, A. and Brinkerhoff, D. J.: Calibrated Mass Loss Predictions for the Greenland Ice Sheet, Geophys. Res. Lett., 49, e2022GL099058, https://doi.org/10.1029/2022GL099058, 2022. Aschwanden, A., Fahnestock, M. A., and Truffer, M.: Complex Greenland Outlet Glacier Flow Captured, Nat. Commun., 7, 10524, https://doi.org/10.1038/ncomms10524, 2016. Aschwanden, A., Fahnestock, M. A., Truffer, M., Brinkerhoff, D. J., Hock, R., Khroulev, C., Mottram, R., and Khan, S. A.: Contribution of the Greenland Ice Sheet to Sea Level over the next Millennium, Science Advances, 5, eaav9396, https://doi.org/10.1126/sciadv.aav9396, 2019. Aschwanden, A., Bartholomaus, T. C., Brinkerhoff, D. J., and Truffer, M.: Brief communication: A roadmap towards credible projections of ice sheet contribution to sea level, The Cryosphere, 15, 5705–5715, https://doi.org/10.5194/tc-15-5705-2021, 2021. Badgeley, J. A., Steig, E. J., Hakim, G. J., and Fudge, T. J.: Greenland temperature and precipitation over the last 20 000 years using data assimilation, Clim. Past, 16, 1325–1346, https://doi.org/10.5194/cp-16-1325-2020, 2020. Bagherbandi, M., Amin, H., Wang, L., and Shirazian, M.: Mantle Viscosity Derived From Geoid and Different Land Uplift Data in Greenland, J. Geophys. Res.-Sol. Ea., 127, e2021JB023351, https://doi.org/10.1029/2021JB023351, 2022. Bamber, J. L., Griggs, J. A., Hurkmans, R. T. W. L., Dowdeswell, J. A., Gogineni, S. P., Howat, I., Mouginot, J., Paden, J., Palmer, S., Rignot, E., and Steinhage, D.: A new bed elevation dataset for Greenland, The Cryosphere, 7, 499–510, https://doi.org/10.5194/tc-7-499-2013, 2013. Bindschadler, R. A., Nowicki, S., Abe-Ouchi, A., Aschwanden, A., Choi, H., Fastook, J., Granzow, G., Greve, R., Gutowski, G., Herzfeld, U., Jackson, C., Johnson, J., Khroulev, C., Levermann, A., Lipscomb, W. H., Martin, M. A., Morlighem, M., Parizek, B. R., Pollard, D., Price, S. F., Ren, D., Saito, F., Sato, T., Seddik, H., Seroussi, H., Takahashi, K., Walker, R., and Wang, W. L.: Ice-Sheet Model Sensitivities to Environmental Forcing and Their Use in Projecting Future Sea Level (the SeaRISE Project), J. Glaciol., 59, 195–224, https://doi.org/10.3189/2013JoG12J125, 2013. Braithwaite, R. J.: Calculation of Degree-Days for Glacier-Climate Research, Zeitschrift für Gletscherkunde und Glazialgeologie, 20/1984, 1–8, 1985. Briner, J. P., Cuzzone, J. K., Badgeley, J. A., Young, N. E., Steig, E. J., Morlighem, M., Schlegel, N.-J., Hakim, G. J., Schaefer, J. M., Johnson, J. V., Lesnek, A. J., Thomas, E. K., Allan, E., Bennike, O., Cluett, A. A., Csatho, B., De Vernal, A., Downs, J., Larour, E., and Nowicki, S.: Rate of Mass Loss from the Greenland Ice Sheet Will Exceed Holocene Values This Century, Nature, 586, 70–74, https://doi.org/10.1038/s41586-020-2742-6, 2020. Bueler, E. and Brown, J.: Shallow Shelf Approximation as a “Sliding Law” in a Thermomechanically Coupled Ice Sheet Model, J. Geophys. Res., 114, F03008, https://doi.org/10.1029/2008JF001179, 2009. Bueler, E. and van Pelt, W.: Mass-conserving subglacial hydrology in the Parallel Ice Sheet Model version 0.6, Geosci. Model Dev., 8, 1613–1635, https://doi.org/10.5194/gmd-8-1613-2015, 2015. Bueler, E., Lingle, C. S., and Brown, J.: Fast Computation of a Viscoelastic Deformable Earth Model for Ice-Sheet Simulations, Ann. Glaciol., 46, 97–105, https://doi.org/10.3189/172756407782871567, 2007. Clark, P. U., He, F., Golledge, N. R., Mitrovica, J. X., Dutton, A., Hoffman, J. S., and Dendy, S.: Oceanic Forcing of Penultimate Deglacial and Last Interglacial Sea-Level Rise, Nature, 577, 660–664, https://doi.org/10.1038/s41586-020-1931-7, 2020. Couette, P.-O., Lajeunesse, P., Ghienne, J.-F., Dorschel, B., Gebhardt, C., Hebbeln, D., and Brouard, E.: Evidence for an Extensive Ice Shelf in Northern Baffin Bay during the Last Glacial Maximum, Commun. Earth Environ., 3, 225, https://doi.org/10.1038/s43247-022-00559-7, 2022. Dansgaard, W., Johnsen, S. J., Møller, J., and Langway, C. C.: One Thousand Centuries of Climatic Record from Camp Century on the Greenland Ice Sheet, Science, 166, 377–381, https://doi.org/10.1126/science.166.3903.377, 1969. England, J., Atkinson, N., Bednarski, J., Dyke, A., Hodgson, D., and Ó Cofaigh, C.: The Innuitian Ice Sheet: Configuration, Dynamics and Chronology, Quaternary Sci. Rev., 25, 689–703, https://doi.org/10.1016/j.quascirev.2005.08.007, 2006. Eyring, V., Gillett, N., Achuta Rao, K., Barimalala, R., Barreiro Parrillo, M., Bellouin, N., Cassou, C., Durack, P., Kosaka, Y., McGregor, S., Min, S., Morgenstern, O., and Sun, Y.: Human Influence on the Climate System, in: Climate Change 2021: The Physical Science Basis. Contribution of Working Group I to the Sixth Assessment Report of the Intergovernmental Panel on Climate Change, edited by: Masson-Delmotte, V., Zhai, P., Pirani, A., Connors, S., Péan, C., Berger, S., Caud, N., Chen, Y., Goldfarb, L., Gomis, M., Huang, M., Leitzell, K., Lonnoy, E., Matthews, J., Maycock, T., Waterfield, T., Yelekçi, O., Yu, R., and Zhou, B., Cambridge University Press, Cambridge, United Kingdom and New York, NY, USA, 423–552, https://doi.org/10.1017/9781009157896.005, 2021. https://doi.org/10.5194/tc-19-3599-2025 The Cryosphere, 19, 3599–3622, 2025 3620 M. L. Lauritzen et al.: Modeled Greenland Ice Sheet evolution constrained by Holocene elevation histories Flanders Marine Institute (VLIZ), Belgium: Maritime Boundaries Geodatabase: Maritime Boundaries and Exclusive Economic Zones (200NM), Version 12, Flanders Marine Institute (VLIZ) [data set], https://doi.org/10.14284/632, 2023. Franke, S., Bons, P. D., Westhoff, J., Weikusat, I., Binder, T., Streng, K., Steinhage, D., Helm, V., Eisen, O., Paden, J. D., Eagles, G., and Jansen, D.: Holocene Ice-Stream Shutdown and Drainage Basin Reconfiguration in Northeast Greenland, Nat. Geosci., 15, 995–1001, https://doi.org/10.1038/s41561-022-01082-2, 2022. GEBCO Bathymetric Compilation Group: The GEBCO_2023 Grid – a Continuous Terrain Model of the Global Oceans and Land, British Oceanographic Data Centre [data set], https://doi.org/10.5285/F98B053B-0CBC-6C23-E0536C86ABC0AF7B, 2023. Gkinis, V., Simonsen, S., Buchardt, S., White, J., and Vinther, B.: Water Isotope Diffusion Rates from the NorthGRIP Ice Core for the Last 16,000 Years – Glaciological and Paleoclimatic Implications, Earth Planet Sc. Lett., 405, 132–141, https://doi.org/10.1016/j.epsl.2014.08.022, 2014. Gladstone, R. M., Payne, A. J., and Cornford, S. L.: Parameterising the grounding line in flow-line ice sheet models, The Cryosphere, 4, 605–619, https://doi.org/10.5194/tc-4-605-2010, 2010. Huybrechts, P.: Sea-Level Changes at the LGM from Ice-Dynamic Reconstructions of the Greenland and Antarctic Ice Sheets during the Glacial Cycles, Quaternary Sci. Rev., 21, 203–231, https://doi.org/10.1016/S0277-3791(01)00082-8, 2002. Imbrie, J. D. and McIntyre, A.: SPECMAP Time Scale Developed by Imbrie et al., 1984 Based on Normalized Planktonic Records (Normalized O-18 vs Time, Specmap.017), PANGAEA [data set], https://doi.org/10.1594/PANGAEA.441706, 2006. Jakobsson, M., Mayer, L. A., Bringensparr, C., Castro, C. F., Mohammad, R., Johnson, P., Ketter, T., Accettella, D., Amblas, D., An, L., Arndt, J. E., Canals, M., Casamor, J. L., Chauché, N., Coakley, B., Danielson, S., Demarte, M., Dickson, M.-L., Dorschel, B., Dowdeswell, J. A., Dreutter, S., Fremand, A. C., Gallant, D., Hall, J. K., Hehemann, L., Hodnesdal, H., Hong, J., Ivaldi, R., Kane, E., Klaucke, I., Krawczyk, D. W., Kristoffersen, Y., Kuipers, B. R., Millan, R., Masetti, G., Morlighem, M., Noormets, R., Prescott, M. M., Rebesco, M., Rignot, E., Semiletov, I., Tate, A. J., Travaglini, P., Velicogna, I., Weatherall, P., Weinrebe, W., Willis, J. K., Wood, M., Zarayskaya, Y., Zhang, T., Zimmermann, M., and Zinglersen, K. B.: The International Bathymetric Chart of the Arctic Ocean Version 4.0, Scientific Data, 7, 176, https://doi.org/10.1038/s41597-020-0520-9, 2020. Jansen, D., Franke, S., Bauer, C. C., Binder, T., Dahl-Jensen, D., Eichler, J., Eisen, O., Hu, Y., Kerch, J., Llorens, M.-G., Miller, H., Neckel, N., Paden, J., Riese, T., Sachau, T., Stoll, N., Weikusat, I., Wilhelms, F., Zhang, Y., and Bons, P. D.: Shear Margins in Upper Half of Northeast Greenland Ice Stream Were Established Two Millennia Ago, Nat. Commun., 15, 1193, https://doi.org/10.1038/s41467-024-45021-8, 2024. Johnsen, S. J., Clausen, H. B., Dansgaard, W., Gundestrup, N. S., Hansson, M., Jonsson, P., Steffensen, J. P., and Sveinbjørnsdottir, A. E.: A “Deep” Ice Core from East Greenland, Meddelelser om Grønland, Geoscience, 29, 1–22, https://doi.org/10.7146/moggeosci.v29i.140329, 1992. Johnsen, S. J., Dahl-Jensen, D., Dansgaard, W., and Gundestrup, N.: Greenland Palaeotemperatures Derived from GRIP Bore Hole Temperature and Ice Core Isotope Profiles, Tellus B, 47, 624– 629, https://doi.org/10.3402/tellusb.v47i5.16077, 1995. Johnsen, S. J.: GRIP Oxygen Isotopes, PANGAEA [data set], https://doi.org/10.1594/PANGAEA.55091, 1999. Khrulev, C., Aschwanden, A., Bueler, E., Brown, J., Maxwell, D., Albrecht, T., Reese, R., Mengel, M., Martin, M., Winkelmann, R., Zeitz, M., Levermann, A., Feldmann, J., Garbe, J., Haseloff, M., Seguinot, J., Hinck, S., Kleiner, T., Fischer, E., Damsgaard, A., Lingle, C., van Pelt, W., Ziemen, F., Shemonski, N., Mankoff, K., Kennedy, J., Blum, K., Habermann, M., DellaGiustina, D., Hock, R., Kreuzer, M., Degregori, E., Schoell, S.: Parallel Ice Sheet Model (PISM) (v2.1), Zenodo [code], https://doi.org/10.5281/zenodo.10202029, 2023. Lambeck, K., Rouby, H., Purcell, A., Sun, Y., and Sambridge, M.: Sea Level and Global Ice Volumes from the Last Glacial Maximum to the Holocene, P. Natl. Acad. Sci. USA, 111, 15296– 15303, https://doi.org/10.1073/pnas.1411762111, 2014. Lauritzen, M.: Dataset belonging to the article; Modeled Greenland Ice Sheet evolution constrained by ice-core-derived Holocene elevation histories [Data set], Zenodo [data set], https://doi.org/10.5281/zenodo.15681862, 2025. Lauritzen, M., Aðalgeirsdóttir, G., Rathmann, N., Grinsted, A., Noël, B., and Hvidberg, C. S.: The Influence of Inter-Annual Temperature Variability on the Greenland Ice Sheet Volume, Ann. Glaciol., 64, 1–8, https://doi.org/10.1017/aog.2023.53, 2023. Lauritzen, M. L., Solgaard, A., Rathmann, N., Vinther, B. M., Grindsted, A., Noël, B., Aðalgeirsdóttir, G., and Hvidberg, C. S.: Modeled Greenland Ice Sheet evolution constrained by ice-corederived Holocene elevation histories, TIB AV-Portal [video], https://doi.org/10.5446/68337, 2024. Lecavalier, B. S., Milne, G. A., Vinther, B. M., Fisher, D. A., Dyke, A. S., and Simpson, M. J.: Revised Estimates of Greenland Ice Sheet Thinning Histories Based on Ice-Core Records, Quaternary Sci. Rev., 63, 73–82, https://doi.org/10.1016/j.quascirev.2012.11.030, 2013. Lecavalier, B. S., Fisher, D. A., Milne, G. A., Vinther, B. M., Tarasov, L., Huybrechts, P., Lacelle, D., Main, B., Zheng, J., Bourgeois, J., and Dyke, A. S.: High Arctic Holocene Temperature Record from the Agassiz Ice Cap and Greenland Ice Sheet Evolution, P. Natl. Acad. Sci. USA, 114, 5952–5957, https://doi.org/10.1073/pnas.1616287114, 2017. Leger, T. P. M., Clark, C. D., Huynh, C., Jones, S., Ely, J. C., Bradley, S. L., Diemont, C., and Hughes, A. L. C.: A Greenlandwide empirical reconstruction of paleo ice sheet retreat informed by ice extent markers: PaleoGrIS version 1.0, Clim. Past, 20, 701–755, https://doi.org/10.5194/cp-20-701-2024, 2024. Lingle, C. S. and Clark, J. A.: A Numerical Model of Interactions between a Marine Ice Sheet and the Solid Earth: Application to a West Antarctic Ice Stream, J. Geophys. Res.-Oceans, 90, 1100– 1114, https://doi.org/10.1029/JC090iC01p01100, 1985. Liu, Z., Otto-Bliesner, B. L., He, F., Brady, E. C., Tomas, R., Clark, P. U., Carlson, A. E., Lynch-Stieglitz, J., Curry, W., Brook, E., Erickson, D., Jacob, R., Kutzbach, J., and Cheng, J.: Transient Simulation of Last Deglaciation with a New Mechanism for Bølling-Allerød Warming, Science, 325, 310–314, https://doi.org/10.1126/science.1171041, 2009. The Cryosphere, 19, 3599–3622, 2025 https://doi.org/10.5194/tc-19-3599-2025 M. L. Lauritzen et al.: Modeled Greenland Ice Sheet evolution constrained by Holocene elevation histories 3621 MacGregor, J. A., Colgan, W. T., Fahnestock, M. A., Morlighem, M., Catania, G. A., Paden, J. D., and Gogineni, S. P.: Holocene Deceleration of the Greenland Ice Sheet, Science, 351, 590–593, https://doi.org/10.1126/science.aab1702, 2016. Millan, R., Mouginot, J., Rabatel, A., and Morlighem, M.: Ice Velocity and Thickness of the World’s Glaciers, Nat. Geosci., 15, 124–129, https://doi.org/10.1038/s41561-021-00885-z, 2022. Morlighem, M.: IceBridge BedMachine Greenland, Version 5, National Snow and Ice Data Center Distributed Active Archive Center [data set], https://doi.org/10.5067/GMEVBWFLWA7X, 2022. Morlighem, M., Bondzio, J., Seroussi, H., Rignot, E., Larour, E., Humbert, A., and Rebuffi, S.: Modeling of Store Gletscher’s Calving Dynamics, West Greenland, in Response to Ocean Thermal Forcing, Geophys. Res. Lett., 43, 2659–2666, https://doi.org/10.1002/2016GL067695, 2016. Morlighem, M., Williams, C. N., Rignot, E., An, L., Arndt, J. E., Bamber, J. L., Catania, G., Chauché, N., Dowdeswell, J. A., Dorschel, B., Fenty, I., Hogan, K., Howat, I., Hubbard, A., Jakobsson, M., Jordan, T. M., Kjeldsen, K. K., Millan, R., Mayer, L., Mouginot, J., Noël, B. P. Y., O’Cofaigh, C., Palmer, S., Rysgaard, S., Seroussi, H., Siegert, M. J., Slabon, P., Straneo, F., van den Broeke, M. R., Weinrebe, W., Wood, M., and Zinglersen, K. B.: BedMachine v3: Complete Bed Topography and Ocean Bathymetry Mapping of Greenland From Multibeam Echo Sounding Combined With Mass Conservation, Geophys. Res. Lett., 44, 11051–11061, https://doi.org/10.1002/2017GL074954, 2017. Mouginot, J. and Rignot, E.: Glacier Catchments/Basins for the Greenland Ice Sheet, Dryad [data set], https://doi.org/10.7280/D1WT11, 2019. Nias, I. J., Nowicki, S., Felikson, D., and Loomis, B.: Modeling the Greenland Ice Sheet’s Committed Contribution to Sea Level During the 21st Century, J. Geophys. Res.-Earth, 128, e2022JF006914, https://doi.org/10.1029/2022JF006914, 2023. Nielsen, L. T., Aðalgeirsdóttir, Gu., Gkinis, V., Nuterman, R., and Hvidberg, C. S.: The Effect of a Holocene Climatic Optimum on the Evolution of the Greenland Ice Sheet during the Last 10Kyr, J. Glaciol., 64, 477–488, https://doi.org/10.1017/jog.2018.40, 2018. Noël, B., van de Berg, W. J., van Meijgaard, E., Kuipers Munneke, P., van de Wal, R. S. W., and van den Broeke, M. R.: Evaluation of the updated regional climate model RACMO2.3: summer snowfall impact on the Greenland Ice Sheet, The Cryosphere, 9, 1831–1844, https://doi.org/10.5194/tc-9-1831-2015, 2015. Noël, B., van de Berg, W. J., Lhermitte, S., Wouters, B., Schaffer, N., and van den Broeke, M. R.: Six Decades of Glacial Mass Loss in the Canadian Arctic Archipelago, J. Geophys. Res.-Earth, 123, 1430–1449, https://doi.org/10.1029/2017JF004304, 2018. Noël, B., Van De Berg, W. J., Lhermitte, S., and Van Den Broeke, M. R.: Rapid Ablation Zone Expansion Amplifies North Greenland Mass Loss, Science Advances, 5, eaaw0123, https://doi.org/10.1126/sciadv.aaw0123, 2019. North Greenland Ice Core Project Members: 50 year means of oxygen isotope data from ice core NGRIP, PANGAEA [data set], https://doi.org/10.1594/PANGAEA.586886, 2007. Nowicki, S., Goelzer, H., Seroussi, H., Payne, A. J., Lipscomb, W. H., Abe-Ouchi, A., Agosta, C., Alexander, P., Asay-Davis, X. S., Barthel, A., Bracegirdle, T. J., Cullather, R., Felikson, D., Fettweis, X., Gregory, J. M., Hattermann, T., Jourdain, N. C., Kuipers Munneke, P., Larour, E., Little, C. M., Morlighem, M., Nias, I., Shepherd, A., Simon, E., Slater, D., Smith, R. S., Straneo, F., Trusel, L. D., van den Broeke, M. R., and van de Wal, R.: Experimental protocol for sea level projections from ISMIP6 stand-alone ice sheet models, The Cryosphere, 14, 2331–2368, https://doi.org/10.5194/tc-14-2331-2020, 2020. Peltier, W. R.: GLOBAL GLACIAL ISOSTASY AND THE SURFACE OF THE ICE-AGE EARTH: The ICE-5G (VM2) Model and GRACE, Annu. Rev. Earth Pl. Sc., 32, 111–149, https://doi.org/10.1146/annurev.earth.32.082503.144359, 2004. Rasmussen, S. O., Andersen, K. K., Svensson, A. M., Steffensen, J. P., Vinther, B. M., Clausen, H. B., Siggaard-Andersen, M.-L., Johnsen, S. J., Larsen, L. B., Dahl-Jensen, D., Bigler, M., Röthlisberger, R., Fischer, H., Goto-Azuma, K., Hansson, M. E., and Ruth, U.: A New Greenland Ice Core Chronology for the Last Glacial Termination, J. Geophys. Res.-Atmos., 111, D06102, https://doi.org/10.1029/2005JD006079, 2006. RGI Consortium: Randolph Glacier Inventory – A Dataset of Global Glacier Outlines, Version 7, National Snow and Ice Data Center Distributed Active Archive Center [data set], https://doi.org/10.5067/F6JMOVY5NAVZ, 2023. Robinson, A. and Goelzer, H.: The importance of insolation changes for paleo ice sheet modeling, The Cryosphere, 8, 1419– 1428, https://doi.org/10.5194/tc-8-1419-2014, 2014. Schumacher, M., King, M. A., Rougier, J., Sha, Z., Khan, S. A., and Bamber, J. L.: A New Global GPS Data Set for Testing and Improving Modelled GIA Uplift Rates, Geophys. J.l Int., 214, 2164–2176, https://doi.org/10.1093/gji/ggy235, 2018. Shapiro, N.: Inferring Surface Heat Flux Distributions Guided by a Global Seismic Model: Particular Application to Antarctica, Earth Planet Sc. Lett., 223, 213–224, https://doi.org/10.1016/j.epsl.2004.04.011, 2004. Solgaard, A., Kusk, A., Merryman Boncori, J. P., Dall, J., Mankoff, K. D., Ahlstrøm, A. P., Andersen, S. B., Citterio, M., Karlsson, N. B., Kjeldsen, K. K., Korsgaard, N. J., Larsen, S. H., and Fausto, R. S.: Greenland ice velocity maps from the PROMICE project, Earth Syst. Sci. Data, 13, 3491–3512, https://doi.org/10.5194/essd-13-3491-2021, 2021. Tabone, I., Robinson, A., Montoya, M., and Alvarez-Solas, J.: Holocene Thinning in Central Greenland Controlled by the Northeast Greenland Ice Stream, Nat. Commun., 15, 6434, https://doi.org/10.1038/s41467-024-50772-5, 2024. Tang, B.: Orthogonal Array-Based Latin Hypercubes, J. Am. Stat. Assoc., 88, 1392–1397, https://doi.org/10.2307/2291282, 1993. The IMBIE Team: Mass Balance of the Greenland Ice Sheet from 1992 to 2018, Nature, 579, 233–239, https://doi.org/10.1038/s41586-019-1855-2, 2020. van Calcar, C. J., van de Wal, R. S. W., Blank, B., de Boer, B., and van der Wal, W.: Simulation of a fully coupled 3D glacial isostatic adjustment – ice sheet model for the Antarctic ice sheet over a glacial cycle, Geosci. Model Dev., 16, 5473–5492, https://doi.org/10.5194/gmd-16-5473-2023, 2023. Vinther, B. M., Buchardt, S. L., Clausen, H. B., Dahl-Jensen, D., Johnsen, S. J., Fisher, D. A., Koerner, R. M., Raynaud, D., Lipenkov, V., Andersen, K. K., Blunier, T., Rasmussen, S. O., Steffensen, J. P., and Svensson, A. M.: Holocene Thinning of the Greenland Ice Sheet, Nature, 461, 385–388, https://doi.org/10.1038/nature08355, 2009. https://doi.org/10.5194/tc-19-3599-2025 The Cryosphere, 19, 3599–3622, 2025 3622 M. L. Lauritzen et al.: Modeled Greenland Ice Sheet evolution constrained by Holocene elevation histories Winkelmann, R., Martin, M. A., Haseloff, M., Albrecht, T., Bueler, E., Khroulev, C., and Levermann, A.: The Potsdam Parallel Ice Sheet Model (PISM-PIK) – Part 1: Model description, The Cryosphere, 5, 715–726, https://doi.org/10.5194/tc-5-715-2011, 2011. Yokoyama, Y., Esat, T. M., Thompson, W. G., Thomas, A. L., Webster, J. M., Miyairi, Y., Sawada, C., Aze, T., Matsuzaki, H., Okuno, J., Fallon, S., Braga, J.-C., Humblet, M., Iryu, Y., Potts, D. C., Fujita, K., Suzuki, A., and Kan, H.: Rapid Glaciation and a Two-Step Sea Level Plunge into the Last Glacial Maximum, Nature, 559, 603–607, https://doi.org/10.1038/s41586-018-0335-4, 2018. The Cryosphere, 19, 3599–3622, 2025 https://doi.org/10.5194/tc-19-3599-2025