the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
| 28 Nov 2025
Quantifying Resilience in Non-Autonomous and Stochastic Earth System Dynamics with Application to Glacial-Interglacial Cycles
Abstract. Understanding Earth resilience—the capacity of the Earth system to absorb and regenerate from perturbations—is key to assessing risks from anthropogenic pressures and sustaining a safe operating space for humanity within planetary boundaries. Classical resilience indicators are designed for autonomous systems with fixed attractors, but the Earth system is fundamentally non-autonomous and out of equilibrium, calling for new ways of defining and quantifying resilience.
Here, we introduce a path-resilience approach that assesses how perturbations deviate from and return to a reference trajectory of a conceptual climate model replicating the glacial-interglacial cycles of the Late Pleistocene. We generate two types of perturbation ensembles: a stochastic ensemble and a single-event ensemble, and compute two complementary metrics: the Reference Adherence Ratio (RAR)—defined as the fraction of stochastic trajectories that remain within a narrow band around the unperturbed trajectory—and return time—defined as the time a single perturbed trajectory takes to return to the reference path. Together, these metrics reveal strong temporal variation in resilience across the glacial-interglacial cycles. We find that RAR increases markedly during deglaciations and peaks in interglacial periods, while return times generally shorten as the system approaches interglacial conditions—indicating that certain phases of the cycles act as convergence zones and potential anchors of Earth system stability. As the Earth system departs from such stable interglacial regimes under ongoing anthropogenic forcing, understanding the resilience of these trajectories—and what it may take to return to them—becomes increasingly important.
These results highlight that resilience in non-autonomous systems is inherently path-dependent and illustrate a promising first step toward its quantification. Further research is needed to develop more general resilience indicators suitable for complex, forced dynamical systems.
Competing interests: Jonathan Donges is a member of the editorial board of Earth System Dynamics.
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 paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors. Views expressed in the text are those of the authors and do not necessarily reflect the views of the publisher.-
Notice on discussion status
The requested preprint has a corresponding peer-reviewed final revised paper. You are encouraged to refer to the final revised version.
-
Preprint
(2266 KB)
-
The requested preprint has a corresponding peer-reviewed final revised paper. You are encouraged to refer to the final revised version.
- Preprint
(2266 KB) - Metadata XML
- BibTeX
- EndNote
- Final revised paper
Journal article(s) based on this preprint
Interactive discussion
Status: closed
-
RC1: 'Comment on egusphere-2025-5058', Takahito Mitsui, 16 Dec 2025
- AC1: 'Reply on RC1', Jakob Harteg, 18 Mar 2026
-
RC2: 'Comment on egusphere-2025-5058', Anonymous Referee #2, 23 Jan 2026
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 4th 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”
Citation: https://doi.org/10.5194/egusphere-2025-5058-RC2 - AC2: 'Reply on RC2', Jakob Harteg, 18 Mar 2026
Interactive discussion
Status: closed
-
RC1: 'Comment on egusphere-2025-5058', Takahito Mitsui, 16 Dec 2025
- AC1: 'Reply on RC1', Jakob Harteg, 18 Mar 2026
-
RC2: 'Comment on egusphere-2025-5058', Anonymous Referee #2, 23 Jan 2026
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 4th 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”
Citation: https://doi.org/10.5194/egusphere-2025-5058-RC2 - AC2: 'Reply on RC2', Jakob Harteg, 18 Mar 2026
Peer review completion
Journal article(s) based on this preprint
Model code and software
Simulation and model analysis code Jakob Harteg https://doi.org/10.5281/zenodo.16603222
Viewed
| HTML | XML | Total | BibTeX | EndNote | |
|---|---|---|---|---|---|
| 1,778 | 574 | 130 | 2,482 | 116 | 122 |
- HTML: 1,778
- PDF: 574
- XML: 130
- Total: 2,482
- BibTeX: 116
- EndNote: 122
Viewed (geographical distribution)
| Country | # | Views | % |
|---|
| Total: | 0 |
| HTML: | 0 |
| PDF: | 0 |
| XML: | 0 |
- 1
Nico Wunderling
Jonathan Donges
The requested preprint has a corresponding peer-reviewed final revised paper. You are encouraged to refer to the final revised version.
- Preprint
(2266 KB) - Metadata XML
The comments are enclosed as a PDF file.