The authors explore tipping (both bifurcation- and rate-induced) in a three-dimensional model of the Greenland ice sheet coupled to a regional energy-moisture balance model. They find that the ice sheet will tip to a much lower ice state depending on the rate of the forcing and the overall level of the forcing. Considering the rate-induced transitions, they find that the time of tipping appears to relate in a non-monotonic way to the forcing. They conjecture the existence of an intermediate (chaotic edge) state or a boundary crisis to be the reason for this unpredictable time of tipping. Oscillatory behaviour of the system is also explained through the response of the Humboldt and Petermann glaciers and is removed to try to isolate the rate-induced tipping.

I have a couple comments regarding the presentation of the results.

1. Initially I was suspect of the chaotic behaviour claims in the paper. At the end of page 7 it is claimed that there is sensitive dependence on initial conditions, however this is not rigorously justified. Moreover, the model response is then continuously referred to as chaotic. Appendix B, however, has a much stronger argument for the conjecture of chaos in the system, including calculations which aim to verify (reject) the existence of a chaotic attractor (chaotic saddle). While I do not think all of this appendix should be moved into the main text, I strongly feel that some critical conclusions do need to be. In the first mention of chaos. I recommend adding a short discussion detailing your evidence and conjecture for the boundary crisis. Some of the more detailed discussions can still be left in the appendix.

2. The authors claim that there is no critical rate for rate-induced tipping. I find this very questionable. If the system does not tip when forcing is constant and then tips when it is not, there is some critical rate (albeit could be very small). Unless that the authors can prove that the system continues to tip as the rate limits to zero, I recommend softening such statements of no critical rates and acknowledging the critical rate could be smaller than the smallest value tested.

3. The forcing in the model is not well described. I imagine it is some kind of regional temperature anomaly based on its units, but I could not find an explicit description or expression. Additionally, the authors mention that the forcing increase is somehow adapted to equilibration time. This would be nice to see written out explicitly.

4. In general section 2 need to be rewritten so that it is much clearer. I found it very hard to follow the introduction of parameters and governing equations. There is a lot of repetition where a parameter is defined mathematically and then introduced later for its physical interpretation. This can be easily fixed with a careful rephrasing.

5. The terms "tipping time" and "tipping point" seem to be used interchangeably in the text. This is a pedantic comment but important because often they mean different things. Choose one term and clearly state how you are defining it.

6. On page 8 you mention the oscillatory behaviour for the first time.  If there is oscillatory behaviour there could be existence of an unstable periodic orbit as a boundary that could potentially explain the irregular transition times which are perhaps more related to crossing this threshold. This also would strengthen your boundary crisis argument as the chaotic attractor can then collide with an unstable periodic orbit at the crisis.

7. Figure 6 appeared to me to have "edge state" or chaotic saddle like dynamics. Again, perhaps a short discussion of the results in the appendix would be useful in the text to clarify why this is not believed to be the case.

This paper investigates the conditions of emergence of tipping points of the Greenland ice sheet (GrIS). The work is based on an ensemble of simulations performed with an ice sheet model (YELMO) coupled to an energy balance-moisture model (REMBO). The runs consist in different warming scenarios with respect to the present-day period with various amplitudes of temperature anomalies and rates of change. The authors explore in particular b-tipping and r-tipping bifurcations and highlight the emergence of oscillations in the Humbold (HIS) and Petermann (PIS) ice streams (northwestern Greenland) which they attribute to chaotic variability leading potentially to a delay of the tipping time of the GrIS (large-scale ice loss). This study suggests that modest forcings may induce complex and unpredictable responses of the ice sheets. It represents an interesting and innovative contribution despite it remains difficult to see how it can be applied to the real world. The objective of the paper is clearly presented. The manuscript presents plausible processes but, unfortunately, not formally demonstrated. What is interesting are the numerous hypotheses/ideas put forward to interpret the results. These provide a new perspective on the behaviour of the climate-ice sheet system despite the simplicity of the atmospheric model. However, the results require a more in-depth analysis to be more convincing and some sections need to be further developed. This is why I recommend the publication of this paper after the authors have addressed a number of revisions listed below.

Major Comments

