the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
| 25 Jun 2026
Ensemble reconstruction of the Greenland Ice Sheet evolution through the last deglaciation
Abstract. The last deglaciation offers valuable insights into ice-climate interactions, as extensive paleoclimatic records document the retreat of ice sheets through a period of major climate changes. During this interval, the Greenland Ice Sheet (GrIS) retreated from its extensive Last Glacial Maximum (LGM) configuration to its present state, passing through the Holocene Thermal Maximum (HTM), when temperatures exceeded present-day values. Despite the large amount of paleoclimatic data available, ice-sheet models struggle to reproduce key aspects of the observational record, and the magnitude of the GrIS contribution to sea level throughout this period, in particular during the LGM and the HTM, remains highly uncertain. In this study, we evaluate an ensemble of 3,000 simulations of the GrIS performed with the Yelmo ice-sheet model against different observational constraints. These include: (1) the LGM ice-sheet extent, (2) ice-core-derived surface elevations, (3) an ice-extent retreat chronology based on the recent PaleoGrIS dataset, and (4) the present-day ice-sheet configuration (ice thickness, ice cover, ice-surface velocity, and bedrock elevation). We characterize the impact of the parameters perturbed along the ensemble on the GrIS evolution using an emulator based on the XGBoost algorithm combined with the SHAP (SHapley Additive exPlanations) framework. This analysis reveals that the climatic parameters (surface melting corrections and ocean thermal sensitivity) dominate the impact. By identifying the simulation that best matches these observables, we provide a constrained reconstruction of the GrIS during the last deglaciation that substantially improves upon previous reconstructions. We obtain a GrIS contribution to global sea level with respect to present of -5.75 m of sea-level equivalent (SLE) at the LGM (with an uncertainty range of -3.93 to -6.31 m) and +0.45 m after the HTM warming (with an uncertainty range of 0.43 to 1.18 m). This is in the upper range of previously existing estimates, indicating a comparatively mid-to-high GrIS sensitivity to climate changes over the last deglaciation.
Competing interests: At least one of the (co-)authors is a member of the editorial board of Climate of the Past.
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. Views expressed in the text are those of the authors and do not necessarily reflect the views of the publisher.- Preprint
(17175 KB) - Metadata XML
-
Supplement
(9577 KB) - BibTeX
- EndNote
Status: final response (author comments only)
- RC1: 'Review egusphere-2026-3331', Nicolas Young, 23 Jul 2026
-
RC2: 'Review of Gutierrez-Gonzalez Greenland Model Ensemble since the LGM', Anonymous Referee #2, 01 Oct 2026
Gutierrez-Gonzalez and others report and analyze the output of a large ensemble of ice sheet model runs with Yelmo for Greenland since the LGM. The ensemble output was compared to both modern ice sheet configuration, the PaleoGRIS extent reconstruction, and ice-core based thinning histories. Thirteen different parameters were adjusted in the ensemble and the analysis focused on which parameters controlled the best-fit members of the ensemble. The authors use the reconstruction to assess the contribution of the Greenland ice sheet to sea level rise since the LGM.
The manuscript has a strong layout and proceeds logically. The ensemble itself has been presented before (Tabone et al., 2024), although the prior analysis focused on the Northeast Greenland ice stream and interior thinning, rather than Greenland as a whole. This work is a logical extension of the ensemble and a useful contribution. The method of assessing the ensemble members was interesting and well described.
My largest concern about the manuscript is the climate forcing. The authors find that the atmospheric corrections terms are quite important to ice sheet evolution, but I believe they represent a single value for all locations on the ice sheet. I’m surprised that only one combination of surface temperature and precipitation was used. The ensemble used Buizert et al. (2018) for the temperature and Badgeley et al. (2020) for the accumulation. The Buizert and Badgeley reconstructions have large differences in the spatial pattern of climate change. The lack of diversity in the climate forcing prevents fuller assessment of ice sheet model parameters. I recognize that this is not something easy to remedy after the ensemble is has been run, and don’t suggest that the authors need to do it for this paper; but I do think the authors need to highlight this limitation much more strongly in the methods, discussion, and throughout the results. Although likely unintended, my primary takeaway from this paper is that the climate forcing likely matters more than any of the ice sheet model parameters.
I also have a confession to make – I spilled a drink on my paper copy that I was editing on a flight, destroying my notes and specific comments on the text. I was added as a referee late to the review process and am now overdue myself, so I think it’s more important to get the review back and keep the process rolling, than wait for me to have additional time to recreate my notes. I am sorry that this review is light on specific comments.
General Comments
Section 3.3.1 on the Emulator Sensitivity doesn’t seem to add much to the paper. It’s short, confusing, relies on “black box” functions, and seems largely detached from the rest of the paper. I don’t think it is described at all in the methods sections, and it also doesn’t seem like there is enough detail to reproduce the methodology. I understand you want to look at more complex multi-variate relationships, but this doesn’t come through. So my vote is to just remove it to help keep what is a pretty long paper better focused and accessible.
Figure 11 seems like the most important figure for the paper, but I find it confusing. First, I don’t understand why there are both positive and negative values. Is this simply indicating positive or anti correlation for the direction a parameter is varied? If so, that doesn’t seem very important to me and is just a result of the decision on how to formulate the parameter – instead, I think the relevant information is the absolute value of the correlation. Second, did I miss the assessment of statistical significance? Third, and here is where I reveal my ignorance and that I’m interpreting the whole thing wrong, but is a simple correlation a useful statistic. Doesn’t this have a bias for parameters that have a linear relationship with the S scores? If I understand, you are essentially comparing the S values for each ensemble member by the parameter value used in that member. But what happens when the optimal parameter value is the middle, and either an increase or decrease in the parameter causes a decrease in the fit? The correlation would be bad, but there would be an important relationship between the score and the parameter that goes undetected. Would it be worth putting plots with the individual relationships in the supplement?
Fp – I find it confusing that this parameter is only tuned in one basin. It’s hard to keep track of this when the info is not in Table 2 and not called out in later figures. I have two concerns: 1) the paper gives a false sense that precipitation scaling isn’t important and 2) there are other important details of other parameters that I didn’t catch and should be highlighted in Table 2.
Weighting of Model Scores – I think there needs to be more discussion of how to combine the different model scores. The present day ones make sense to me overall. It’s the paleo ones that I’m a bit concerned about. Of the five paleo scores, there are ice sheet thinning history. Of these three, two are essentially the same (GRIP and NGRIP) which is noted in the paper, but we know a priori they are going to be incredibly similar given their proximity at the summit of the ice sheet. The other, Dye3, is wildly more uncertain. It’s easy to forget that Dye3 is not at a divide, so the ice has flowed through strong climate gradients, and the LGM ice is deep and highly thinned. I feel like having these three be separate adds too much weight to them. It seems worth testing how different weighting of the scores going into the total score affects the best model. Looking at Figure 9, it seems possible that you could get a big shift in the parameters for the best fitting model.
Figures 2-4: I’m unclear if “best simulation” is referring to the overall best, or the best for that sub-score. Can you clarify that in the captions and text.
Citation: https://doi.org/10.5194/egusphere-2026-3331-RC2
Data sets
Video and dataset of the best simulation Lucía Gutiérrez-González, Ilaria Tabone, Alexander Robinson, Jorge Álvarez-Solas, Jan Swierczek-Jereczek, Britta Grusdt, and Marisa Montoya https://zenodo.org/records/20544171
Model code and software
Yelmo model Alexander Robinson, Jorge Alvarez-Solas, Marisa Montoya, Heiko Goelzer, Ralf Greve, and Catherine Ritz https://github.com/palma-ice/yelmo/
Interactive computing environment
Sripts for the ensemble analysis and figures Lucía Gutiérrez-González https://github.com/luciagutierrezg/GrIS_deglaciation_yelmox
Video supplement
Video and dataset of the best simulation Lucía Gutiérrez-González, Ilaria Tabone, Alexander Robinson, Jorge Álvarez-Solas, Jan Swierczek-Jereczek, Britta Grusdt, and Marisa Montoya https://zenodo.org/records/20544171
Viewed
| HTML | XML | Total | Supplement | BibTeX | EndNote | |
|---|---|---|---|---|---|---|
| 249 | 130 | 29 | 408 | 79 | 24 | 24 |
- HTML: 249
- PDF: 130
- XML: 29
- Total: 408
- Supplement: 79
- BibTeX: 24
- EndNote: 24
Viewed (geographical distribution)
| Country | # | Views | % |
|---|
| Total: | 0 |
| HTML: | 0 |
| PDF: | 0 |
| XML: | 0 |
- 1
This manuscript by Gutiérrez-González presents a large ensemble of simulations of Greenland Ice Sheet deglaciation and regrowth from the Last Glacial Maximum (LGM) through present. From the large ensemble (original simulations are actually from Tabone et al., 2024 which many of the authors here are also on), the authors single out the best-fit simulation against a suite a metrics used to assess to model performance while at the same appropriately highlighting the range of ensemble results and the specific model parameters controlling the full range of results. The authors use their best-fit simulation and exploration of model parameter space to highlight 1) surface melting corrections and ocean thermal sensitivity as having the most dominant control on model output, 2) estimate of the ice sheet’s sea-level equivalent (SLE) at the LGM, and 3) SLE at its minimum Holocene extent. These latter points, arguably, being of particular importance as they are values that are often discussed and debated in the literature.
This is an excellent paper. A very thorough piece of work that was a pleasure to read. Paper is well written and the figures are excellent. Figures do a great job at conveying information while being aesthetically pleasing. And, as a paleo person (I am not a modeler) I found the LGM and mid-Holocene reconstructions particularly interesting and commend the authors for doing a great job and discussing their results within the context of the geological record. Excellent stuff.
This paper is obviously suitable for Climate of the Past and I can’t imagine the authors receiving much resistance in publishing this. Below I have a number of light to moderate comments and suggestions for the authors but, again, this is an excellent manuscript.
---
Lines 25-27: I found these two sentences oddly structured. Perhaps the word “nevertheless” is making this seem like a bigger mystery than it really is? The authors could even think about just starting with the second sentence minus “nevertheless”
Line 49: I don’t think you need to start a new paragraph here.
Line 71-73: Yes, I agree, but this also seems like a slightly too negative way to cast this. There are tradeoffs with every simulation and no simulation hits all the benchmarks, including the author’s own here.
Modelling framework section:
The ensemble combines atmospheric temperature anomalies from Buizert et al. (2018) with precipitation
anomalies from Badgeley et al. (2020). Could the authors elaborate on the rationale for selecting the Buizert temperature reconstruction rather than the Badgeley temperature reconstruction, given that the precipitation forcing is taken from Badgeley? The two temperature reconstructions differ during climatically important intervals, including the Younger Dryas and the early Holocene, which could influence the simulated timing and magnitude of deglaciation. To what extent might the inferred ice-sheet evolution and best-fit simulation depend on the choice of temperature reconstruction? Related, what is the rationale behind the “high-precipitation” scenario from Badgeley vs. the moderate or low scenarios?
The best-performing parameter set is identified using a single prescribed climate forcing (Buizert temperature anomalies and Badgeley precipitation anomalies). Because several perturbed parameters directly modulate the ice sheet response to climate (e.g., melt corrections, precipitation scaling, and ocean sensitivity), to what extent are the inferred optimal parameter values specific to this forcing? If an alternative paleoclimate reconstruction were adopted, would the calibrated parameter values shift substantially? Additionally, because the SHAP analysis identifies parameter importance within the context of the prescribed forcing, it is unclear whether the inferred parameter sensitivities and perhaps even the SHAP importance rankings themselves would be transferable across alternative climate reconstructions. A brief discussion of this limitation would help clarify the robustness and generality of the conclusions. More broadly, I think it’s worth mentioning or highlighting that the authors here have essentially taken fixed climatology and then sampled model parameter space to achieve a best-fit to the observations while there is also the near-opposite approach of sampling climate space to see of that alone can produce a good fit to the observations (i.e., Briner et al., 2020; Cuzzone et al., pre-print/in review; https://doi.org/10.21203/rs.3.rs-9962320/v1)
Figure 1: Did you want the HTM blue box to come all the way down to the x-axis? Double check
Line 252: Slightly confused here. I don’t think the fit to the Leger product was weighted more heavily, but here makes it sound like it is? Not entirely sure, very possible it is weighted more towards a Leger fit (which is fine) and I missed it.
Line 264-265 and Figure 3: I found the relative poorness of fit in SW Greenland extremely interesting as 1) SW Greenland has the most paleo observations in all of Greenland by a comfortable margin (authors might actually mention this) and 2) it’s probably the region of Greenland least affected by dynamical processes (primarily land-based, KNS region excluded, at least relative to other regions) and mostly controlled by SMB (which is partly why the fit in Briner et al., 2020 is good)…and yet the fit here is relatively poor. Authors could highlight this point, possibly expanded on this with a few sentences. Or maybe at the very least point to section below where they speculate that the climate input is too warm in this region?
Table 4: Likely due to my lack of model knowledge…how “reasonable” are these parameter fits? Is there any single value(s) that are pushing the limits of what is typically accepted/used for each of these parameters?
Figure 5: Why not change the timing of panel G to the actual minimum you reconstruct? I know the difference between the 5 ka and the actual minimum at 4.8 ka is negligible, but it seems odd to me to have 5 ka instead of the 4.8 ka that most people are going to care about.
Line 331: Could add the Sbarra reference in here, and probably Jennings et al., 2017 (v. 472; EPSL) and Ó Cofaigh et al., 2013 (v. 41 (2) Geology). All of these pretty strongly support LGM ice at the shelf edge.
Line 353-355: I don’t think this is necessarily the best way to cast this. I think it’s fairly well known that the minimum extent, whenever it is, falls after peak HTM temperatures but…of course it does? The ice sheet is still so large and has so much inertia that the minimum almost must fall after the HTM. And, of course the ice sheet was indeed retreating/losing mass during the HTM. A more appropriate way to think about this might be looking at the ***rate*** of mass loss. Rates are almost certainly higher during the actual HTM.
Line 359: You can just start this paragraph with ‘At 4.8 ka..”
Line 400-403: Does this suggest the best-fit LGM sea-level equivalent is not a robust estimate?
Line 529-530: Could add Cuzzone et al., 2022 (The Cryosphere, v. 16) here. They investigated this very idea in the KNS system and still couldn’t achieve a great fit with the observations.
Line 563-564: See the last 1/3 of my second modeling framework comment above.
Figure 14 (and related text): It seems like the authors’ model is perhaps more sensitive than other set-ups? Their LGM SLE is towards the higher end of the range of estimates and, conversely, their SLE for the Holocene minimum is also towards the high end? Best fit = both “big” LGM ice and “small” mid Holocene ice.
Line 595-597 and again down ~619-622 and ~635: I might couch all of this a bit differently. I think the authors reconstructions are consistent with all other simulations, that is, they tend to show quite a bit of retreat in SW Greenland during the Holocene. But I also think all of the geological benchmarks over the last ~15 years point to rather limited retreat. So if the question is how far exactly did the ice margin retreat, then it’s a bit unknown, but if the question is did it retreat “a lot” (models) or “not that much” (paleo observations) then I’m not sure I’d call it an “open question” (its not that much).
Holocene minimum SLE – if the best-fit simulation likely overshoots SW Greenland inland retreat, the same region thought to be the primary contributor to SLE during this interval (which the authors fully acknowledge), is the 0.45m SLE itself too high then? Should this (and the full range presented by the authors) be considered a maximum value?