The comments are enclosed as a PDF file.

The paper by Harteg et al “Quantifying resilience in non-autonomous and stochastic Earth system dynamics with application to glacial-interglacial cycles” considers modelled output of the earlier introduced glaciation model by Talento and Ganopolski (2021), which the authors modify with stochastic perturbations to investigate its response and apply two measures of resilience.

The paper is interesting, and I suggest that it can be published after a revision. My comments are as follows.

The modelled ice volume anomaly is normalised to zero at preindustrial level and unity at the last glacial maximum. Then the parameters of tolerance, etc. are set to proportional values to serve this specific dataset. I think if the authors like to achieve generality with the novel measures they introduce, they need to scale them with parameters of an arbitrary time series that may be considered in future. This generality would help adoption of the measures in the research community. What the authors demonstrate in Fig.2 is the choice of the parameter value based on the fluctuations of the data. Can this be derived in general terms based on the properties of the input data?

I disagree with the statement (lines 210-213) that what rapidly happens during deglaciations has anything to do with resilience. Resilience means stability of a system state, when the system trajectory reliably fluctuates inside sufficiently deep potential well. Deglaciation is not that kind of behaviour; rather it is a transition from one state to another. If you drop one or two or three stones from a tower, they all will do something very similar, but this will not be resilience. I understand that the authors would like to interpret resilience as “the capacity of the ice-climate-carbon system to remain close to its characteristic large-scale path in face of external shocks” (lines 35-36). But if the planet suddenly becomes hot, massive ice melting is alike a falling stone, and there are no less or more “resilient” alternatives to that.

I think what the resilience function offers in Fig.3 is a clear marking of deglaciation episodes (except N13) under the specific stochastic perturbation, which is interesting. The model was perturbed with a minor stochastic term, thus producing sufficiently similar trajectories driven by the planetary forcing – in particular, reproducing massive deglaciation when the coupled variables create conditions for that. In the context of the definition quoted above (lines 35-36), this stochastic perturbation is indeed resilient for a large number of simulations. Does it mean that the Earth system is resilient during the deglaciation?

In the plot with the resilience function (Fig.3), I can understand the meaning of most of the maxima of the function except the one at 420-400 kyr. There, perturbed and reference datasets (upper panel) are in discord, and it is difficult to see these as staying within a small tolerance. Is there an explanation for that?

In the definition of the return time (Eq.9), in case of very long return times, the value of infimum will be the length of the considered time series. In the text, you mention that you run simulations 1000 kyr longer than 800-length series (line 230), and this explains why in Fig.4, lower panel, at 100kyr the return time in blue shaded area reaches 300 or so in value, while the remaining series is much shorter than that. To avoid confusion, I think the caption should have a note that the simulations were 1800-long rather than 800-long.

FURTHER COMMENTS

Expression “single-event ensemble” is confusing, as it suggests that an ensemble may consist of a single event (which is not what the authors mean). It is better to say “single-event-perturbed ensemble”

What the authors call “Reference Adherence Ratio” is actually not a ratio (a value) but a function of time. Therefore, it is better to name it Reference Adherence Ratio Function (RARF).

In Figure 2, panels b with spectra – in the 4^th from top spectrum, there is an appearing state at about 80kyr, which later disappears. Can the authors explain this effect?

Line 195, it would be helpful to mention that both ensemble and reference are modelled datasets (the reference is not a reconstruction but rather a deterministic model trajectory).

The preprint by Lucarini and Chekroun has already been published: https://www.nature.com/articles/s42254-023-00650-8

MINOR COMMENTS

LaTeX typesetting of text after formulas in display mode: the authors apparently leave a blank line, which leads to indentation of the next line, although it continues the sentence started by the equation. This happens in lines 102, 111, 115, 126, 172, 199. This issue will disappear if the authors remove the blank lines after the display formulas preceding these lines.

Line 24: I am not sure that “revitalise” is the right term for a planetary resilience. “Reinforce” may be a better term.

Line 117, paleo-reconstruction of what? (state the variable).

Line 175, “We apply this” – what this?

In Fig.2, vertical axes are missing

Line 316, I think system state cannot be “novel”; I would call them “previously unsampled”. Or simply “new”

Line 370, “RAR sensitivity to the tolerance value”