Revisions of “GREB-ISM v0.3: A coupled ice sheet model for the Global Resolved Energy Balance model for global simulations on time-scales of 100 kyr”

Dear Editor and referees, We would like to thank the two referees and the editor for their time spent on reviewing this manuscript, and for the many very helpful comments they provided. We think the referee comments have helped us to substantially improve the presentation of this work. For the new manuscript and model, we adopted the suggestion from Referee #2 and changed our standalone transition experiment boundary condition, which has greatly improved this sensitivity experiment. Meanwhile, we fixed several minor bugs in the model code that did have a minor effect on the coupled model simulations. Below we give a point-to-point response to all referee comments, hoping the revised manuscript has now been improved in clarity and is ready for publication.

the authors hint that the coarse vertical resolution might impact the strong differences in basal temperature. The model being very cheap, I would appreciate seeing a sensitivity to rule this out.
Response: Please note, that the aim of this study is to provide a model tool to understand iceage cycle on 100kyrs time-scale, for the global climate and the corresponding physical mechanisms. Our study focusses on illustrating the fully coupled model for this aim. The EISMINT experiments are one element in testing this, but given they only focus on a higherresolution problem for only a small region they are not the most important ones. We therefore have to strike a balance and have to focus more on the global coupled system for which we have presented a number of different sensitivity studies.
Nevertheless, we tried to improve the discussion of the EISMINT experiments. In order to verify our assumptions that the basal temperature bias comes from vertical resolution, we rerun the EISMINT experiments with 10 layers. For EISMINT I (see Table R1), the basal temperature with 10 layers from ISM is much closer to the results in H96 than that with 4 layers. As for EISMINT II (Table R2), increase of vertical resolution leads to a higher percentage melt fraction, closer to results from P2000. These results provide good support for our assumptions.
We further also did some test with the horizontal resolution, which does improve the mass flux at the midpoint. We have now mention these results in the manuscript to give the reader some idea of what is causing the GREB-ISM model limitations in the EISMINT experiments.

Additionally, the horizontal resolution of the ISM is increased in the EISMINT experiments to
better compare the results to the one published (I assume). Typically, good results at coarse resolution translates to as good or better results at higher resolutions. However, you plan on using your ISM on a much coarser global grid than the one used in the EISMINT experiments.
I would have liked to see a comparison at that resolution as well.