1. The dispersion of ice volume curves is interpreted as a signature of chaos (and this seems very plausible), but this has not been formally proven. I recommend therefore to provide additional figures showing visual evidence of the chaotic behaviour. These could be a plot of the ice volume trajectories (i.e. Vice(t+dt) vs Vice(t)) and/or a plot showing the dispersion around the mean of the simulated ice volumes (taking all simulations into account) as a function of time.
2. There is evidence for r-tipping demonstrated with the collapse below the bifurcation threshold for certain ramping rates. However, it remains unclear whether a critical ramping rate exists or whether the tipping results from an internal chaotic variability. This issue is discussed later (Section 4.2), but I think it should be emphasized in the comments of Figures 5 and 6. One way to distinguish pure r-tipping from internal chaotic variability is to plot the final ice volume as a function of the ramping rate for a given DT forcing. This would make it possible to determine whether tipping points occur abruptly beyond a certain threshold or, on the contrary, in a chaotic mode.
3. As I am not specialist of chaos dynamics, I found the Appendix B1 too technical. It took me a while to fully understand what you meant. The first part (lines 347-358) sounds like a lesson on non-linear dynamic systems and I regret that the explanations do not refer more to the physical systems under study, despite this is partially rectified in the second part. For example, I suggest to define what are the monostable or the bistable parameter regimes. By the way, some of these explanations could be moved in the main text. This would allow to clarify the link between local fluctuations and the large-scale ice loss.
4. I have also a couple of questions regarding the origin of the oscillations.

i/ How the initial states A, B and C differ? Are the differences due to the vertical temperature profile? sub-glacial topography, sub-glacial hydrology?

ii/ You show that oscillations can be removed by lowering the basal velocities resulting in an increase of the basal friction. This raises questions regarding the sensitivity to the choice of the sliding law and this could be commented in Section 4.1 (see for example Brondex et al. 2017).

iii/ As oscillations are claimed to be due to thermomechanical coupling, a special attention may be given to the initial ice sheet topography. In this study, model initialization consists in a 400 ka spin-up experiment to arrive to an initial state corresponding to the present-day GrIS. Although such approaches are widely used, particularly for paleo experiments, they generally provide a too large GrIS compared to its actual present-day geometry (with the ice volume above flotation estimated at around 2.9 10⁶ km³). In general, the longer the spin-up, the further the final state is moved away from observations. I acknowledge that the difference between the GrIS initial state and the present-day GrIS geometry is of secondary importance within the scope of the paper and does not question the main findings. However, I think it would be interesting to perform additional simulations to investigate the impact of the GrIS geometry on the chaotic fluctuations and add a dedicated comment in Section 4.1.

1. How do you explain why HIS and PIS behave differently? Are these differences related to the geothermal heat flux, bed topography, basal velocity?

L180-181: The fact that the ice thickness variations of PIS are influenced by the oscillations of HIS in the unretreated configurations sounds a bit speculative. Can you provide other arguments in addition to the proximity of the ice streams? In the same way, they report an influence of PIS on HIS in the retreated case (L194-195). What is the coupling mechanism that would cause one of the oscillators to influence the other?

You mention oscillations in the HIS/PIS sector, but you do not specify whether similar oscillations have been reproduced by the model in other marine-terminating glaciers and ice streams. Should these oscillations be specific to the HIS/PIS ice complex only, can you explain why or provide suggestions? Have such oscillations been detected in the real world with satellite observations in other GrIS basins or in Antarctica?

1. In the introduction, you point out the importance of the melt-elevation and temperature-albedo feedbacks. However, the results are never discussed as a function of these feedbacks. Another effect which has a significant influence through the melt-elevation feedback is the global isostatic adjustment (GIA). Owing to the relaxation time (3000 years) in the ELRA model and the time scales considered in this study, GIA should not be ignored in the analysis/discussion of the results.
2. The presentation of the method requires some additional details concerning the REMBO model, its performance and its resolution. Furthermore, I am not sure that all potential readers of the manuscript are familiar with the ITM method, and I believe it would be necessary to provide the equation linking the melt rate to insolation, temperature and albedo. Finally, what downscaling method is used between REMBO and Yelmo?
3. In the abstract, you mention that ice stream oscillations may complicate the anticipation of ice-sheet tipping points. I guess you refer to near future projections. On the other hand, you explain in the conclusion that the consequences of chaotic transients on anthropogenic climate change are phenomenological rather sociologically relevant because of the longer time scales addressed in this study. Both statements seem to be a bit contradictory. My feeling is that this paper demonstrates that chaotic fluctuations may theoretically occur in ice sheets. However, it is difficult to conclude about potential implications on decadal to centennial (or multi-centennial time scales) without explaining how to reconcile both time scales. However, I would be very interested in having your comments on this issue.
4. The limitations of the modelling approach are briefly mentioned in Section 5 and should be developed further developed either in the Conclusion section or in a dedicated section addressing the main sources of uncertainties. I have noted several points (see below) but there may be others of interest.

