This work investigates the evolution of both the Greenland Ice Sheet (GrIS) and the Canadian Arctic Archipelago during the last deglaciation using simulations of the PISM ice-sheet model. It explores how Holocene surface elevation histories, derived from ice-core data across the ice sheet, can be used to validate model results and refine model parameters. A well-known issue in the field is the difficulty of 3D ice sheet models in accurately simulating Holocene thinning curves across Greenland: this work successfully replicates these curves and attributes the previous data-model discrepancies to the limited modeled ice-sheet advance and the lack of a connection between the Innuitian ice sheet (IIS) and the GrIS during the Last Glacial Maximum (LGM).

Given the scarcity of information on the GrIS evolution away from the margin, surface elevation histories can be an additional useful metric for model validation helping to constrain the evolution of the ice-sheet interior and the overall volume. By considering a domain that includes the IIS and by adjusting the bedrock elevation at the LGM, the authors can replicate the Holocene thinning curves with good agreement upon having calibrated 20 key model parameters.

I find this a very interesting work. It shows how the Holocene ice-core derived surface elevation changes can be used to validate the modeled ice-sheet evolution, as long as the response of the GrIS to non-local changes in the ice load and/or ice dynamics are considered. However, given the high uncertainty in deriving such thinning curves from ice cores and the availability of both more reliable paleodata sources, such as exposure dates of moraines, and present day observations, I think one must be careful in using these elevation histories as - almost - the only metric for validating model simulations. I strongly agree with Jessica Badgeley, who commented in the discussion. While I won't reiterate her points, I believe she raised several important considerations. Additionally, there is a significant discrepancy between the modeled retreat history and that reconstructed from proxy data, particularly regarding the timing of the retreat. This issue needs to be addressed before considering the elevation history of Greenland, as the magnitude and timing of bedrock uplift and ice thickness changes are influenced by the timing of ice margin retreat.

Major concerns I’d like to point out:

• Wrong timing and spatial variability of ice-sheet retreat: The GrIS evolution simulated throughout the deglaciation does not match the timing of the ice margin retreat inferred from several proxies (most of all moraines) and geomorphological reconstructions. Leger et al. 2023 did an amazing work in collecting and processing most of data available across Greenland to provide a detailed information on area change and timing of the deglaciation across the ice sheet. Most of the GrIS deglaciated already before 9 kyr ago (see Figure 15 of Leger et al., 2023) while your figure 6b still shows a glacial-expanded Greenland. This suggests that simulations presented here are systematically late in modeling the ice-sheet shrinkage during the last deglaciation and I am wondering how this might affect the simulated ice-core thinning curves. I think this is a key aspect that has to be solved before investigating the thinning curves in detail, otherwise you might find a good elevation history match but for the wrong reasons. Surface elevation changes occur because of changes in bedrock elevation, but also because of changes in ice thickness, which are in turn related to changes in the mass balance and ice dynamics. Ice dynamics and local changes in bedrock elevation are associated with the retreat history of Greenland, the former mostly due to dynamic reorganization occurring during the deglaciation, and the latter to changes in the ice load. I do believe that simulating a correct bedrock elevation change during the last deglaciation is key to correctly replicate the Holocene surface elevation change. However we should first be able to simulate the right retreat history of Greenland before trying to replicate the Holocene surface elevation curves from ice cores with a model. 
• Importance of bedrock adjustment: I find the procedure for adjusting the bedrock quite interesting but I think I don’t fully understand the motivation/significance of it. What does this strategy mean conceptually: why do you need to refine the bedrock at the LGM, if the model simulates the glaciostatic adjustment? Is this a necessary step to well replicate the modeled surface elevation history throughout the Holocene or is it a strategy to reduce the elevation mismatch only at the present? Have you run an ensemble of simulations without such an adjustment to evaluate its effect on the ice-sheet shrinkage during the last deglaciation? I think that a more detailed justification of this procedure and a discussion about the implications of such a procedure on your results is needed. 
• Outdated paleo climatology: using a spatially homogeneous temperature and precipitation signal to force the model during the last deglaciation is an outdated procedure. There are several products available with high temporal resolution (e.g. Badgeley et al, 2020, Buizert et al., 2018), based on the trace21k experiment (Liu et al., 2009, He et al., 2013) and improved with ice-core derived information that can do the job better. These paleoclimate reanalysis have enhanced our knowledge on the Holocene climate in Greenland and make the anomaly method for this period completely outdated. These products have already been used in several works tackling the GrIS evolution during the last deglaciation (Briner et al., 2020, Cuzzone et al., 2022, very recently in Tabone et al, 2024) and they clearly show that considering a spatially homogeneous temperature across Greenland is incorrect. I believe that using such data would also help to reproduce the spatially variable retreat suggested for the last deglaciation (Leger et al., 2023). 
• Discussion missing on drivers of Holocene thinning: I am missing a clear discussion on what is driving the Greenland elevation change following your simulations. You claim that by including the IIS and an ice shelf in Baffin Bay you can reproduce the Holocene elevation history in Greenland, but which is the glaciological explanation for it? Is it the bedrock response to non-local changes in the ice load, or is it the dynamic effect induced by the loss in buttressing upon ice shelf retreat in Baffin Bay, or both of them or what is responsible for the Holocene elevation drop in Greenland? If the loss of buttressing is the preponderant mechanism, why do we see a drop in elevation already several thousand of years prior to the ice shelf collapse? Please, provide a clearer explanation of your findings as this is key to improve our knowledge in the field.
• Usage of other present-day observables: I understand that the estimation of model uncertainty/model parameter calibration has been done by applying a Bayesian approach only to the ice-core surface elevation histories, whilst the others “observables” of Table 2 (present day thickness, surface velocity) are estimated from the resulting probability density functions. Why not including such “observables” already in the bayesian inference to validate your model simulations? This would allow to not rely only on surface elevation change (which is much more uncertain than satellite inferred observables) to evaluate model results.
• Unclear glacial spin up: Is the model run with a fixed LGM climatology for 80 kyrs or is it run forced with a temporal variable temperature and precipitation signal? If it is the latter, which is the glacial climatology used? Moreover, how are the parameters set up during the spin up? If I understand correctly, only one spin up is done for all the simulations in the ensemble, then the response of the model to the 20 key parameters is explored only since the LGM. If this is so, this might lead to a certain “shock” of the model in the first years after the initialization (after the LGM). Actually, I think I see this shock in Figure 7b, where the volume suddenly drops after the LGM for the “constrained to present geometry” case. I don’t think this procedure undermines your results, since your work focuses on the Holocene (which starts after 8000 years from this “new” initialization), but it would be good either to do one spin up for each ensemble member to avoid this inconsistency or at least describe the first years of your ensemble simulations as part of the spinup. 

Specific comments.

Please check the citation format throughout the text: The sentence “as demonstrated by (Adalgeirsdóttir et al., 2014)” should be “as demonstrated by Adalgeirsdóttir et al. (2014)”, for instance.

Lines 6-8: How can we have “confidence in the modeled GrIS historical evolution” if the model does not reproduce the retreat history correctly?

Lines 38-41: this is not entirely true. What about the constraints on past ice-sheet extent given by e.g. moraines, marine sediment cores, coastal organic material? This is a reliable information we do have on past GrIS retreat.

Lines 44-45: I think here a more detailed description about previous work is needed. It is true that previous work did not consider the effect of the IIS on the deglaciation but it was still successful in reducing the data-model misfit at specific ice-core sites. Lecavalier et al., 2017 was actually successful in reproducing the magnitude of the Camp Century elevation change by correcting the Holocene climate forcing in North Greenland, as reconstructed from the Agassiz record, to force an ice-sheet model coupled to a GIA model of relative sea level change. Another very recent work (Tabone et al., 2024) was able to reduce the mismatch at the NGRIP site by considering the effect of the paleo NEGIS dynamics on the ice-sheet interior using a 3D ice-sheet model. 

Lines 60-70, 114-116: as written in the general comments, using a homogeneous temperature reconstruction from ice cores to force the model is unnecessary here, as there are already several deglaciation climate products (Badgeley et al., 2020, Buizert et al., 2018) which do consider the spatial climatic variability that the GrIS experienced during the last deglaciation. I think using this reanalysis instead of a uniform climatology would allow to better replicate the asynchronous ice-sheet retreat during the last deglaciation (see below and Leger et al., 2023). 

