the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
| 18 Feb 2025
A comparison of the last two glacial inceptions (MIS 7/5) via fully coupled transient ice and climate modeling
Abstract. Little is known about the evolution of continental ice sheets through the last two glacial inceptions (Marine isotope stages, MIS 7d and MIS 5d). Here, we present the results of a perturbed parameter ensemble of transient simulations of the last two glacial inceptions and subsequent interstadials (MIS 7e-7c, 240–215 ka and MIS 5e-5c, 122–98 ka) with the fully coupled ice/climate model LCice. LCice includes all critical direct feedbacks between climate and ice. As shown herein, it can capture the inferred sea level change (of up to 80 m) of the last two glacial inceptions within proxy uncertainty. One key underlying question of paleoclimate dynamics is the non-linear state dependence of the climate system. Concretely, in a model-centric context, to what extent does the capture of one climate interval in an Earth systems model guarantee capture of another interval? For LCice, the capture of present-day climate is insufficient to predict capture of glacial inception climate, as only a small fraction of ensemble members that performed "well" for present-day captured inception. Furthermore, the capture of inferred sea level change in one inception has weak correlation with the same outcome for the other.
After partial history matching against present-day and past sea level constraints, the resultant NROY (not ruled out yet) ensemble of simulations have a number of features of potential interest to various paleo communities, including the following. (i) In correspondence with the inferred last glacial maximum configuration, the simulated North American ice sheets are substantially larger than the Eurasian ice sheet throughout MIS 5d-MIS 5c and MIS 7d-MIS 7c. (ii) Hudson Bay can transition from an ice-free state to full ice cover (grounded ice) within 2000 years. (iii) The North American and Eurasian ice sheets advanced southward with rates well above 100 m/yr during the penultimate glacial inception and over 70 (Eurasia) and 90 (North America) m/kyr during last glacial inception. (iv) the Laurentide and Cordilleran ice sheets merge in their northern sectors in 13 out of 14 NROY simulations for MIS 7d, contrary to what is assumed from limited geological data. (v) larger ice sheets display a larger lag in the timing of stadial maximum ice volume compared to that of the insolation minimum; the North American ice sheet maximum lags 5.3 ± 0.5 kyrs behind the MIS 7d insolation minimum. Supplemental resources include a dynamic display of ice advance and subsequent retreat for a sub-ensemble of 14 NROY simulations from MIS 5d-5c and MIS 7d-7c.
Competing interests: One author (Lev Tarasov) is a member of the editorial board of "Climate of the Past". The other authors declare that they have no conflict of interest.
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 preprint. The responsibility to include appropriate place names lies with the authors.- Preprint
(3495 KB) - Metadata XML
- BibTeX
- EndNote
Status: final response (author comments only)
-
RC1: 'Comment on egusphere-2025-495', Anonymous Referee #1, 19 Mar 2025
Geng et al. present numerical experiments of the last two glacial inceptions with a coupled ice sheet – climate model. The experiments are thoroughly performed but I found the manuscript uneasy, long and mostly descriptive, with not-so-pretty figures. I found no clear take-home messages from this impressive amount of simulations.
Major comments
1) My main comment is about the general purpose of the manuscript. At present the manuscript is mostly descriptive with almost inexistent discussions on the physical mechanisms at play. Given the multiple sources of uncertainty, a model is not expected to perfectly reproduce past changes and the authors seem fully aware of this. Instead, a model is based on physical assumptions which can be tested in multiple ways and confronted to palaeo data. Geng et al. manuscript mostly discuss simulated ice sheet changes and how it compares with the palaeo data. What can we learn from this? It is not particularly useful to know that some experiments work better than other, what is useful to know is “why some experiments work better and why others do not”. The main justification for using forward coupled ice sheet – climate model experiments instead of an inverse approach is to be able to discuss mechanisms. For example, from the title, I would expect to see some discussions on the difference in terms of mechanisms for the last two glacial inceptions. Differences in oceanic heat transport for the two inceptions and how they are linked (or not) with the Eurasian ice sheet growth (or else)? How can we link the differences in ice sheet geometry with the orbital forcing? Different patten of precipitations and/or temperature for the two inceptions or simply a shift driven by the insolation anomalies? Why the Antarctic ice sheet is so stable for the two time periods? What makes that one experiment is successful and not the other?...etc. I think there is material for many interesting discussions but the aim of the paper has first to be clarified.
2) A less important comment is that we have very few information of the coupling is done. Earlier work is referenced but the manuscript lacks some important clarifications in the methods. Please explain briefly how the downscaling works and how the surface mass balance is computed. Same for the sub-shelf melt (quadratic dependency? Near neighbour interpolation?…). What kind of ice sheet model is GSM? What is its spatial resolution? How friction and calving are computed? Also, it seems that a relative sea level solver is used but how does it work? Does it assume an eustatic homogeneous value at the start of the transient experiments or a perturbation of the present-day geoid is used? It seems that the bathymetry (ocean grid model) is unchanged, but then what happened for the air-sea fluxes and albedo in the Hudson Bay area when it becomes glaciated? This last point is highly relevant since Hudson Bay area results are highlighted in the manuscript.
3) A time-dependent bias correction is used to correct the climate data used to force the ice sheet model. This correction is maximal when the total ice volume is close to present-day and decrease with increasing ice volume. This is motivated by the fact that we cannot guarantee that the biases remain constant in time. While it sounds appealing I find this approach not really justified. The biases can also, in principle, increase with increasing ice volume. For example, a too wet mountainous area might become even wetter if the surface topography increases (same for temperature). Right now I think the reduction of the bias correction with increasing ice volume is unjustified and seems simply like tuning strategy. It would be nice to show simulations with constant bias correction and to properly discuss what are the implications of the chosen methodology on the results. It strikes as odd to claim to have explored thoroughly the parametric uncertainty with a large ensemble when ad-hoc tuning strategies are used (bias dependent correction and additional warming in certain locations).
4) I feel a bit frustrated to have no information on the spun-up ice sheets and climate at 122 ka and at 240 ka. What the initial ice sheets look like – in terms of ensemble mean and standard deviation? Have you evaluated these initial states against palaeo data? Is part of the transient response can be explained by the state of the ice sheet and the climate at the beginning of the simulation? I expect much larger ice sheet after the 7 ka spin-up at 240 ka compared to the ones at 122 ka since both CO2 and Northern Hemisphere summer insolation are lower.
5) I can find no coupled modern simulations in the manuscript. It seems that the calibration against the present-day climate has been run with LOVECLIM only, not including the interactive ice sheets component. This seems illogical since it is plausible that amongst the 90 NROY members some of them would produce unrealistic present-day ice sheets. Perhaps the members that produce the largest ice mass gain for the inceptions does not reproduce observed present-day ice sheets under modern boundary conditions. It would have been useful to run modern / Holocene ice sheet – climate simulations to rule out the members that do not manage to maintain present-day ice sheets.
6) Is albedo downscaled? If not, the temperature provided by ECBilt to GSM has probably a strong imprint of the original albedo. Also if the albedo is indeed computed on the native T21 grid, it can explains the stepwise ice sheet advances and/or retreat, which is then an artificial feature of the model. Please comment / discuss the albedo in the manuscript.
7) One improvement with respect to Bahadory and Tarasov (2018) is the representation of an interactive Antarctic ice sheet component. However there is virtually no information on Antarctica in the paper. Was it expected to have a stable Antarctic ice sheet for the time periods considered? Don’t we expect a retreated ice sheet at 122 ka? Do the coupling affect the large-scale climate through freshwater and/or salt rejection?
8) The quality of the figures can be largely improved. Most of them are blurry (Fig. 1 or 3 in particular) or contain too much information (Fig. 1 can be probably split into two sub-panels).
Minor comments
- P4L100: “4 x acceleration” what is meant here?
- P4L105: ECBilt contains ageostrophic correction terms.
- P4L108: NA is defined only P7. More generally, do we need these abbreviations?
- P5 Fig.2: It would be useful to have the frequency of exchanges (annual, monthly,...)
- P5L110: I am not sure how robust are the results of Lofverstrom and Liakka (2018). We cannot use a model at different resolutions without retuning some parameters. It is likely that the base model of Lofverstrom and Liakka (2018) – meaning the one that has been calibrated – was the one at high resolution.
- P5L123-124: If the ocean cannot be turned to land, what happens to the Hudson Bay area in terms of albedo, air-sea fluxes,...?
- P6L129-131: The daily variability is not taken into account if monthly means are used – or do you use some kind of a paremetrisation for daily variability?
- P6L134: Salinity is not used in the melt equation? For ice sheet points with no corresponding ocean points, do you use some kind of spatial interpolation?
- P6L136-137: If this is a novelty from earlier work it has to be a bit more described… How the ice shelves are handled? Does it make a difference to have the interactive Antarctic ice sheet? Does it impact the Northern Hemisphere simulated climate?
- P6L139: It is not enough to cite a paper in preparation. We need to know how the surface mass balance is computed since it is a critical point for glacial inception.
- P6L141-148: What formula is used to correct the climatic variables? A simple delta method? Also for evaporation and winds? Please provide the equations used (in supplement) to correct the different variables.
- P6L152: What can we learn from this? Event though a bias correction method is used, we still need to impose an ad-hoc correction. This additional correction is here to make the results “pretty” or does it really change climate and ice sheet dynamics? This point deserves further discussion.
- P6L155: The value 1 is not included in this range while it is what has been used as standard in earlier LOVECLIM works. How compatible are these values with previous model evaluations? Also, the climate sensitivity is low with standard LOVECLIM parameters but the ensemble used here is very large and it is possible that it does not need this artificial climate sensitivity increase.
- P7L172: Refer to appendix A2.
- Figure A1: Precipitation units?
- P7L186: How these dates have been selected?
- P7L187-189: With this methodology, the climate and ice sheets are supposed to be in equilibrium while the sea level record shows large changes, especially for MIS7. Assuming equilibrium is a practical choice but it would be useful to have a discussion on how alternative climate and ice sheet at the start of the transient experiments can impact the results.
- P7L190-191: 7 kyr is not enough to spin up ice sheet internal variables (temperature / viscosity).
- P7L191-192: What freshwater flux is given to the ocean in accelerated experiments? Conserving mass would lead to overestimated fluxes (x4 too high).
- P9 Fig. 3: The ensemble members seem to start with very similar ice volumes at 240 and 120 ka. This is very surprising given the fact that the climate is necessary very different (lower Northern Hemisphere insolation and GHG at 240 ka).
- P11 L253-257: This is the only section mentioning the Antarctic ice sheet. What is the geometry of the spun-up ice sheet? The largest shelves are simulated? I find it not convincing to invoke a direct correlation between Southern Hemisphere maximum in insolation and minimum ice volume. Surface melt is low in Antarctica and most of SMB changes are linked to precipitation and oceanic sub-shelf melt. Do you imply that the sub-surface oceanic temperatures are correlated with the local insolation? There is probably a better (positive) correlation between global mean temperature and Antarctic volume (wetter climate). In its current form the Antarctic results are not informative.
- P12 Fig. 4: After an initial ice growth the model melts the Eurasian and North American ice sheets, which is not supported by proxy data. Does it mean that the albedo – melt feedback is too strong?
- P12L266: Specify in the figure that this is the summer isotherm.
- P13 Fig. 6: why not adding the Antarctic ice sheet in this plot – a correlation with Northern Hemisphere insolation is not excluded.
- P13L270: How do you evaluate the sea-level temperature? Which lapse rate is used?
- P16L312-316: Do we have a control of sub-shelf melt? Please discuss.
- P18 End of section 3.3.2: The model simulates the “right” volume but underestimates the extent. Is it because the model is too wet despite the precipitation correction? Or sliding/deformation is underestimated?
- P19 Fig. 10: Dates for ice geometry?
- P20 Sec. 4.1: But then, what is the purpose of Sec. 3.3 where model results are confronted to uncertain geological data? Again for me the strength of a model is to understand mechanisms and processes. Uncertain models are not meant to reproduce uncertain data…
- P21L417 vs. P21L424: Full ice cover in 2 kyr or 1 kyr?
- P21L426: It is really possible to discuss results in this area given the fact that a fixed bathymetry is used and as such the surface mass balance model is evaluated using atmospheric fields above an ocean when it is land.
- P21L436: This is probably due to an artefact of the model since the ice mask is evaluated at the atmospheric model resolution. A whole grid cell can artificially switch from unglaciated to glaciated causing discontinuity in the albedo (hence near-surface air temperature). This result should not be discussed as if it is a real thing.
- P22L456-464: I am not sure that it is related to atmospheric circulation. Albedo is not downscaled which can cause an artificially homogeneous albedo over the whole ice sheet. A sub-grid parametrisation of albedo could also help to reproduce this feature.
- P22L466-472: It looks like a simple repetition of a previous section.
- P23L496: “model’s capability to simulate present-day climate” – it seems to me that this is not completely justified. The coupled ice sheet – climate model has not been used for present-day.
- P23L505-506: It is not really a new results. If some modellers have used deglacial ice sheets for pre-LGM it was for simplicity, but the effect of glacial isostasy is well known and accepted for a while now.
- P24L520: The code is not publicly available, meaning that it is not possible for other scientists to have a look of the coupling strategy in LCice.
Technical comments
- Abstract: problem with the units - 90 m/kyr
- Abstract: remove last sentence
- P2L45: “between the two” → 4 sea level reconstructions are discussed.
- Fig. 1: show x- and y-axis units. Figure is poor quality and really hard to read. To improve the readibility you could to separate CO2 and insolation from the sea level reconstructions. It is very blurry right now.
- Fig. 4 & 5: Capital “m”
- P15L308: Typo – inception
- P25 Fig. A1: Units of precip? Change the color palette so that it appears while when there is no change.
- P27 Tab. B1: Add a short description of these parameters and their units.
Citation: https://doi.org/10.5194/egusphere-2025-495-RC1 -
AC3: 'Reply on RC1', Marilena Geng, 28 May 2025
Thank you for the thorough review and giving us lots of food for thought. Here are our responses:
Summary of key changes we are proposing to address RC1’s comments:
- Including a section analysing the difference between runs that pass our inception filters and those that didn’t for a couple metrics (North American vs Eurasian ice sheet size, AMOC strength, ...)
- Submitted new runs with full bias correction to compare to our transient bias correction simulations
- Submitted new runs without dynamic Antarctic ice sheet coupling to analyse influence of Antarctic freshwater on Northern Hemisphere ice sheets
- Inclusion of more details of the GSM
In Detail:
Major:
1) My main comment is about the general purpose of the manuscript. At present the manuscript is mostly descriptive with almost inexistent discussions on the physical mechanisms at play. Given the multiple sources of uncertainty, a model is not expected to perfectly reproduce past changes and the authors seem fully aware of this. Instead, a model is based on physical assumptions which can be tested in multiple ways and confronted to palaeo data. Geng et al. manuscript mostly discuss simulated ice sheet changes and how it compares with the palaeo data. What can we learn from this? It is not particularly useful to know that some experiments work better than other, what is useful to know is “why some experiments work better and why others do not”. The main justification for using forward coupled ice sheet – climate model experiments instead of an inverse approach is to be able to discuss mechanisms. For example, from the title, I would expect to see some discussions on the difference in terms of mechanisms for the last two glacial inceptions. Differences in oceanic heat transport for the two inceptions and how they are linked (or not) with the Eurasian ice sheet growth (or else)? How can we link the differences in ice sheet geometry with the orbital forcing? Different patten of precipitations and/or temperature for the two inceptions or simply a shift driven by the insolation anomalies? Why the Antarctic ice sheet is so stable for the two time periods? What makes that one experiment is successful and not the other?...etc. I think there is material for many interesting discussions but the aim of the paper has first to be clarified.
Some of these questions can not be answered in this study. Questions about mechanism and whether differences between the two inceptions are linked to insolation anomalies or different patterns in precipitation/temperature (which would be triggered by insolation changes, since orbital and GHG forcing are the only differences in model input) will need to be answered in a sensitivity analysis (in progress) where forcings and feedbacks are isolated.
The purpose of the paper is as stated in the introduction: 1) Is there a correlation between PD performance and inception performance? Between the last and penultimate glacial inception performance? 2) Is ice evolution during an inception always following the same pattern? Or are there differences between the last and penultimate inception? Are there differences within the ensemble for each inception? 3) Do the simulations match the limited data? Where do they not match, and do these non-matches teach us anything new?
We will consider clarifying these points at the end of the introduction, so no false expectations are raised. We will also add an extra section on what differentiates good from bad runs (with the answer being that there is no clear answer): bad runs have both little NA and EA ice volume (i.e. limited EA ice growth isn’t the only key issue), there is no significant difference in AMOC strength between good and bad runs.
2) coupling: explain briefly how the downscaling works and how the surface mass balance and sub-shelf melt is computed. “what kind of ice sheet model” is GSM? Spatial resolution? Friction, calving? How does relative sea level solver work? What happened for the air-sea fluxes and albedo in the Hudson Bay area when it becomes glaciated (since land-sea mask can’t change)?
We have added a few more details about the GSM, but full details are already provided in the cited discussion paper (which was erroneously listed as in preparation but was actually already under open review in GMD): “Since the first publication using LCice (Bahadory and Tarasov, 2018), we have updated the coupler, the most significant update of which is the inclusion of a dynamic Antarctic ice sheet. The addition of the latter is subject to the same coupling as described in Bahadory and Tarasov (2018) for the Northern hemispheric ice sheets. However, given the strong (and turbulent) circumpolar ocean current around Antarctica and the relatively coarse resolution of CLIO, we assume the spatial distribution of freshwater discharge from Antarctica will have little impact within the model. Therefore, discharge is uniformly spread around the perimeter of the ice sheet.”
“The GSM has also been updated, including the conversion from pure shallow ice approximation to hybrid shallow ice/shallow shelf ice dynamics (Pollard et al., 2015), and the introduction of a novel (and physically motivated) accounting for the impact of changing insolation forcing on surface mass-balance (Tarasov et al., 2025, with surface insolation computed indepently by the GSM and accounting for monthly and orbital dependence of solar angle). The GSM includes full thermomechanical coupling, global visco-elastic glacial isostatic adjustment with a 0-order approximation to compute spatial variations in the Geoid, warm-based basal drag dependence on surface sediment cover and roughness, a pro-glacial lake resolving surface drainage solver, and both marine and lacustrine ice calving. Marine calving accounts for meltwater enhancement of stress-balanced crevasse propagation and accumulated strain, with a marine ice cliff instability enabled for ice marginal grid cells. For the experiments herein, the GSM is run at 0.5° longitude by 0.25° latitude grid resolution except for a 20 km polar stereographic grid for Antarctica.”
“The constant land mask is a major limitation for many advanced paleo climate models (e.g. Meccia and Mikolajewiez, 2018). The dynamic Bering Strait parametrization of LOVECLIM at least captures the most significant ocean gateway change during a glacial cycle. However other changes in ocean gateways (e.g. with the likely next most critical being Nares Strait) are not captured. The impact of constant land mask on surface fluxes is at least partially compensated for by sea-ice cover.”
3) [...] the reduction of the bias correction with increasing ice volume is unjustified and seems simply like tuning strategy. It would be nice to show simulations with constant bias correction and to properly discuss what are the implications of the chosen methodology on the results. It strikes as odd to claim to have explored thoroughly the parametric uncertainty with a large ensemble when ad-hoc tuning strategies are used (bias dependent correction and additional warming in certain locations).
The bias correction dependence on glacial state is under ensemble parameter control, and many of the simulations have moderate to no reduction in bias correction for the full glacial (LGM) state (which unfortunately was not made clear in the original submission). We would also counter that there is no justification for preservation of full bias correction through a glacial cycle with Loveclim. As we explicitly state near the end, bias corrections are a bitter pill. We did now set up simulations that use the same parameter vectors as the NROY ensemble but have full constant bias correction. Their results will be compared to the transient bias correction and added to our analysis in our revised submission.
4) no information on the spun-up ice sheets and climate at 122 ka and at 240 ka. What the initial ice sheets look like – in terms of ensemble mean and standard deviation? Have you evaluated these initial states against palaeo data? Is part of the transient response can be explained by the state of the ice sheet and the climate at the beginning of the simulation? I expect much larger ice sheet after the 7 ka spin-up at 240 ka compared to the ones at 122 ka since both CO2 and Northern Hemisphere summer insolation are lower
Summer insolation (mid-July at 65 N) at 240 ka is higher than at 122 ka, but CO2 concentration is lower at 240 ka than at 122 ka. Based on this, there is no basis we see to a priori judge whether there should be larger or smaller ice sheets after the 240 ka vs 122 ka spin-up. In Figures 4 and 5 (time series of ensemble mean and standard deviation of ice volume and area by individual ice sheet), the initial ice sheet size is shown. Based on the Figures, North American and Eurasian ice sheets start with close to zero ice, and there is virtually no difference between the last and penultimate inception of ice sheets. Greenland and Antarctica increase their volume and area over the course of the spin-up.
5) I can find no coupled modern simulations in the manuscript. It seems that the calibration against the present-day climate has been run with LOVECLIM only, not including the interactive ice sheets component. This seems illogical since it is plausible that amongst the 90 NROY members some of them would produce unrealistic present-day ice sheets. Perhaps the members that produce the largest ice mass gain for the inceptions does not reproduce observed present-day ice sheets under modern boundary conditions. It would have been useful to run modern / Holocene ice sheet – climate simulations to rule out the members that do not manage to maintain present-day ice sheets
What would be the purpose of this? The glacial cycle configuration of LCice imposes full bias correction at present-day. If the intent is that this should be turned off, we do not see what value this would add here. Our present-day filtering is for LOVECLIM parameters. To couple with the GSM would entail at least a 10-fold increase in the ensemble size, given that GSM parameters would also need to be probed. Furthermore, our approach mirrors the context for state-of-the-art climate models that are tuned for present-day climate without interacting ice sheets.
6) how is albedo downscaled form EC Bilt to GSM?
It is not downscaled and would not make sense given the large difference in grid resolutions. Figure 2 is a complete diagram of couplings and does not show any albedo nor short wave couplings. Instead, the GSM independently computes top of the atmosphere insolation as a function of month and year, and resultant surface insolation as a function of solar angle and assumed mean atmospheric transmissivity. To make this clear, we’ve added the following text: “The GSM has also been updated, including the conversion from pure shallow ice approximation to hybrid shallow ice/shallow shelf ice dynamics (Pollard et al., 2015), and the introduction of a novel (and physically motivated) accounting for the impact of changing insolation forcing on surface mass-balance (Tarasov et al., 2025, with surface insolation computed independently by the GSM and accounting for monthly and orbital dependence of solar angle).”
7) One improvement with respect to Bahadory and Tarasov (2018) is the representation of an interactive Antarctic ice sheet component. However there is virtually no information on Antarctica in the paper. Was it expected to have a stable Antarctic ice sheet for the time periods considered? Don’t we expect a retreated ice sheet at 122 ka? Do the coupling affect the large-scale climate through freshwater and/or salt rejection?
As we will make clearer in the revised document, Antarctica has high grounding-line sensitivity to the LCice ensemble parameters. Yes, ideally, we would want some AIS retreat at 122 ka. However, GSM parameter tuning attempts persistently resulted in present-day Antarctica having excessive grounding line retreat when parameters were tuned for 122 ka mass loss.
8) figure quality/overload of info in figure 1
We’ve split Figure 1 into 2 figures as suggested to avoid overload. In our version, the figures are not blurry.
Minor:
- P4L100: “4 x acceleration” what is meant here?
Added a clarification to the end of section 2.2: “In this context, "acceleration" refers to asynchronous coupling between the climate and ice sheet components: under a 4x acceleration, five years pass in the climate model while twenty years pass in the ice sheet model. The orbital and greenhouse gas forcings applied to the climate are effectively condensed, i.e., 20 years of forcing are applied over just 5 years of climate simulation.”
- P4L105: ECBilt contains ageostrophic correction terms.
Not sure what the referee means here? Goosse et al. (2010), like us, refer to ECBilt as “quasi-geostrophic” in the description of the model. This should imply that the model is not purely geostrophy but does contain ageostrophic terms, as the referee notes.
- P4L108: NA is defined only P7. More generally, do we need these abbreviations?
Fixed
- P5 Fig.2: It would be useful to have the frequency of exchanges (annual, monthly,...)
The GSM/LOVECLIM coupling timesteps are already described in the text. We’ve added the LOVECLIM internal coupling timesteps: “ECBilt and CLIO are coupled once a day, VECODE is coupled to ECBilt once a year.”
- P5L110: I am not sure how robust are the results of Lofverstrom and Liakka (2018). We cannot use a model at different resolutions without retuning some parameters. It is likely that the base model of Lofverstrom and Liakka (2018) – meaning the one that has been calibrated – was the one at high resolution.
Fair concern, however, we don’t think that retuning would affect the resolution of atmospheric wave dynamics that Lofverstrom and Liakka (2018) found to be critical for this resolution dependence. Furthermore, they indicate that only the T21 version they tested had not been (fully?) retuned.
- P5L123-124: If the ocean cannot be turned to land, what happens to the Hudson Bay area in terms of albedo, air-sea fluxes,...?
added to model description: “The constant land mask is a major limitation for many advanced paleo climate models (e.g. Meccia and Mikolajewicz, 2018). The dynamic Bering Strait parametrization of LOVECLIM at least captures the most significant ocean gateway change during a glacial cycle. However, other changes in ocean gateways (e.g. with the likely next most critical being Nares Strait) are not captured. The impact of constant land mask on surface fluxes is at least partially compensated for by sea-ice cover.”
and to results:
“Ice appears first in the North of Hudson Bay, including the connection to Hudson Strait, after which the rest of Hudson Bay can glaciate within 2 kyr. Since the land-sea mask cannot change during the simulation, the ocean does not see the ice sheet. However, sea ice appears in glaciated areas in the ocean component and limits atmosphere-ocean interaction.”
- P6L129-131: The daily variability is not taken into account if monthly means are used – or do you use some kind of a paremetrisation for daily variability?
As stated in Figure 2, the standard deviation of monthly mean temperature is taken into account. We have corrected this now to indicate that this is also the case for wind velocities. We now make clear that this standard deviation accounts for daily variations within a given month over the ice sheet coupling interval. Daily variations in precipitation will have much less impact and therefore are ignored. Lapse rates are known to significantly vary diurnally, but the fundamental problem is that the required lapse rate (2 m air temperature response to changes in surface elevation) is not readily available and is crudely approximately by the free air lapse rate from Loveclim. The revised section now reads: “The pre-existing LOVECLIM components are coupled to the Glacial Systems Model GSM (Tarasov et al. 2025). The coupler includes all important feedbacks between climate and ice for a glacial cycle context (Figure 2) except potentially dust (Willeit et al., 2022). It passes the monthly mean and standard deviation of temperature and wind, and monthly means of precipitation, evaporation and lapse rate from ECBilt to the GSM. The standard deviation is computed on the LOVECLIM timestep to account for both diurnal and daily variability within a month.”
- P6L134: Salinity is not used in the melt equation? For ice sheet points with no corresponding ocean points, do you use some kind of spatial interpolation?
The submarine melt model is detailed in the GSM paper and assumes a constant salinity. We would argue that uncertainties in subshelf circulation and temperature swamp out the missing impact of changes in salinity. All potentially marine sites are attached to an upstream associated ocean temperature profile extracted from CLIO output. Though the original submissions mentions “profiles” it wasn’t clear enough about this coupling, and has now been remedied: “To facilitate dealing with complicated ocean grids as well as enable ocean temperature forcing where the land mask has a terrestrial value, ocean temperature profiles from CLIO for computing submarine melt in the GSM are extracted at key index sites and then applied to an assigned downstream region (Bahadory and Tarasov 2018).”
- P6L136-137: If this is a novelty from earlier work it has to be a bit more described… How the ice shelves are handled? Does it make a difference to have the interactive Antarctic ice sheet? Does it impact the Northern Hemisphere simulated climate?
We have now submitted our NROY parameter vectors without dynamic Antarctic ice sheet coupling. Given the fixed ocean land mask of LOVECLIM, we expect the main influence of the Antarctic ice sheet in LCice to be freshwater input. We will discuss the role of this on NH ice evolution when the runs are complete. We do not understand the query about how ice shelves are handled. The GSM paper fully describes how ice shelves are handled, but this should already be clear from the stated “hybrid shallow ice/shallow shelf ice dynamics”. With respect to the coupling, this follows from the stated constant land/sea mask. Freshwater coupling from GSM to LOVECLIM is the same as for other ice sheets; otherwise, we would have stated so. However, as the freshwater injection site for Antarctica isn’t included in the original LCice description paper, we added a brief description already presented in response to major comment 2.
- P6L139: It is not enough to cite a paper in preparation. We need to know how the surface mass balance is computed since it is a critical point for glacial inception.
This was in error. The cited paper was already under discussion on GMD, we’ve fixed the citation now: Tarasov, L., Lecavalier, B. S., Hank, K., and Pollard, D.: The glacial systems model (GSM) Version 24G, Geosci. Model Dev. Discuss. [preprint], https://doi.org/10.5194/gmd-2024-175, in review, 2025.
- P6L141-148: What formula is used to correct the climatic variables? A simple delta method? Also for evaporation and winds? Please provide the equations used (in supplement) to correct the different variables.
Added in supplements as requested for temperature and precipitation. There is no correction for evaporation or winds.
T_bias_corrected(t) = T_LOVECLIM(t) + max(0,(1 - I_s^(rBiasT)) x (T_ERA5(0 ka) – T_LOVECLIM(0 ka))
P_bias_corrected(t) = P_LOVECLIM(t) x (P_avg_ERA5(0 ka)/P_LOVECLIM(0 ka))^(rBiasP*max(0,1-I_s))
with sea level index I_s and bias parameter rBiasT/rBiasP
- P6L152: What can we learn from this? Even though a bias correction method is used, we still need to impose an ad-hoc correction. This additional correction is here to make the results “pretty” or does it really change climate and ice sheet dynamics? This point deserves further discussion.
We state that excessive ice over Alaska and Siberia “ would otherwise distort results for adjacent regions”. Large excess ice will buttress adjacent ice and impact atmospheric wave dynamics and summertime surface albedo, so compensating for this model limitation seems sensible. We have now also added a more detailed discussion on the regional impacts of limited atmospheric model resolution with Alaska and Siberia (southern Laurentide) being explicitly identified as problematic regions in the methods section: “[...] Lofverstrom and Liakka (2018) find that T42 is the minimum resolution to grow inferred LGM EA ice extent, though LGM NA ice extent can be reasonably captured at T21 (except for excessive glaciation of Alaska, which is not present in the T31 simulation). However, even at T21 there is excessive JJA (June to August) heat flux convergence by stationary waves over the southern half of the Laurentide ice sheet in comparison to their T85 simulations (Lofverstrom and Liakka, 2018). As such, the southern glacial ice sheet extents for the Laurentide portion of the North American ice sheet complex are likely to be underestimated in LCice. Simulated LGM Siberian summertime air temperatures anomalies (relative to present) have also been shown to be highly sensitivity to changes in the atmospheric model physics (from CAM4 to CAM5) in experiments with CESM version 1.2 as well a having large variance among simulation from PMIP 2 and 3 paleoclimate model intercomparisons (Bakker et al., 2020).”
- P6L155: The value 1 is not included in this range while it is what has been used as standard in earlier LOVECLIM works. How compatible are these values with previous model evaluations? Also, the climate sensitivity is low with standard LOVECLIM parameters but the ensemble used here is very large and it is possible that it does not need this artificial climate sensitivity increase.
None of the parameter vectors using a GHG factor of 1 grow enough ice during inception to pass our filtering.
- P7L172: Refer to appendix A2
Done
- Figure A1: Precipitation units?
cm/year, added to figure
- P7L186: How these dates have been selected?
Based on sea level low stand, rephrased to: “Penultimate glacial inception simulations start at 240 ka, last glacial inception at 122 ka, the timing of respective sea level low-stands close to PD (Spratt and Lisiecki, 2016)”
- P7L187-189: With this methodology, the climate and ice sheets are supposed to be in equilibrium while the sea level record shows large changes, especially for MIS7. Assuming equilibrium is a practical choice but it would be useful to have a discussion on how alternative climate and ice sheet at the start of the transient experiments can impact the results.
This methodology is not intended for equilibration, and nowhere is that indicated in the text. The issue of initial state dependence will be much more appropriately addressed in a subsequent submission, where we will compare last glacial cycle simulations against the last two glacial cycle simulations. We will add a brief discussion of this issue to the revised text.
- P7L190-191: 7 kyr is not enough to spin up ice sheet internal variables (temperature / viscosity).
The GSM has its own internal thermodynamic spin-up routine as documented in Tarasov et al. (2025) for simulation startup, and this is now stated in the revised text. However, the challenge is that the ice/climate system is never in equilibrium. Thus, we do not aim for equilibration. That is why we explicitly state “spin-up” as opposed to “equilibration”.
- P7L191-192: What freshwater flux is given to the ocean in accelerated experiments? Conserving mass would lead to overestimated fluxes (x4 too high).
Mass is not conserved. The flux is either preserved or enhanced by a factor of 2 as an intermediate between conserved flux and conserved mass. We have added the following text to the revision: “Freshwater volumes are not conserved as resulting fluxes would likely trigger excessive response of the meridional overturning circulation. Instead, depending on the ensemble parameter vector, freshwater flux from the GSM is either conserved or enhanced by a factor of up to 2 (in which case freshwater volume is still reduced by 50% given the 4X acceleration of LOVECLIM). Furthermore, this enhancement is only applied to the signed freshwater flux due to changes in ice mass. Given challenges with LOVECLIM internal treatment of freshwater routing, the internal routing within LOVECLIM of precipitation minus evaporation is retained. As such, this component of the flux is conserved.”
- P9 Fig. 3: The ensemble members seem to start with very similar ice volumes at 240 and 120 ka. This is very surprising given the fact that the climate is necessary very different (lower Northern Hemisphere insolation and GHG at 240 ka).
See major comment 4: insolation is higher but CO2 is lower at 240 ka than at 122 ka. Under these spin-up conditions, the ice NA and EA ice sheets show only minimal to no growth during spin-up for both periods.
- P11 L253-257: This is the only section mentioning the Antarctic ice sheet. What is the geometry of the spun-up ice sheet? The largest shelves are simulated? I find it not convincing to invoke a direct correlation between Southern Hemisphere maximum in insolation and minimum ice volume. Surface melt is low in Antarctica and most of SMB changes are linked to precipitation and oceanic sub-shelf melt. Do you imply that the sub-surface oceanic temperatures are correlated with the local insolation? There is probably a better (positive) correlation between global mean temperature and Antarctic volume (wetter climate). In its current form the Antarctic results are not informative.
The GSM has a hydro-fracturing representation, and therefore, the enhanced summertime insolation might be able to induce warm enough temperatures to drive at least partial shelf collapse. However, we will carry out a detailed analysis to determine what process(es) directly drove the WAIS retreat, especially given the strong role of submarine melt for the current AIS. The revised submission will have more analysis of the AIS. This will be facilitated by analysis of simulations of NROY parameter vectors with Antarctic ice sheet coupling turned off that are currently in progress.
- P12 Fig. 4: After an initial ice growth the model melts the Eurasian and North American ice sheets, which is not supported by proxy data. Does it mean that the albedo – melt feedback is too strong?
We are curious what proxy data the referee is referring to, especially for interstadial states. The ice melt happens during the interstadial MIS 7c/5c and is expected.
- P12L266: Specify in the figure that this is the summer isotherm.
Fixed in figure descriptions.
P13 Fig. 6: why not adding the Antarctic ice sheet in this plot – a correlation with Northern Hemisphere insolation is not excluded.
Fair enough, this will be added in the revised submission.
- P13L270: How do you evaluate the sea-level temperature? Which lapse rate is used?
As already stated in the submission, “ The temperature downscaling uses the LOVECLIM vertical atmospheric temperature gradient (lapse rate, Bahadory and Tarasov, 2018)”. This has now been clarified to include that there is a 3:1 weighting of the LOVECLIM lapse rate field (with lon, lat, and month dependence) with a nominal mean global value of 6.5K/km.
- P16L312-316: Do we have a control of sub-shelf melt? Please discuss.
We’ve added a tabular list of GSM parameters used to the supplement and their ranges. That will now make clear that there are four GSM ensemble parameters controlling submarine melt for the Northern Hemisphere ice sheets (along with 6 more regional ocean temperature bias parameters for the major ice shelves of Antarctica). This will also include a short description of how the GSM parameter ranges and values were chosen (strongly informed by history matchings of the four last glacial cycle major ice sheets).
- P18 End of section 3.3.2: The model simulates the “right” volume but underestimates the extent. Is it because the model is too wet despite the precipitation correction? Or sliding/deformation is underestimated?
The model doesn’t necessarily underestimate the extent given geological age uncertainties (as stated in text “the measure of confidence in this ice reconstruction is low owing to a shortage of geological constraints”, to which we’ve now appended “and associated age uncertainties”). As now described in the revised supplement, the ensemble parameter ranges for hard and soft bed sliding are biased towards the high end of values from those in the NROY (not ruled out yet) set of the history matchings for the corresponding ice sheets. Furthermore, given the large range of sea level uncertainty for glacial inception, it is unclear what the “right ice volume” is. But the issue of primary controls and sources of uncertainty on southern extent deserves more attention. In general the foremost control on terrestrial ice extent is physically 2 meter air temperature (as we’ve shown by the bounding of ice margins by sea-level isotherms) . The core issue is limitations from the atmospheric grid resolution. To make this clearer, we’ve have expanded the discussion of resolution limitations based on the cited Lofverstrom and Liakka (2018) article in the methods section:
“[...] Lofverstrom and Liakka (2018) find that T42 is the minimum resolution to grow inferred LGM EA ice extent though LGM NA ice extent can be reasonably captured at T21 (except for excessive glaciation of Alaska, which is not present in the T31 simulation). However, even at T21 there is excessive JJA (June to August) heat flux convergence by stationary waves over the southern half of the Laurentide ice sheet in comparison to their T85 simulations (Lofverstrom and Liakka, 2018). As such, the southern glacial ice sheet extents for the Laurentide portion of the North American ice sheet complex are likely to be underestimated in LCice. Simulated LGM Siberian summertime air temperatures anomalies (relative to present) have also been shown to be highly sensitivity to changes in the atmospheric model physics (from CAM4 to CAM5) in experiments with CESM version 1.2 as well a having large variance among simulation from PMIP 2 and 3 paleoclimate model intercomparisons (Bakker et al., 2020).”
- P19 Fig. 10: Dates for ice geometry?
Added clarification in figure title: “Two example NROY ensemble members with similar global ice volume at MIS 5d (~ 112 ka, top left) but different ice volume distribution among the ice sheets (top right) and ice sheet geometry (bottom, at 112 ka)”
- P20 Sec. 4.1: But then, what is the purpose of Sec. 3.3 where model results are confronted to uncertain geological data? Again for me the strength of a model is to understand mechanisms and processes. Uncertain models are not meant to reproduce uncertain data…
see major comment 1. If the interest is purely mechanisms and processes, then why bother with paleo? Why not idealized high-resolution sensitivity experiments? Yes, disentangling mechanisms and feedback loops is one key role of modelling, but it is not the sole role. Just because data has uncertainty (as is the case for all paleo data), doesn’t mean it has no constraint or comparative value. And some aspects of the Lcice are likely to be more robust than others.
- P21L417 vs. P21L424: Full ice cover in 2 kyr or 1 kyr?
2 kyrs, fixed.
- P21L426: It is really possible to discuss results in this area given the fact that a fixed bathymetry is used and as such the surface mass balance model is evaluated using atmospheric fields above an ocean when it is land.
The ocean is covered in sea ice when glaciated. We’ve added: “Ice appears first in the North of Hudson Bay, including the connection to Hudson Strait, after which the rest of Hudson Bay can glaciate within 2 kyr. Since the land-sea mask cannot change during the simulation, the ocean does not see the ice sheet. However, sea ice appears in glaciated area in the ocean component and limits atmosphere-ocean interaction.”
- P21L436: This is probably due to an artefact of the model since the ice mask is evaluated at the atmospheric model resolution. A whole grid cell can artificially switch from unglaciated to glaciated causing discontinuity in the albedo (hence near-surface air temperature). This result should not be discussed as if it is a real thing.
We were trying to describe another (likely real) phenomenon that will lead to sudden jumps. We now clarify: “Advance and retreat rates are derived from the latitude of the southernmost ice extent during growth and melt phase along a single longitude transect. If the ice advances (retreats) primarily in the east–west direction rather than north–south, it can suddenly appear (disappear) along the evaluated longitude. This may result in an artificially high calculated advance or retreat rate.”
P22L456-464: I am not sure that it is related to atmospheric circulation. Albedo is not downscaled which can cause an artificially homogeneous albedo over the whole ice sheet. A sub-grid parametrisation of albedo could also help to reproduce this feature.
Fair point. We’ve added to the discussion: “[...] Furthermore, the atmosphere's coarse resolution can render the narrow gap between the Cordilleran and Laurentide ice sheets sub-grid. As a result, albedo changes in this region aren't resolved, potentially leading to artificially high albedo and spurious ice growth.” A future sensitivity analysis could test the influence of (sub-grid) albedo changes on ice evolution
- P22L466-472: It looks like a simple repetition of a previous section.
Not fully, but agreed, too much. It will be condensed and perhaps consolidated with another discussion topic.
- P23L496: “model’s capability to simulate present-day climate” – it seems to me that this is not completely justified. The coupled ice sheet – climate model has not been used for present-day.
None of the state-of-the-art Earth System Models that participated in the last CMIP model intercomparison were run with coupled ice sheet models, and yet they are being assessed and intercompared with respect to their simulation of present-day climate, as we have done for LOVECLIM, though using an assessment geared for glacial cycle contexts. As such, running without ice sheet-model coupling is more appropriate, especially given the tentative conclusion of our results, i.e. just because a climate model captures present-day climate (in whatever sense), doesn’t imply it can capture glacial inception when coupled to a climate model.
- P23L505-506: It is not really a new results. If some modellers have used deglacial ice sheets for pre-LGM it was for simplicity, but the effect of glacial isostasy is well known and accepted for a while now.
When the GIA community continues to largely ignore this issue, then the claim “is well known and accepted” does not seem to be the case. Furthermore, we do not understand the singling out of glacial isostasy. The history dependence of the coupled system will depend on glacial isostasy, but also on the ocean state, the spatial pattern of insolation, and the thermal state of the ice sheet. And no one has previously shown this result with a fully coupled ice and climate model that resolves both atmospheric and ocean circulation for this interval.
- P24L520: The code is not publicly available, meaning that it is not possible for other scientists to have a look of the coupling strategy in Lcice.
The coupling is documented in Bahadory and Tarasov (2018), which includes a linked code archive. We will create an updated archive to capture the few changes to the coupling and the associated version of the GSM.
Technical comments
- Abstract: problem with the units - 90 m/kyr
fixed: 90 m/yr
- Abstract: remove last sentence
Done.
- P2L45: “between the two” → 4 sea level reconstructions are discussed.
fixed
- Fig. 1: show x- and y-axis units. Figure is poor quality and really hard to read. To improve the readibility you could to separate CO2 and insolation from the sea level reconstructions. It is very blurry right now.
We’ve split Figure 1 into two and added axis labels
- Fig. 4 & 5: Capital “m”
fixed
- P15L308: Typo – inception
fixed
- P25 Fig. A1: Units of precip? Change the color palette so that it appears white when there is no change.
Units are added.
- P27 Tab. B1: Add a short description of these parameters and their units.
Description of parameters is in Tab. 2 in supplements. The large majority of parameters are unitless. Only “ahu” -ocean horizontal viscosity – has (expected) unit [m^2/s].
Citation: https://doi.org/10.5194/egusphere-2025-495-AC3
-
AC3: 'Reply on RC1', Marilena Geng, 28 May 2025
-
RC2: 'Comment on egusphere-2025-495', Anonymous Referee #2, 31 Mar 2025
Review of "A comparison of the last two glacial inceptions (MIS 7/5) via fully
coupled transient ice and climate modeling" by Geng et al.The manuscript presents the result of two glacial inception periods using climate and ice sheet coupled model LCice version 2.0, composed of climate module LOVECLIM and ice sheet module GSM. The model simulations were focused on uncertain parameters, primarily in a climate model, and evaluated the results based on simulated sea level changes during two glacial inceptions. The authors analyzed the extent of ice sheets, primarily in the Northern Hemisphere, changes during two glacial inceptions and retreats in the following insolation changes, and they analyzed the speed of ice sheet migration and retreats.
While I think the topic of the study and the methodology are well-suited for Climate of the Past, I could not recommend accepting the manuscript in the current form because of concerns in the model settings:Concerns about Bias corrections
A significant point in the experimental design is that the climate model's bias correction term depends on the climate state, as in L144-145. While it is not trivial that the biases in climate models are the same during glacial periods, the assumption of zero bias in the LGM state is somewhat subjective. More importantly, the bias correction, depending on climate, can act as feedback that does not exist in the real world. Let us assume there is a warm bias of 6 K in the NA ice sheet margin (Figure A1 top). If the ice sheet were to grow from this state to MIS7d, where ice sheet volumes are about half of the LGM (the paper does not state how many meters this is), the weakening of the bias correction would cause the temperature at that point to rise by 3°C. Given that the global mean temperature at the LGM differs by about 4.5±0.9°C (e.g. Annan et al., 2022), wouldn't this feedback from the bias correction significantly affect comparable radiative forcing of insolation and Greenhouse gases? Another issue related to this bias correction is ad-hoc temperature increases in these two regions (L152) It is unclear from the text whether this can be achieved simultaneously with the above bias correction.
This climate-dependent bias correction makes it challenging to make comparisons with other studies, such as the MIS7 (Choudhury et al. 2020), glacial inception experiments without bias correction (we recommend citing Kageyama et al. 2004, Gregory et al. 2012), and glacial cycle experiments performed with bias correction throughout simulations (Ganopolski et al. 2010, Abe-Ouchi et al. 2013). One common explanation for using identical bias correction independent from climate state is that the model bias comes from limited horizontal resolution or the limited climate processes unique to the climate model.
I recommend the authors to have a set of simulations with identical bias correction throughout the simulations as a reference simulation.Specific comments
L2-3: The sentence might be accurate, but as this study uses sea level changes during MIS7 and MIS5 as constraints of the NROY simulations, this statement can be more specific based on the method of the simulations.L4: It can be clarified that insolation and GHG changes are used as transient forcing in the experiments.
L6-10: What does "model-centric context" mean in this sentence? According to Table B1, most of the parameters after MIS7 and MIS5 constraints are within range of the original parameters adapted in the original LOVECLIM.
L16: please indicate the uncertainty range as in Table 2.
L66: Studies on glacial inception using coupled ice sheet-climate models (Kageyama et al., 2004; Gregory et al., 2012) should be referred to.
L 96: As the model name is "LCice2.0," please clarify the major differences from LCice version 1.0 (Bahadory Tarasov 2018; 2021).
L98: which insolation data is used?
L99: Bereiter et al. (2015) include only CO2 time-series. What about other greenhouse gas components? For example, Choudhury et al., (2020) used three GHG components (CO2, CH4, N2O) from Dome C ice core (Loulergue et al., 2008; Lüthi et al., 2008; Schilt et al., 2010)
L126-133: How are bedrock changes due to ice sheet loading? For example, does GSM include Glacial Isostatic Adjustment component?
L136-140: The equations and parameters of the Surface mass balance scheme should be clarified because they are essential for replicating the results. I understand the model is still in development, but the equations and parameters used in this study should be clarified (e.g. section 2.2 of Choudhury et al., 2020).
L145: Please clarify the value of the LGM sea level used in this study.
L150: Please clarify the areas of "continental scale" when defining spatial mean anomaly of precipitation bias correction.
L152: "we impose further ad-hoc temperature increases in these two regions (ranging from +1K to +9K)" What does this mean? Firstly, please clarify the area of temperature anomaly by using a figure. Secondly, does it mean that the present-day temperature over Alaska and Siberia in the experiments is much warmer than the ERA climatology? The discussion section should discuss how this ad-hoc temperature anomaly affects the results. In addition, it is unclear from the text why it has to be 1 to 9°C range.
Method: What about the equations or parameters for Antarctic ice sheets? For example, the equation for Surface mass balance is the same between the Northern Hemisphere ice sheets and the Antarctic ice sheets. There are some sub-shelf melts in Bahadory and Tarasov (2018), but I could not find a reference article on the evaluation of Antarctic ice sheets in the LCice model. While this article focused on NH, the minimum assessment of Antarctic ice sheet changes is necessary when discussing Antarctic ice sheet changes as in Figures 4 and 5.
L190-192: What does acceleration mean? The transient insolation and GHG changes accelerated by four times? As discussed by Choudhury et al. (2020), depending on the setup of the experiments, the acceleration can break the conservation of freshwater volume.
L202-204: Please clarify the definition of the period of MIS7d, 7c, 5d, 5c, respectively. The same applies to Table 2.
L224-225: According to Table B1, almost all parameters in the "all members" column and "pass all filters" columns are within the uncertainty range. Is it true? I assume "all members" corresponds to 2000 simulations and "pass all filters" corresponds to 14 NROY simulations, but it seems strange that almost all parameters are the same within uncertainty.
L254: What is the method of tuning of sub-shelf ocean temperature bias corrections for Antarctica?
L282, L317, and L323 : missing uncertainty range information; please make it consistent with the Tables.
L445-455: According to Figure A1 (left bottom), LOVECLIM has an excessive precipitation bias over Western Canada. As the precipitation bias correction is defined as the continental-scale spatial mean anomaly (L149-150), this precipitation bias might have contributed to the pre-LGM merging of the northern Laurentide and Cordilleran ice sheets in the simulations, which can be discussed.
L466-472: This section discusses -2, 0, 4℃ isotherm as indicators of ice sheet margin. Is there any objective way to extract these numbers as a representative? In addition, do these numbers correspond to the ablation scheme of the LCice model?
Table 1:
Firstly, please clarify that the areal-mean value is used for seven regions for 2-m temperature, precipitation and ocean temperature.
Secondly, the sign of the change for the filter criteria MIS7c/5c would be the opposite: "decreased ice sheet volume compared to MIS7d" or "increased sea level compared to MIS7d"Figure 4: insolation axis missing
Figure 7: this figure can be moved to SI because it is unrelated to method or result.
Figures 10 and 11: Are the selected NROY ensemble members common between the two figures?
Figure 11: please clarify "selected three NROY ensemble members"
References:
Kageyama, M., Charbit, S., Ritz, C., Khodri, M., and Ramstein, G.: Quantifying ice-sheet feedbacks during the last glacial inception, Geophys. Res. Lett., 31, L24203, doi:10.1029/2004GL021339, 2004.
Gregory, J. M., Browne, O. J. H., Payne, A. J., Ridley, J. K., and Rutt, I. C.: Modelling large-scale ice-sheet–climate interactions following glacial inception, Clim. Past, 8, 1565–1580, https://doi.org/10.5194/cp-8-1565-2012, 2012.
Annan, J. D., Hargreaves, J. C., and Mauritsen, T.: A new global surface temperature reconstruction for the Last Glacial Maximum, Clim. Past, 18, 1883–1896, https://doi.org/10.5194/cp-18-1883-2022, 2022.
Abe-Ouchi, A., Saito, F., Kawamura, K., Raymo, M. E., Okuno, J., Takahashi, K., and Blatter, H.: Insolation-driven 100,000-year glacial cycles and hysteresis of ice-sheet volume, Nature, 500, 190–194, https://doi.org/10.1038/nature12374, 2013.
Ganopolski, A., Calov, R., and Claussen, M.: Simulation of the last glacial cycle with a coupled climate ice-sheet model of intermediate complexity, Clim. Past, 6, 229–244, https://doi.org/10.5194/cp-6-229-2010, 2010.
Choudhury, D., Timmermann, A., Schloesser, F., Heinemann, M., and Pollard, D.: Simulating Marine Isotope Stage 7 with a coupled climate-ice sheet model, Climate of the Past, 16, 2183–2201, https://doi.org/10.5194/cp-16-2183-2020, 2020Citation: https://doi.org/10.5194/egusphere-2025-495-RC2 -
AC2: 'Reply on RC2', Marilena Geng, 28 May 2025
Thank you for the thorough review! Here are our responses to your comments:
Major comments:
1) Concerns about bias correction: I recommend the authors to have a set of simulations with identical bias correction throughout the simulations as a reference simulation.
Fair point! Simulations with constant bias correction are set up and their analysis will be included (simulations should be done in a few weeks).
Specific comments:
L2-3: The sentence might be accurate, but as this study uses sea level changes during MIS7 and MIS5 as constraints of the NROY simulations, this statement can be more specific based on the method of the simulations.
Rephrased: “Here, we present the results of a sea-level constrained perturbed parameter ensemble of transient simulations of the last two glacial inceptions and subsequent interstadials (MIS 7e-7c, 240-215 ka and MIS 5e-5c, 122-98 ka) with the fully coupled ice/climate model LCice.”
L4: It can be clarified that insolation and GHG changes are used as transient forcing in the experiments.
Added: “LCice includes all critical direct feedbacks between climate and ice and uses only greenhouse gas and insolation forcing.”
L6-10: What does "model-centric context" mean in this sentence? According to Table B1, most of the parameters after MIS7 and MIS5 constraints are within range of the original parameters adapted in the original LOVECLIM.
We have removed “in a model-centric context” as it’s superfluous. The core point is how sensitive the Earth's dynamical system (as represented by always imperfect computational models) is to uncertainties in the configuration (as represented by the parameter vector).
L16: please indicate the uncertainty range as in Table 2.
Done: “The North American and Eurasian ice sheets advanced southward with rates well above 100 m/yr (168 ± 35 m/yr and 158 ± 45 m/yr respectively) during the penultimate glacial inception and over 70 (73 ± 84 m/yr, Eurasia) and 90 (97 ± 35 m/yr North America) m/yr during last glacial inception.
L66: Studies on glacial inception using coupled ice sheet-climate models (Kageyama et al., 2004; Gregory et al., 2012) should be referred to.
We added Kageyama 2004 in with the group of models that struggle to simulate EA ice sheet extent and added a sentence to include Gregory 2012: “Uncoupled steady-state experiments (climate- or ice-only) lack key feedbacks and they therefore are not able to capture expected rates of sea level change (Khodri et al., 2001; Yoshimori et al., 2002; Vettoretti and Peltier, 2003; Otieno and Bromwich, 2009; Born et al., 2010; Colleoni et al., 2014). Fully-coupled simulations using constant forcing cannot explain the transition in and out of glacial states (Gregory et al., 2012). Early transient simulations of the entire last glacial cycle used two-dimensional energy balance or quasi-geostrophic climate models coupled to simple two or three-dimensional ice sheet models (Gallée et al., 1992; Peltier and Marshall, 1995; Tarasov and Peltier, 1997). While these models could capture the overall structure of ice growth and decay, they were unable to capture the minimum ice volume required to explain proxy-based inferences for the last glacial inception sea level low-stand. More recent work has employed Earth system models of intermediate complexity (EMICs) coupled to ice sheet models. Such studies to date that tested inception and subsequent ice retreat did not include an interactive Antarctic ice sheet. Most such studies also have significant discrepancies between simulated and geologically-inferred ice extent. For example, Alaska tends to be nearly fully ice-covered (e.g. Bahadory et al., 2021; Bonelli et al., 2009; Willeit et al., 2023) contradicting geological records (Kauman and Manley, 2004; Kaufman et al., 2011). Meanwhile, there is often not enough ice simulated over Quebec and Eurasia (EA) (e.g. Bonelli et al., 2009; Ganopolski et al., 2010; Kageyama et al., 2004, at least for LGM where geological data is available for comparison).
L 96: As the model name is "LCice2.0," please clarify the major differences from LCice version 1.0 (Bahadory Tarasov 2018; 2021).
Revised text: “Since the first publication using LCice (Bahadory and Tarasov, 2018, LCice1.0), we have updated the model to Lcice2.0. LCice2.0 now includes a dynamic Antarctic ice sheet, updates were made to the GSM, and a climate bias correction and a radiative factor were added to the parameter vectors.”
L98: which insolation data is used?
Added citation: Berger, 1978
L99: Bereiter et al.(2015) include only CO2 time-series. What about other greenhouse gas components? For example, Choudhury et al., (2020) used three GHG components (CO2, CH4, N2O) from Dome C ice core (Loulergue et al., 2008; Lüthi et al., 2008; Schilt et al., 2010)
We also use a CH4 time series. The citation is added now: “LCice is forced by orbital parameters (Berger, 1978) and greenhouse gas chronologies (Bereiter et al., 2015; Loulergue et al., 2008, CO2 and CH4 respectively).”
L126-133: How are bedrock changes due to ice sheet loading? For example, does GSM include Glacial Isostatic Adjustment component?
We’ve adjusted model structure figure (Figure X, relative sea level solver -> visco-elastic GIA) and added the following sentence to fill in some gaps: “The GSM includes full thermomechanical coupling, global visco-elastic glacial isostatic adjustment with a 0-order approximation to compute spatial variations in the Geoid, a pro-glacial lake resolving surface drainage solver, and both marine and lacustrine ice calving.”
L136-140: The equations and parameters of the Surface mass balance scheme should be clarified because they are essential for replicating the results. I understand the model is still in development, but the equations and parameters used in this study should be clarified (e.g. section 2.2 of Choudhury et al., 2020).
The cited “in preparation” was incorrect, as the GSM description paper had already been under GMD discussion in Jan/25 (Tarasov, L., Lecavalier, B. S., Hank, K., and Pollard, D.: The glacial systems model (GSM) Version 24G, Geosci. Model Dev. Discuss. [preprint], https://doi.org/10.5194/gmd-2024-175, in review, 2025.). The discussion version covers most of the surface mass balance scheme details, but the revised version is more complete and now posted on Tarasov’s research webpage. As such, we see no value to duplicate the description in our current submission.
L145: Please clarify the value of the LGM sea level used in this study.
The nominal scalar value for this index is 120 m, now clarified in text: “At LGM sea level (nominal 120 m ESL), the correction is 0, and at PD sea level, the full correction is applied.”
L150: Please clarify the areas of "continental scale" when defining spatial mean anomaly of precipitation bias correction.
Not sure how to make this more clear, but perhaps “scale” is the problem, so we have just removed that word (i.e. the default is a single separate bias correction for each separate ice sheet).
L152: "we impose further ad-hoc temperature increases in these two regions (ranging from +1K to +9K)" What does this mean? Firstly, please clarify the area of temperature anomaly by using a figure. Secondly, does it mean that the present-day temperature over Alaska and Siberia in the experiments is much warmer than the ERA climatology? The discussion section should discuss how this ad-hoc temperature anomaly affects the results. In addition, it is unclear from the text why it has to be 1 to 9°C range.
We will add supplemental figures of both the ad-hoc warming field and the PD July bias corrections (or actually the negative of the bias correction: LOVECLIM(0ka) - ERA05) for the Northern Hemispheric ice sheets for the NROY LOVECLIM parameter vectors. There are a number of issues here that will be addressed in a revised discussion. For instance, it has long been known that keeping Alaska from fully glaciating over is a challenge for typical paleo ice sheet model resolutions as the deep warm valleys are hardly or not at all resolved (e.g. Marshall, 2002).
We already make clear that this is to suppress excessive ice in those regions. However, an inaccurate statement in section 4.3 has now been corrected from “[..] simulation of the expected NA ice sheet geometry, at least in part due to the imposition of ice sheet wide present-day bias corrections along with glacial temperature bias corrections for Alaska and Siberia.“ to: “[...] simulation of the expected NA ice sheet geometry, at least in part due to the imposition of ice sheet wide present-day bias corrections along with the imposed enhanced warming for Alaska.”
Method: What about the equations or parameters for Antarctic ice sheets? For example, the equation for Surface mass balance is the same between the Northern Hemisphere ice sheets and the Antarctic ice sheets. There are some sub-shelf melts in Bahadory and Tarasov (2018), but I could not find a reference article on the evaluation of Antarctic ice sheets in the LCice model. While this article focused on NH, the minimum assessment of Antarctic ice sheet changes is necessary when discussing Antarctic ice sheet changes as in Figures 4 and 5.
The surface mass balance, submarine melt, and calving models are all fully described in the GSM submission under open review. These same models are applied to all ice sheets. We will be adding some discussion about Antarctica and have submitted the NROY parameter vectors without dynamic Antarctic ice sheet coupling for a direct comparision. The main take-away is that this is going to be a long-term challenge for fully coupled ice and climate modelling given the high sensitivity of grounding line positions to errors in the ocean temperature forcing (and the associated uncertainties in the submarine melt model). As a case in point, no one else has yet published the results of such a fully coupled model.
L190-192: What does acceleration mean? The transient insolation and GHG changes accelerated by four times? As discussed by Choudhury et al. (2020), depending on the setup of the experiments, the acceleration can break the conservation of freshwater volume.
We added a clarification to the end of section 2.2: “In this context, "acceleration" refers to asynchronous coupling between the climate and ice sheet components: under a 4 x acceleration, five years pass in the climate model while twenty years pass in the ice sheet model. The orbital and greenhouse gas forcings applied to the climate are effectively condensed, i.e., 20 years of forcing are applied over just 5 years of climate simulation. Freshwater volumes are not conserved as resulting fluxes would likely trigger excessive response of the meridional overturning circulation. Instead, depending on the ensemble parameter vector, freshwater flux from the GSM is either conserved or enhanced by a factor of up to 2 (in which case freshwater volume is still reduced by 50% given the 4X acceleration of LOVECLIM). Furthermore, this enhancement is only applied to the signed freshwater flux due to changes in ice mass. Given challenges with LOVECLIM internal treatment of freshwater routing, the internal routing within LOVECLIM of precipitation minus evaporation is retained. As such, this component of the flux is conserved."
L202-204: Please clarify the definition of the period of MIS7d, 7c, 5d, 5c, respectively. The same applies to Table 2.
Added filtering periods in the table and in the text: “Ensemble members are ruled out if their simulated sea level minima does not reach at least 25 m at MIS 7d (filtering interval 235-217 ka) and 20 m at MIS 5d (filtering interval 115-105 ka). Furthermore, ensemble members need to display an increased sea-level from MIS 7d to MIS 7c (filtering interval 217-205 ka) and from MIS 5d to MIS 5c (filtering interval 105-95).”
L224-225: According to Table B1, almost all parameters in the "all members" column and "pass all filters" columns are within the uncertainty range. Is it true? I assume "all members" corresponds to 2000 simulations and "pass all filters" corresponds to 14 NROY simulations, but it seems strange that almost all parameters are the same within uncertainty.
The given ranges are not uncertainty ranges but ensemble statistics. Given the large number of varied parameters and the non-linear nature of the model/climate system, we find the behaviour not strange. There is no clearly favoured parameter value for any parameter as a high value for one parameter can be offset by a low value for another.
L254: What is the method of tuning of sub-shelf ocean temperature bias corrections for Antarctica?
An adhoc combination of tuning for present-day mean inferred submarine melt for each major ice shelf along with results from a recent glacial cycle history matching for the AIS (Lecavalier and Tarasov, 2025).
L282, L317, and L323 : missing uncertainty range information; please make it consistent with the Tables.
Done.
L445-455: According to Figure A1 (left bottom), LOVECLIM has an excessive precipitation bias over Western Canada. As the precipitation bias correction is defined as the continental-scale spatial mean anomaly (L149-150), this precipitation bias might have contributed to the pre-LGM merging of the northern Laurentide and Cordilleran ice sheets in the simulations, which can be discussed.
Fair point! We will add a discussion of the implications of the present day precipitation bias.
L466-472: This section discusses -2, 0, 4℃ isotherm as indicators of ice sheet margin. Is there any objective way to extract these numbers as a representative? In addition, do these numbers correspond to the ablation scheme of the LCice model?
We do not understand how a temperature can correspond to an ablation scheme? The ice sheet margin is from the GSM output of Lcice which necessarily uses the GSM’s non-trivial surface mass-balance scheme (see GSM paper, Tarasov et al., 2025) which depends on both positive degree days as well as surface insolation when air temperature is greater than or equal to 0C. As for an objective extraction of these numbers, this would be non-trivial and the visual analysis would still be needed to comprehend the regional variations. Nevertheless, we will examine whether the software tools we have can address this.
Table 1: Firstly, please clarify that the areal-mean value is used for seven regions for 2-m temperature, precipitation and ocean temperature.
Changed “mean” to “spatio-temporal mean” to clarify
Table 1: Secondly, the sign of the change for the filter criteria MIS7c/5c would be the opposite: "decreased ice sheet volume compared to MIS7d" or "increased sea level compared to MIS7d"
Fixed.
Figure 4: insolation axis missing
Added.
Figure 7: this figure can be moved to SI because it is unrelated to method or result.
Done.
Figures 10 and 11: Are the selected NROY ensemble members common between the two figures?
Not intentionally. Ensemble members in Figure 10 (now 9) were picked based on having the same ice volume at MIS 5d. Ensemble members in Figure 11 (now 10) were independently from that picked based on having a large difference between their MIS 7d and MIS 5d maximum extent, as well as having a distinct extent from each other.
Figure 11: please clarify "selected three NROY ensemble members"
Done.
Citations:
Berger, A.: Long-Term Variations of Daily Insolation and Quaternary Climatic Changes, Journal of Atmospheric Sciences, 1978
Loulerluge, L.; Schilt, A.; Spahni, R.; Masson-Delmotte, V.; Blunier T.; Lemieux, B.; Barnola, J.; Raynaud, D.; Stocker, T.; Chappellaz, J.: Orbital and millennial-scale features of atmospheric CH4 over the past 800,000 years, Nature, 2008
Marshall, S. J.: Modelled nucleation centres of the Pleistocene ice sheets from an ice sheet model with subgrid topographic and glaciologic parameterizations, Quaternary International, 2002, https://doi.org/10.1016/S1040-6182(02)00033-2
Lecavalier, Benoit S.; Tarasov, Lev: A history-matching analysis of the Antarctic Ice Sheet since the Last Interglacial – Part 1: Ice sheet evolution. The Cryosphere, 2025, 10.5194/tc-19-919-2025
Citation: https://doi.org/10.5194/egusphere-2025-495-AC2
-
AC2: 'Reply on RC2', Marilena Geng, 28 May 2025
-
RC3: 'Comment on egusphere-2025-495', Anonymous Referee #3, 03 Apr 2025
Summary
Geng and colleagues have conducted PPE simulations using a coupled ice sheet–climate model to study glaciation and subsequent ice-sheet variations during MIS 7 and MIS 5. In these simulations, 18 uncertain climate model parameters were varied and constrained based on their ability to reproduce the seasonality of the modern climate. Using the constrained members, the authors simulated two glaciation events and found that 14 of them satisfied sea-level constraints. These 14 members did not exhibit a strong preference for specific parameter spaces. The simulated ice-sheet configurations in Eurasia (EA) and North America (NA) from the best-performing members were then compared against existing paleoenvironmental records. Lastly, the authors showed that the relative sizes of individual ice sheets can vary among members, even when the total ice volume is similar.
General comments
I find the simulation results presented here interesting and relevant to the readers of Climate of the Past. However, I believe there is room for improvement in the analysis of the modeling results, as outlined below. Therefore, I recommend major revisions before the manuscript can be accepted for publication.
My main comment focuses on analyzing the causes of the spread in the ensemble simulations. While this study emphasizes paleo data–model comparisons, the authors acknowledge that proxy data are sparse. Given this limitation, I believe it would be more valuable to investigate the reasons behind the model spread in the ensemble simulations, providing readers with a clearer physical understanding of the coupled ice sheet–climate system's dynamics.
For instance, it would be particularly interesting to explore why some ensemble members successfully simulate glacial inception while others do not, despite similarly representing the modern climate. In Section 3.1, the authors note the importance of nonlinear interactions among parameters, but I would appreciate further elaboration. Previous studies (e.g., Gregory et al. 2012; Sagoo et al. 2021) have highlighted the critical roles of ice albedo and clouds in glacial inception and others (e.g., Gandy et al. 2023; Sherriff-Tadano et al. 2024) have emphasized their influence on ice sheet maintenance. In the context of these PPE simulations, could these effects be canceling each other out? Or are other nonlinear compensatory processes at play? A pair plot, similar to those presented in Gandy et al. (2023, Fig. 4) and Quiquet et al. (2018, Fig. 9), might help clarify these relationships. Additionally, is there a link between the magnitude of ice sheet inception and global cooling in the model? This aspect is often difficult to investigate using more complex Earth system models, making it a potential advantage of the modeling approach used here.
Secondly, this is a minor comment, but I would appreciate more explanation regarding the motivation for comparing the two glacial inceptions. As far as I understand, the timing of changes in obliquity and precession differs between them, which could lead to different sensitivities of inception to uncertain parameters. Adding an explanation in the paragraph around L74–79 might help clarify this point, though I leave it to the authors to decide where best to include it.
Specific comments
L20-22: I suggest revising the concluding sentence to better highlight the key findings of this study or suggest potential directions for future research.
L59: The issue of not including Antarctica was unclear at first. While this is discussed later in Section 3.4, a brief explanation here would help clarify it.
L79-80: Agreed!
L116-117: Please provide the exact value here for clarity.
L134: Could you specify the depth? Does it vary by location?
L144: What is the rationale for this choice? Please add a brief explanation.
L159 &L225-226: It would be helpful to reference Appendix B1 here. Additionally, how did the authors determine the initial range of parameters in the ensemble simulations? If there is a relevant reference, please cite it; otherwise, consider adding an explanation in Appendix B1.
L165: What about the actual summer temperature, which is crucial for ice sheet development?
L212: Is this due to the weaker increase in summer insolation (Figs. 4 and 5)?
L215-217: To avoid confusion, it would be better to ensure consistency with the values shown in Table 2.
L224-225: Why is the greenhouse gas radiative factor slightly constrained? Could this be related to the magnitude of global cooling during the inception?
L232-234: This finding is consistent with Abe-Ouchi et al. (2013, Fig. 2) and may be linked to warmer conditions over Eurasia compared to North America. Please consider adding a sentence or discussion on this point where appropriate.
L302-308: Repetition?
L308: “MIS7d inception”?
L318-319: Is there any relation between the speed of ice sheet advance and the model parameters used for the ice sheet?
L337-338 & L456-464: Given the large uncertainties on both the proxy and modeling sides, I would prefer to see a more detailed analysis of the causes of the model spread, especially since this is a modeling study. For instance, Fig. 9 shows one member without the merger of the two ice sheets. What explains this? Is it related to low albedo parameters, weak global cooling, or another factor?
L496-498: What is the underlying reason for this result?
Figs. 4 & 5: Please add a label for the insolation.
References
Abe-Ouchi, A., Saito, F., Kawamura, K., Raymo, M. E., Okuno, J., Takahashi, K., & Blatter, H. (2013). Insolation-driven 100,000-year glacial cycles and hysteresis of ice-sheet volume. Nature, 500(7461), 190–193. https://doi.org/10.1038/nature12374
Gandy, N., Astfalck, L. C., Gregoire, L. J., Ivanovic, R. F., Patterson, V. L., Sherriff-Tadano, S., Smith, R. S., Williamson, D., & Rigby, R. (2023). De-Tuning Albedo Parameters in a Coupled Climate–Ice Sheet Model to Simulate the North American Ice Sheet at the Last Glacial Maximum. Journal of Geophysical Research: Earth Surface, 128(8), e2023JF007250. https://doi.org/10.1029/2023JF007250
Gregory, J. M., Browne, O. J. H., Payne, A. J., Ridley, J. K., & Rutt, I. C. (2012). Modelling large-scale ice-sheet–climate interactions following glacial inception. Climate of the Past, 8(5), 1565–1580. https://doi.org/10.5194/cp-8-1565-2012
Quiquet, A., Dumas, C., Ritz, C., Peyaud, V., & Roche, D. M. (2018). The GRISLI ice sheet model (version 2.0): calibration and validation for multi-millennial changes of the Antarctic ice sheet. Geoscientific Model Development, 11(12), 5003–5025. https://doi.org/10.5194/gmd-11-5003-2018
Sagoo, N., Storelvmo, T., Hahn, L., Tan, I., Danco, J., Raney, B., & Broccoli, A. J. (2021). Observationally Constrained Cloud Phase Unmasks Orbitally Driven Climate Feedbacks. Geophysical Research Letters, 48(6), e2020GL091873. https://doi.org/10.1029/2020GL091873
Sherriff-Tadano, S., Ivanovic, R., Gregoire, L., Lang, C., Gandy, N., Gregory, J., Edwards, T. L., Pollard, O., & Smith, R. S. (2024). Large-ensemble simulations of the North American and Greenland ice sheets at the Last Glacial Maximum with a coupled atmospheric general circulation–ice sheet model. Climate of the Past, 20(7), 1489–1509. https://doi.org/10.5194/cp-20-1489-2024
Citation: https://doi.org/10.5194/egusphere-2025-495-RC3 -
AC1: 'Reply on RC3', Marilena Geng, 28 May 2025
Thank you for lots of good comments and ideas how to improve this submission! Here are our answers:
General comments:
1)It would be particularly interesting to explore why some ensemble members successfully simulate glacial inception while others do not, despite similarly representing the modern climate. In Section 3.1, the authors note the importance of nonlinear interactions among parameters, but I would appreciate further elaboration. Previous studies (e.g., Gregory et al. 2012; Sagoo et al. 2021) have highlighted the critical roles of ice albedo and clouds in glacial inception and others (e.g., Gandy et al. 2023; Sherriff-Tadano et al. 2024) have emphasized their influence on ice sheet maintenance. In the context of these PPE simulations, could these effects be canceling each other out? Or are other nonlinear compensatory processes at play? A pair plot, similar to those presented in Gandy et al. (2023, Fig. 4) and Quiquet et al. (2018, Fig. 9), might help clarify these relationships. Additionally, is there a link between the magnitude of ice sheet inception and global cooling in the model? This aspect is often difficult to investigate using more complex Earth system models, making it a potential advantage of the modeling approach used here.
We will add a subsection to the revised submission analyzing any possible differences in ensemble members passing the filtering and those that don’t. Preliminary results indicate that runs that do not pass filtering have both little NA and EA ice volume (i.e. limited EA ice growth isn’t the only key issue), and that there is no significant difference in AMOC strength between good and bad runs. Given the quasi-geostrophic atmosphere, we are more hesitant to consider global statistics that will be heavily influenced by equatorial region response in the model. As we will make clearer in the text, the purpose of this submission is not to decipher mechanistic controls and feedback loops. That is much more clearly done with appropriate sensitivity experiments which will be presented in an upcoming submission. Aside from the GHG radiative factor, the large dimension of the ensemble parameter space, parameter interactions, and nonlinearity means that even pair plots are unlikely to tease out simple relationships, as born out by preliminary checks. We will do a more detailed examination before finalizing revisions, but parametric controls are much better addressed with a different experimental design beyond the bounds of this submission.
2) I would appreciate more explanation regarding the motivation for comparing the two glacial inceptions. As far as I understand, the timing of changes in obliquity and precession differs between them, which could lead to different sensitivities of inception to uncertain parameters. Adding an explanation in the paragraph around L74–79 might help clarify this point.
The inceptions were chosen as they are the two most recent periods of rapid ice evolution starting with similar-to-present ice configuration (to avoid high uncertainties in initial conditions). While data is limited, at least there is some data available to compare model results to (rather then going even further back in time). The introduction will be revised to make clearer motivation and research questions of this paper, including the reason for this comparison.
Specific comments:
L20-22: I suggest revising the concluding sentence to better highlight the key findings of this study or suggest potential directions for future research.
Will do.
L59: The issue of not including Antarctica was unclear at first. While this is discussed later in Section 3.4, a brief explanation here would help clarify it.
Agreed. And as the role of Antarctica has been brought up by all referees, we have submitted the NROY ensemble without dynamic Antarctic ice sheet coupling to analyse the role of Southern Hemisphere freshwater on Northern Hemisphere ice sheets. Antarctica will get more attention in the revised submission.
L79-80: Agreed!
:)
L116-117: Please provide the exact value here for clarity.
We base this statement on our reference (Goossee et al. 2010) here, and they don’t give one global value. What we are referencing is: “After 4000 years, the climate reached a quasi-equilibrium state characterized by a huge cooling of more than 25 ◦C over the Laurentide and Fennoscandian ice sheets [...]. The model also simulates a large cooling in the Southern Ocean associated with a large increase in the sea-ice extent there. The cooling is larger over the Atlantic than over the Pacific, in particular northward of 45◦ N. In the tropics, the signal is weaker. In some regions, such as North Australia, the changes are very close to zero. Those results are similar to the ones of other simulations performed in the framework of the PMIP2 project (Braconnot et al., 2007), except in the Southern Ocean where the signal obtained in LOVECLIM is larger than the one given by most other models.”
L134: Could you specify the depth? Does it vary by location?
The profile is taken to a depth of 594 m. The temperature at this depth is used for anything deeper. Based on ref1’s comments we have also added: “To facilitate dealing with complicated ocean grids as well as enable ocean temperature forcing where the land mask has a terrestrial value, ocean temperature profiles from CLIO for computing submarine melt in the GSM are extracted at key index sites and then applied to an assigned downstream region (Bahadory and Tarasov 2018).”
L159 &L225-226: It would be helpful to reference Appendix B1 here. Additionally, how did the authors determine the initial range of parameters in the ensemble simulations? If there is a relevant reference, please cite it; otherwise, consider adding an explanation in Appendix B1.
Agreed! We added: “We use fit to present-day (PD) climate to select initial LCice parameter vectors for glacial inception simulations. To do so, plausible prior distributions for 18 LOVECLIM parameter values were defined. The ranges of parameter values were inspired by Shi et al. (2019) (see Supplement Table 1)”
L165: What about the actual summer temperature, which is crucial for ice sheet development?
Even state-of-the-art GCM climate models have significant biases. Given the bias corrections invoked, capture of present-day summer temperature is not critical. As stated in the text, with bias corrections, the amplitude of the seasonal cycle is reasonable present-day metric for sensitivity to radiative forcing (though with the limitation of different time scale, which will be most relevant to deep ocean constraint) is a more critical metric to capture for selecting parameter vectors more likely to give an adequate glacial cycle response.
L212: Is this due to the weaker increase in summer insolation (Figs. 4 and 5)?
We’ve added: “MIS 5d sea level poses the strongest constraint on the ensemble. The majority of ensemble members cannot grow enough ice, likely due to the relatively weaker CO2 decrease during inception.”
215-217: To avoid confusion, it would be better to ensure consistency with the values shown in Table 2.
Fixed.
L224-225: Why is the greenhouse gas radiative factor slightly constrained? Could this be related to the magnitude of global cooling during the inception?
We’re not clear on the intent of this question. Other LOVECLIM ensemble parameters will also affect radiative forcing sensitivity, but apparently to an insufficient extent as required for the capture of glacial inception. Our results show that even given the ensemble parameter exploration we’ve carried out, LOVECLIM needs greater than 1.5 enhancement of the GHG radiative forcing.
L232-234: This finding is consistent with Abe-Ouchi et al. (2013, Fig. 2) and may be linked to warmer conditions over Eurasia compared to North America. Please consider adding a sentence or discussion on this point where appropriate.
Will do
L302-308: Repetition?
Yes, thanks! Fixed.
L308: “MIS7d inception”?
Fixed
L318-319: Is there any relation between the speed of ice sheet advance and the model parameters used for the ice sheet?
No. It’s the combination of the many parameters and cannot be disentangled here. This will depend on horizontal temperature gradients, upstream precipitation fluxes, changes in the seasonality of precipitation,...
L337-338 & L456-464: Given the large uncertainties on both the proxy and modeling sides, I would prefer to see a more detailed analysis of the causes of the model spread, especially since this is a modeling study. For instance, Fig. 9 shows one member without the merger of the two ice sheets. What explains this? Is it related to low albedo parameters, weak global cooling, or another factor?
There is no one parameter that can explain the behaviour. As mentioned in response to the general comment, we will add a section dissecting any characteristics that differentiate “good” from “bad” ensemble members. However, a clear answer to the question “what makes ice grow where” is beyond the scope of this study. It will rather be addressed in a planned sensitivity analysis where individual parameters/feedbacks can be turned on/off.
L496-498: What is the underlying reason for this result?
A good question that will be elaborated upon in the revisions. There are a few facets involved. We have not ruled out that a more complete set of criteria for what constitutes an adequate simulation of present-day climate might improve the chances of capturing glacial inception. There is also the issue of possible equifinality in climate models vis a vis capture of present-day climate not holding for glacial inception capture. The two-way coupling of an ice sheet and climate model adds extra feedback loops into the climate system, which will change the full system sensitivity on longer timescales. Finally, any tuning of a complex model with many parameters will entail trade-offs; the more metrics you have, the harder it is to optimize each one. Note that all the facets are related/overlapping to varying degrees.
Figs. 4 & 5: Please add a label for the insolation.
Done.
Citation: https://doi.org/10.5194/egusphere-2025-495-AC1
-
AC1: 'Reply on RC3', Marilena Geng, 28 May 2025
Viewed
| HTML | XML | Total | BibTeX | EndNote | |
|---|---|---|---|---|---|
| 228 | 49 | 15 | 292 | 14 | 17 |
- HTML: 228
- PDF: 49
- XML: 15
- Total: 292
- BibTeX: 14
- EndNote: 17
Viewed (geographical distribution)
| Country | # | Views | % |
|---|
| Total: | 0 |
| HTML: | 0 |
| PDF: | 0 |
| XML: | 0 |
- 1