i/ The use of a simplified atmospheric model raises questions about its ability to represent properly the accumulation (and its variability) of small spatial scale regions (i.e. the HIS/PIS sector) and, therefore, the surface mass balance.

ii/ The absence of coupling with the ocean results undoubtedly in missing processes that are of key relevance for marine terminating glaciers (see for example Catania et al., 2020, Crowton et al., 2018). This is particularly necessary to mention this point as you refer to Alvarez-Solas and Bassis in Appendix A.

iii/ The resolution of the ice sheet model (16 km) prevents from an explicit representation of small-scale features such ice streams. This could be mentioned

I would like to point out that I am not denigrating your modelling approach, as I am aware that studies are conducted using the available numerical tools. However, I believe it is essential for the climate modelling community to clearly explain the framework within which the results are obtained.

Minor Comments

Several figures are quite difficult to examine (insufficiently detailed captions, lines too thin, colours too similar and font too small).

Figures 1, 2, 3, 5, 6, 7, 8: Increase label font size

Figure 1b caption: surface velocity instead of ice sheet extent

Figure 2a does not resemble an “Equilibrium ice volume”. The term “Equilibrium” should be removed

Figures 3, 7, 8: Increase the size of the colour scale

Fig. 6 Choose another set of colours to make more visible the different curves.

Fig. 7: Panels b and c show other marine-terminating glaciers that are not discussed in the manuscript. Could these maps be restricted to the same region as those displayed In Fig.8?

Fig. 8 caption: Which is the time period considered for the temporal average

L18: Provide also more recent references (e.g. Wunderling et al., 2018)

L43: forcings (instead of forcing)

L48: the ability of what?

L121: saddle-node bifurcation could here be related to melt-elevation and temperature-albedo feedbacks.

L131: Why 120 ramping experiments?

L135-136: Do the authors mean that tipping can occur for ramping rates below 10^-5 K/yr. Did they conduct simulations with lower ramping rates?

L139-142: The sentence seems to contradict the previous sentence (unless I missed something?).

L153: I would say “during the first 25-50 ka of simulation time

L155: Refer here to Figure 7.

L184: “In contrast to the PIS” à “In contrast to the unretreated PIS”

L185: “as the HIS” à “as the unretreated HIS”.

L186: The oscillations are not displayed for the ice extent in panels (b) to (i) of Fig. 8. Replace “extent” with “thickness” or “volume”.

L188-190: If the basal velocity increases, the mass loss should be accelerated. This contradicts the statement “This slows the mass loss”.  Replace “slows” with “accelerates”.

L189-190: This results from the thermomechanical coupling and could be reminded here.

L199-202: The sentence is a bit long. Please split in two.

L201: I did not find in Fig. 8 some features showing periods when both retreated HIS and retreated PIS are in a steady flow. Can you be more specific?

L208: lost à loss

L216-219: Would the results/conclusion be the same if the forcings were between 1.05 and 1.30 K? There is no firm conclusion in section 3.4 concerning the role of the parameter delta.

L243: Please provide references and add comments related to the complexity of models used in other studies. I wonder whether the use of more complex climate and ice-sheet models would enable to reproduce a tipping of the ice sheet on shorter time scales.

L241-242: I suggest that the fact that the occurrence of chaos in this study is an open debate should be mentioned earlier

L246: What is the amplitude of the oscillations in the present work? It is not mentioned in the text? See Major Comments