Lines 130-131 and Table 1: ocean melt onset parameter “tau”: why is this explored between -4 and -8 ka, whilst Clark et al., 2020 show an increase in the oceanic forcing already at the early Holocene (at 45°N, 30°W, figure 1k)? Below in the manuscript it is found that the sub-shelf ocean melt that best matches the Camp Century thinning curve starts around -5 ka. Is there any evidence of this, e.g. from sediment cores in Nares Strait? 

Lines 155-156: could you show a 2D map of the bedrock at -20 ka after it has been adjusted? Also, from Figure 2 it seems that the adjusted LGM bedrock elevation is higher than that of the present day, on average, but I believe that it is lower than the one simulated at the LGM before the iteration…

Line 162: “unwanted deglaciation”: is this because by lowering the bedrock elevation at the margin (I guess without updating the ice thickness), the surface elevation decreases too, therefore the ice sheet surface becomes exposed to higher air temperatures?

Figure 3: why not showing the 2D bedrock elevation at the last iteration too, for comparison?

Line 169: “while the oceanic and atmospheric parameters chosen reflect the change in our model setup”: unclear, please rephrase. 

Line 202: “two simulations were restricted from advancing beyond the present-day GrIS coast…”. To my understanding one simulation did not advance from the present-day GrIS coast, and the other from the Greenlandic continental shelf, isn’t that it? Please rephrase. 

Lines 202-206: As far as I understand, you initialize the model at the LGM prior to start the simulation ensemble, so how can it be that you simulate different GrIS extents at the end of the glacial period depending on the parameters set? I would expect that all your simulations show a well advanced LGM, since this configuration is generated by the same initialization, and then in the ensemble you explore different deglaciation patterns. What am I missing?

Lines 230-231: what about the configuration at the LGM? Why do you simulate a maximum glacial extent only at -12 ka? Can this be a drift of the model as it might still be adapting to the new parameters set after the initialization at the LGM? 

Lines 243-244: “The GrIS rate of change becomes negative at -10.7 ka and peaks at -7.8 ka with a mass loss rate of 548 Gt a−1”, is it? Or does it peak around -5 ka? I don’t see such a mass loss rate either. 

Figure 6: please change the color scales so that the bathymetry can’t be confused with the surface velocity. Also, the “bridge” between the LIS and the GrIS is floating ice, right? Please show clearly where are the ice shelves in these maps. 

Figure 11c: it would be interesting to see your simulated bedrock uplift and the ice thinning for the four ice-core sites as in Figure S2 of Lecavalier et al., 2017 or Figure S15 of Tabone et al., 2024. Let’s consider Camp Century for instance. From Figure 11c you simulate an Holocene uplift of ~400-500 m, so to replicate the elevation drop of ~600 m you model in Figure 4, you’d need a decrease in ice thickness more or less of the same magnitude. Is this what you simulate? 

Line 269: please correct “Smith Ice Stream” to “Smith Sound Ice Stream”.

Lines 272-274: how are your isochrones compared with isochrones from Leger et al., 2023? Also, “The modeled collapse of the ice bridge in Nares Strait occurs at 4.9±0.5 ka before present,…”: this is in contrast with several evidence/modeling work (e.g. Figure 15 from Leger et al, 2023, but also Lecavalier et al., 2017, England et al., 2006, )…I think this is a central point that has to be solved before we can use elevation change histories to constrain models. The timing in the bedrock uplift should reflect the timing in the deglaciation, so how can we trust the Holocene thinning curves if the modeled margin retreat lags several thousand of years the observations?

Figure 9c: please choose a discrete color palette to better show the isochrones of the last deglaciation (one color every 2kyr for example). 