Response:
The reviewer is right, we changed the GREB-ISM grid to match with the EISMINT experiment setting. If we use GREB resolution setting, we have only two grid cells in the EISMINT domain, which limits the extent to which a useful comparison can be made between the GREB-ISM and EISMINT simulations. The EISMINT experiments have only a limited scope and application, which is why we also present a series of other sensitivity experiments.
Our dynamic equilibrium and paleoclimate transition experiments are designed to show the model capability to reproduce realistic ice sheets and related dynamics (transports) in the GREB-ISM grid resolution and parametrization. It turns out that our ice thickness simulation matches well with the observation in Greenland and most of Antarctica (e.g.,Figs. 8,9 and 13).
We have improved the discussion of the EISMINT experiments and also the other experiments to illustrate the skill of the new GREB-ISM in section 4.1 and 4.2. Ultimately, the GREB -ISM real skill can only be evaluated if it is applied to real research problems, while here in the introduction paper we can only give some indication.
There is a lack of discussion about the consequences of using such a coarse horizontal resolution for ice sheet modeling while the literature has shown abundantly that a grid resolution of at least 4km if not higher is necessary to simulate ice sheets (Pattyn et al. (2012), Pattyn et al. (2013), Gladstone et al. (2010), Leguy et al. (2014), Seroussi et al. (2015), Seroussi et al. (2018), Cornford et al. (2020, Leguy et al. (2021), …). This is especially important to simulate marine ice sheets like the West Antarctic Ice Sheet for which GREB-ISM show large differences compared to observation.

Response:
The reviewer points out a number of studies that aim at regional scales like hundreds of kilometers. The reviewer also points at "simulate marine ice sheets like the West Antarctic Ice Sheet". All of these are important aspect of ice sheet modelling, but they are not the primary aim of the GREB-ISM, as we have mentioned in the title, the first sentence of the abstract, the introduction and the summary section: We aim at simulating the time evolution of global-scale ice sheets on time scales of 100 kyr. For these kinds of problem, one has to make compromises to achieve feasible speed on the computing system. Past studies that addressed similar kinds of time scales to the ones we are interested in, also use coarser resolution models (Abe-Ouchi et al., 2007;Ganopolski et al., 2010;Willeit and Ganopolski, 2018). And as a result, the high resolution such as 4 km is not yet feasible on these timescales (Fyke et al., 2011).
We do provide some discussion on the coarse resolution in respect to the simulation of West Antarctica and in more general. We hope the revised text does better highlight this. Please see Sections 4.3 and 6.
Many aspects of the ISM and choices are not clearly stated. For example, you choose to use the basal sliding velocity from Greve 1997. I am not saying you should not, I am simply asking to justify your choice while so many sliding laws are now available and studied (Weertman (1957), Schoof (2005), Schoof (2007), Aschwanden et al. (2013), Leguy et al. (2014, 2021), Tsai (2015, …).

Response:
We have improved the overall discussion of the model development and the sensitivity experiments (see section 4.3). We hope the revised version give a more complete presentation of the model. The sliding law from Greve (1997) is one of Weertman type power law (Weertman, 1957).This type of explicit form sliding law is efficient and widely used in paleoclimate research (Abe-Ouchi et al., 2007;Fyke et al., 2011). In respect to the basal sliding law we added some discussion in Section 4.3 in which we argue that the variations of the basal sliding parameters have no significant impact on the model simulations in the context of our applications.
To test the model sensitivity to sliding law, we run dynamic equilibrium experiments by varying sliding law coefficient (change Csl in Eq. (22)) from 610 3 to 610 5 . Compared with GREB-ISM simulation, only when Csl is larger than 310 5 a -1 , we can see a significant ice thickness reduction on Antarctica ( Fig R1). Increasing sliding coefficient does help reduce the WAIS extra ice growth, but the East Antarctic Ice Sheet loses too much ice at the same time, which worsens slight negative bias in East Antarctica (Fig 13). So our current sliding law coefficient is still a reasonable choice. Hence, our experiments here suggest that stronger ice sheet sliding is unlikely a solution for WAIS positive ice thickness bias.
…Also, you mention in many places in your paper how well simulated calving matches observation, while you never clearly define how you actually calve the ice (unless I missed it in which case I apologize). In section 3.2.2, you mention that a condition for floating ice shelves is that H>=10m. Does it mean that 10m is a thickness threshold below which ice calves?

Response:
We had previously defined calving at the end of Section 3.2.3 (Mass balance). To make this clearer, we now define a subsection 3.2.3.1 (Calving) to better highlight how we diagnose calving.
There is also no discussions of any sub-grid scale parameterizations (for grounding lines, calving fronts, …), does the reader need to understand they aren't any?
I do note that the authors left a more in-depth study of the ISM properties for another publication, but if so, I almost feel like it should have come before publishing this paper to rule out many possibilities leading to large differences in the results.
Response: If we have specific sub-grid parameterization scheme in our model, we will explicitly explain it in the text, such as snowfall rate and sliding law. So, we do not include those sub-grid parameterizations you mentioned yet. The grounding line in our model is decided by the glacier type mask, which is defined in section 3.2.2. We have now added a sentence (section 3.2.2) to explicitly point out that grounding lines of glaciers are defined by shifting points from grounded ice to floating ice and therefore change the method of how the dynamics of the ice flow is simulated. For the simulation of calving, please see our response to the previous comment.
We have also tried to improve our presentation in many other aspects to better illustrate how Response: We rephrased this section to better highlight the long and large-scale ice sheets simulations we are interested in.

P17, l506: I suggest removing this sentence especially since the next sentence refers to results
showing all 3 oscillation periods.

Response:
We think the sentences need to stay, as the following paragraphs indeed focus on the 20kyrs period. This is then followed by paragraphs explaining how things are changing for the longer periods.
P18, l536: this might be a resolution effect even though these places have been covered by ice sheets at some times.
Response: It is interesting that these shallow ocean regions are ice sheet covered in our idealized experiments and as the reviewer points out this may indeed have happened in the past. We made no changes in response to this comment. P19, l566: replace"ISM-GREB" by "GREB-ISM" to be consistent.
Response: It has been corrected.

Please, distinguish in your figure what is sea ice to what is ice sheet.
Response: Please, note Fig. 3 should not be compared to Fig. 1, as Fig 3. focuses on seasonal snow and sea ice cover, but Fig. 1 is focusing on large ice sheets. We did revise this figure to improve the color scheme.

In particular, add what the black box is (it is cumbersome to go back and forth between figure and text to piece all the information together). Same remark as in Fig1 for the color scale.
Response: Panel d is at the end of year 1. This is now mentioned in the figure caption. The color scheme in this picture is mostly distinguishable so we do small modification.

Modify the color scale for panel (a).
Response: The figure has been revised. The Tsurf is limited to the black box is because cooling due to lifting does not lead to strong heat transportation to the surrounding. The potential air temperature, which is used to calculate the heat transport, in lifted region does not significantly change even though the surface temperature change. And air column at the point becomes shorter and thus it is difficult to be transported. Both effects reduce the potential heat transport because of lifting.

Fig6: Define what "R" is in the caption.
Response: It has been added.      Response: We did revise the figure. Please note we also changed the way we present the ice sheet extent. In the previous version we only showed ice covered regions with ice sheets larger than 10m to focus on the ice sheets and not the snow/sea ice cover. However, this was not mentioned in the text and it did not fit well with the discussion. We now show the full ice cover including snow and sea ice cover.

Fig17: increase the font size of color bar labels.
Response: We made these revisions.

Referee #2
The authors introduce a newly developed global ice sheet-climate coupled model which aims to explore the evolutions and variabilities of ice sheet-climate system on very long time scale. Response: We like to thank the referee for the evaluation of our manuscript and the many comments that will help us to improve the manuscript. Our one-by-one response is listed below.
Major comments: 1. My main concern appears in the ambiguity on what kind of insights could the model provide in future studies. This model is clearly not meant to reproduce the ice sheet and climate precisely, but rather meant to reproduce them crudely so that the authors can discuss their interactions on global and long time scale. However, there is still a big problem in the spatial distribution of the simulated ice sheet; thickest glacial ice sheet exists over Siberia (Figs. 10 and 16). Previous studies have shown that the longitudinal locations of ice sheets affect their dynamical characteristics, which is important for interpreting the evolution of glacial-interglacial cycle (e.g. Abe-Ouchi et al. 2013, nature, Fig.2b (Clark et al., 2009;Fyke et al., 2011;Ganopolski et al., 2010;Niu et al., 2019;Velichko et al., 1997;Willeit and Ganopolski, 2018) and the Siberian ice sheet disappears. It is able to capture large European (e.g. Fennoscandia) and North American (Laurentide) ice sheets.
The other oscillation simulations do also create a Siberian ice sheet, but they are idealized short wave forcing simulations that should only illustrate the capability of the GREB-ISM to respond to external forcing in a fully coupled climate system. Future applications in more realistic setups are needed to address this interesting problem. In the final section we highlight what kind of insights the new GREB-ISM model can give.

While the ice sheet model simulates the ice shelf component of the ice sheet, I couldn't find
discussions on the effect of basal melting at ice shelf-ocean margin. Perhaps this effect is partly incorporated through sea ice dynamics, but the authors should discuss possible effects of not including ice shelf basal melting in the model, especially over Antarctica. I assume this is partly causing the too thick west Antarctica in modern climate. Also, it might be related to the small sensitivity of Antarctic ice sheet to climate changes (L466-467).

Response:
We have indeed not included an explicit basal melting scheme but did test it. This is now discussed in the manuscript. The ice shelf velocity in our simulation is parametrized with a relatively large value for viscosity, this implicitly includes some other unresolved physical processes like basal melting effects. We have tested several basal melt schemes but the results do not have a fundamental improvement to our simulations. Where Δ is basal mass balance, , are thermal exchange rate (10 -4 ) and model tuning factor (510 -3 ), respectively. All of the other parameters can be found in Table 1. Δ + in the equation is used to represent temperature difference between ocean and ice sheet basal melt point. This temperature difference is assumed to be a constant offset temperature Δ plus a basal melting temperature modification based on ice thickness .