L247-248: Can you explain why the isostatic adjustment has a much larger role in the Zeitz et al’s (2022) work? Maybe you could also discuss the differences with their results (possibly in terms of the differences in the experimental setups)

L250: Could you add more recent references (e.g., Wood et al., 2021)

L251: “suggesting that…in these areas” should be justified

L253: Note that Joughin et al. (2010) reported more significant changes in southeastern glaciers than in northwestern glaciers due to differences in bed geometries

L282: “the value of the tipping point” à Do you mean here “the value of the forcing”?

L296-298: I am not sure that all potential readers of the manuscript are familiar with the notions of edge-tracking algorithm (like me). Can you define please? Same for ghost attractor and chaotic saddle. Moreover, this is the first time you employ the terms “ghost attractor” and “chaotic saddle”. In a previous comment (see the “Main comments” section), I recommended that part of Appendix B be moved to the main text and that the explanations of these terms refer to the physical system (and not be discussed in excessively theoretical terms).

L365-368: The sentence is too long. This is detrimental for the understanding. Split it in two parts. It would also benefit from a clearer explanation based on terms describing the GrIS and not only on chaos terminology.

L376: I am not sure to understand what you mean with “from the left”?

L378: What do the different parameters actually represent in Eq. B1?

L385: Remove the first occurrence of “either” and the term “otherwise”

L402: What do you mean with “diminished size”?

References

Brondex et al. (2017), Sensitivity of grounding line dynamics to the choice of the friction law, doi: 10.1017/jog.2017.51

Catania G.A et al (2020), Future evolution of Greenland's marine-terminating outlet glaciers. Journal of Geophysical Research: Earth Surface, 125, e2018JF004873. https:// doi.org/10.1029/2018JF004873

Crowton, T.R et al. (), Linear response of east Greenland’s tidewater glaciers to ocean/atmosphere warming, doi: 10.1073/pnas.1801769115.

Wood et al. (2021), Ocean forcing drives glacier retreat in Greenland, https://www.science.org/doi/pdf/10.1126/sciadv.aba7282

Wunderling, et al. (2024), Climate tipping point interactions and cascades: a review, Earth Syst. Dynam., 15, 41–74, https://doi.org/10.5194/esd-15-41-2024, 2024.

The submission by Kypke examine the potentially significant role of

ice stream cycling in the tipping point behaviour of a GRIS model.  To

avoid repetition, my review will largely focus on issues not yet

raised by the other two reviewers.

# major comments

Foremost, I can't adequately assess the significance of this

submission given some missing critical information about the model

configuration, and the lack of any analysis of numerical and input

sensitivities in the model's simulation of ice stream cycling.