Lines 278-280: where do we see this rate of change? Also, do you mean 5000 years (and not 500 years)? At least from your video (https://av.tib.eu/media/68337) I can see by eye that the rate of change seems to follow the timing of the ocean melt scaling, which in your best simulation starts 5000 years ago.  

Lines 281-282: there is a lot of uncertainty around the timing and the magnitude of the oceanic forcing, but we know (1) from evidence inferred from sediment cores, that this might have increased already several thousand of years before the Holocene Thermal Maximum (e.g. Jennings et al., 2017, Lloyd et al., 2023 …) and (2) from paleo model simulations (e.g. Trace-21ka experiments using the NCAR-POP model), that this hasn’t been uniform around Greenland during the last deglaciation. For example, warmer oceanic waters have been suggested to occur in the east/northeast of Greenlandic coasts already at the early Holocene (e.g. Lloyd et al., 2023, Werner et al., 2016, …). In this work, the ocean thermal forcing is scaled depending on the latitude (<71°N, between and > 80°N), but this is a crude simplification. Given the importance of such a forcing in your simulations, I believe you should discuss in more detail the limitations of using a quasi-uniform oceanic forcing and the unrealistic activation of this forcing only in the last 5000 years. Specifically, describe how these factors might affect your results.

Line 283: “the estimated mass loss rate shifts from a prior of -12±40 mm ka−1 to a posterior of -23±26 mm ka−1”: I don’t see the -12±40 mm ka−1 in Table 2. Where does this estimate come from?

Lines 310-312: so, do you think that the assumption made by Lecavalier et al., 2017 was wrong? And if yes, could you explain better why? I think that this whole section (6.3) should be better discussed. Please follow Jessica Badgeley’ comments on the uncertainty in deriving elevation changes from ice cores. I think this is a central point of your discussion: how much can we trust elevation histories only to validate model simulations if we can’t quantify their uncertainties? 

Line 355-358: I agree that an ice-sheet model coupled to an atmospheric one would help to investigate the response of the ice sheet to non-local climatic effects, but there are cheaper solutions that could already improve the representation of the non-uniform temperature and precipitation patterns across Greenland, e.g. the usage of a spatial variable paleo climatology (i.e. Buizert et al., 2018, Badgeley et al., 2020). See my general comments.

Lines 366-367: I believe that “further investigations” are actually needed since paleoclimate has a primary control on the evolution of the GrIS, and potentially on the surface elevation history. See point above. 

Lines 421-427: I find this paragraph really interesting, but it should be better explored in the discussion, not only in the conclusion. For example I wasn’t able to see a proper discussion on the deviation of your modeled uplift rates at the present with respect to the observed ones, besides one sentence in Section 5.4. Again, I think that the delay in the modeled retreat is a central part of the discussion and should be better addressed. 

References:

Badgeley, J. A., Steig, E. J. and Dütsch, M. (2022) Uncertainty in reconstructing paleo-elevation of the Antarctic Ice Sheet from temperature-sensitive ice core records, Geophysical Research Letters, 49, e2022GL100334. https://doi.org/10.1029/2022GL100334

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. (2020) Rate of mass soss 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

Buizert, C., Keisling, B. A., Box, J. E., He, F., Carlson, A. E., Sinclair, G., and DeConto, R. M. (2018) Greenland-wide seasonal temperatures during the last deglaciation, Geophysical Research Letters, 45, 1905–1914, https://doi.org/10.1002/2017GL075601

Cuzzone, Joshua K., et al. "Simulating the Holocene deglaciation across a marine-terminating portion of southwestern Greenland in response to marine and atmospheric forcings." The Cryosphere 16.6 (2022): 2355-2372.

England, J., Atkinson, N., Bednarski, J., Dyke, A., Hodgson, D., and Ó Cofaigh, C.: The Innuitian Ice Sheet: Configuration, Dynamics and Chronology, Quaternary Science Reviews, 25, 689–703, https://doi.org/10.1016/j.quascirev.2005.08.007, 2006.

He, F., Shakun, J. D., Clark, P. U., Carlson, A. E., Liu, Z., Otto-Bliesner, B. L., and Kutzbach, J. E. (2013) Northern Hemisphere forcing of Southern Hemisphere climate during the last deglaciation, Nature, 494, 81–85.

Jennings, Anne E., et al. "Ocean forcing of Ice Sheet retreat in central west Greenland from LGM to the early Holocene." Earth and Planetary Science Letters 472 (2017): 1-13.

Lecavalier, B. S., et al., High Arctic Holocene Temperature Record from the Agassiz Ice Cap and Greenland Ice Sheet Evolution, Proceedings of the National Academy of Sciences, 114, 5952–5957, 2017.

Leger, T. P. M., et al., A Greenland-wide Empirical Reconstruction of Paleo Ice Sheet Retreat Informed by Ice Extent Markers: PaleoGrIS Version 1.0, Climate of the Past, 20, 701–755, https://doi.org/10.5194/cp-20-701-2024, 2024.

Liu, Z., Otto-Bliesner, B., He, F., Brady, E., Tomas, R., Clark, P., Carlson, A., Lynch-Stieglitz, J., Curry, W., Brook, E., Erickson, D., Jacob, R., Kutzbach, J., and Cheng, J. (2009) Transient simulation of last deglaciation with a new mechanism for Bølling-Allerød warming, Science, 325, 310–314

Lloyd, J. M., et al. "Ice-ocean interactions at the Northeast Greenland Ice stream (NEGIS) over the past 11,000 years." Quaternary Science Reviews 308 (2023): 108068.

Tabone I., et al., “Holocene thinning in central Greenland controlled by the Northeast Greenland Ice Stream”. Nat Commun 15, 6434 (2024). 

Werner, Kirstin, et al. "Holocene sea subsurface and surface water masses in the Fram Strait–Comparisons of temperature and sea-ice reconstructions." Quaternary Science Reviews 147 (2016): 194-209.

This manuscript presents simulations of the Greenland ice sheet spanning the last 20,000 years, using elevation histories derived from ice cores as constraints. 20 model parameters were systematically tested using Bayesian inference on an ensemble of 841 simulations. Aside from a glacial spin-up, all simulations use PISM at 10km resolution, which is adequate for this task. The main results is a set of model parameters that was constrained by time-varying reconstructions and the conclusion that a good fit is not possible without allowing the ice sheet to bridge over Nares Strait and connect to the Innuitian ice sheet.

I think this study is timely and highly relevant, even if the data used as constraint is not new (Vinther et al., 2009) and more comprehensive datasets exist. I come back to this point in my comments below. Overall, the manuscript presents work of high quality and the presentation, text and figures, is very good. I agree with the conclusion that long-term, transient trends in ice volume should be considered to accurately and reliably project future mass loss and sea level contribution. This work is a big step in this direction. However, similar points have previously been made MacGregor et al. (2016, doi: 10.1126/science.aab1702) who also point out the importance of the ice bridge across Nares Strait. The manuscript should reference this earlier work.

My criticism is best summarized in four major comments:

1) Use of idealized climate