For instance, in Hank et al (2023,

https://doi.org/10.5194/gmd-16-5627-2023) we carried out a detailed

analysis of simulated ice stream cycling response to grid resolution

along with an assessment of approaches to minimize that

sensitivity. Given the significant resolution sensitivity identified

in that study ( even between 6.25 and 3.125 km grid resolution), along

with the challenge of accounting for the impact of the large km scale

variation in basal topography around much of the GRIS ice margin, I'm

skeptical of coarse resolution analysis of ice stream cycling dynamics

without associated numerical assessment. At the very least, I would

need to see a couple of ice stream cycling sensitivity experiments at

8 km grid resolution (twice as fine as current). Otherwise, I'm

unclear of the extent to which the current results are just numerical

artefacts.

The submission is also missing critical model configuration

information, in particular that of the 2 km deep geothermal heat flux

boundary condition and the details on the climate forcing (as raised

by others). In that regard, given the strong impact of uncertainties

in deep geothermal flux identified in Hank and Tarasov (2024,

https://doi.org/10.5194/cp-20-2499-2024) for Hudson Strait cycling

(for which one of the co-authors of the current submission,

Alvarez-Solas, was a reviewer), this submission also needs a

sensitivity assessment in response to those uncertainties. This

sensitivity is not surprising, since geothermal heatflux will strongly

impact the time required to reach the basal pressure melting point

required for ice stream activation along with subsequent basal

meltwater production.  The lack of any reference to the above two

papers are also examples of insufficient review of relevant

litterature. Other relevant papers not cited include Soucek and

Martinec (2011, https://doi.org/10.3189/002214311798843278) and Sayag

and Tziperman (2011, https://doi.org/10.1029/2010JF001839).

The precipitation aspect of the climate forcing will likely have the

largest uncertainties. It also plays a major role not just in surface

mass-balance, but also in vertical cold advection to the base with

associated impact on basal temperatures. Given that the    uncertainties

in precipitation are not considered, there at least needs more details

in the appendix  or supplement about exactly what REMBO accounts    for

along with plots of a few precipitation  map timeslices.

Given its pivotal role, the choice of basal drag representation and

associated parameter choices will have a major impact on results. To

justify the base value choice, there needs to be a present-day

comparison of simulated surface velocities against those observed.

The analysis needs more depth as to controls on the stream

activation/quiesence. As an example a few basal temperature maps would

convey the spatial range of basal temperatures proximal to the ice

streams.  For context, I'm a bit surprised with the amount of

lateral velocity diffusion in the simulated ice streams.  Is this due

to most or all of the GRIS bed being close to the pressure melting

point? It would also help to have the map plots made in non-fill (ie

no graphical smoothing) mode to show the actual grid cells.

The core take-away (subject to numerical verification) for me is

further affirmation that ice dynamics (as opposed to just

surface/marine mass-balance processes) can really matter in GRIS

future evolution (though here on a relatively long time-scale

herein). A possible elephant in the room for any detailed analysis of

GRIS ice stream dynamics, especially in the context of chaotic or near

chaotic response, is that even if these simulations were done at 4 km

resolution, they would be at least 4 times too large to even partially

resolve the dominant and extremely rough fjord valley/mountain scale

of the GRIS margin. This is another consideration that is never

discussed.

# specific comments #####

# need surface velocity map comparison between a present-day simulation and that observed, eg as in

Joughin et al (2018, https://tc.copernicus.org/articles/12/2211/2018/)

Ice streams in particular are relevant, as they are characterized by

sliding of ice at the base due to a till that is saturated with water.

# The dominance of subglacial till deformation is not the case (unless

past inferences have significantly changed) for Ilulissat Glacier, or

does your definition of ice stream exclude this significant component of GRIS ice

drainage?

The model is first run to equilibrium for 400 ka to arrive at an

initial state (Fig. 1) corresponding to the present-day GrIS

# This choice is never justified, I can see pros/cons, but at least

it should be made clear this is not meant to represent the current

(non-equilibrium) state of the GRIS, even though the converse is claimed.

there is a geothermal heat flux boundary condition imposed 2 km below

the surface

# What is the specification for this boundary condition and how is

this specific choice justified given the large uncertainties for this

field. Aside from a sensitivity test, listing of simulated present-day

basal temperatures at key ice core sites would provide some sense of

how reasonable the chosen boundary condition is.

Studies of oscillatory behaviour in ice sheets include parameterized

models (Oerlemans, 1983; Fowler and Johnson, 1996; Payne, 1995; Robel

et al., 2013) and comprehensive ice-sheet models with both idealized

geometries (Calov et al., 2010; Van Pelt and Oerlemans, 2012; Feldmann

and Levermann, 2017) and realistic topographies (Papa et al., 2006;

Roberts et al., 230 2016; Schannwell et al., 2023).

# insufficient review of relevant litterature

# The basal topography has insufficient colour contrast in all the map plots

# figure 7 b,e and 8 a,j need geographic or km grid axes.

The conclusions are limited by the amounts and types of simulations

conducted. An obvious next step is to repeat experiments using a

different grid size to observe the dependence, if any, of these

oscillations on the model domain.

# This needs to be addressed in this study. I could see repeating the

full set of experiments would be expensive, but at least a few tests

would provide some hint as to what extent the current results may just

be a numerical artefact.

While the variability of the oscillations seen in this study seems

similar to the build-up/surge variability 330 seen in most simulations

of HEs in the Laurentide ice sheet (LIS) of the LGM (Calov et al.,

2002; Papa et al., 2006; Roberts et al., 2016; Ziemen et al., 2019;

Schannwell et al., 2023), they do not match exactly. Most notably, the

Hudson ice stream in those studies displays a more gradual increase in

ice volume followed by a sudden surge.

# Hank and Tarasov (2024) should also be cited.