The simulations use one of five spatially uniform temperature anomalies. Precipitation is based on these anomalies and modified with a simple meridional gradient. I noted, after(!) my own reading, the extensive comments on this issue in the discussions section and will therefore keep my comment brief. However, I believe that a more in-depth discussion of this approach and its implications is needed. For example, most parameter estimates agree within their uncertainty (Table 1), but E_SIA is clearly different for the northern and southern core sites (line 212ff). I believe that this is a symptom of the model fighting systematic bias in the boundary conditions, possibly the climate.

2) Comparison with more recent data

I think a revised version of the manuscript should include a comparison with the independent dataset by Leger et al. 2024.

3) Unclear bedrock optimization

I am not sure if I understood the bedrock optimization routine correctly. Is it correct that it was only performed with one single set of parameters? The modern bedrock topography of the best fit simulation shows a substantially larger deviation from observations (Fig. 11b) and a higher RMSE than the simulation used to initialize the ensemble. Why is this the case and what implications does that have?

4) Conclusions on ice bridge

The two simulations that restrict the GrIS from advancing beyond the present-day coastline or the ECS mask were run with the same parameters as the best fit simulation from the ensemble without any spatial restrictions (line 202f). Can strong conclusions be drawn from such a setup? How can you exclude the possibility that a different climate is compatible with the ice core constraints without the need to limit the ice extent?

Minor comments:

equations 1 and 3: It is not clear how the latitudes were chosen and why they are different.

l 140: Please include an explanation why E_SIA and n_SSA are varied but not E_SSA and n_SIA.

l 221: I think this should read "The northern 'precipitation' parameter, not 'accumulation'.

l 289f: "[..] modeled bedrock topography is very sensitive to the history of the ice load." Please see my comment #3 above.

section 3 could be more explicit about the model resolution.

section 6.6 needs to discuss the implications of using a uniform temperature anomaly forcing, maybe after presenting the findings by Lecavalier et al. (2017) (line 362). Please consider the comment by Jessica Badgeley.

l 355f: I think it is too strong to state that a full coupling of ice sheets to atmosphere and ocean is the only alternative to including spatially non-uniform climate forcing.

figure 5 (and others): The line colors, in particular the two shades of blue, are difficult to distinguish.