Changes in stability and jumps in Dansgaard–Oeschger events: a data analysis aided by the Kramers–Moyal equation
 ^{1}Forschungszentrum Jülich, Institute for Energy and Climate Research  Systems Analysis and Technology Evaluation (IEKSTE), 52428 Jülich, Germany
 ^{2}Institute for Theoretical Physics, University of Cologne, 50937 Köln, Germany
 ^{3}German Aerospace Center (DLR), Institute of Networked Energy Systems, Oldenburg, Germany
 ^{4}Department of Computer Science, OsloMet – Oslo Metropolitan University, N0130 Oslo, Norway
 ^{5}Research Domain IV – Complexity Science, Potsdam Institute for Climate Impact Research, Telegrafenberg A31, 14473 Potsdam, Germany
 ^{6}Department of Mathematics and Computer Science, Freie Universität Berlin, Berlin, Germany
 ^{7}Department of Physics, HumboldtUniversität zu Berlin, Newtonstraße 15, 12489 Berlin, Germany
 ^{8}NordSTAR – Nordic Center for Sustainable and Trustworthy AI Research, N0166 Oslo, Norway
 ^{9}Artificial Intelligence Lab, Oslo Metropolitan University, N0166 Oslo, Norway
 ^{10}Technical University of Munich, Germany; School of Engenieering & Design, Earth System Modelling
 ^{11}Global Systems Institute and Department of Mathematics, University of Exeter, United Kingdom
 ^{}These authors contributed equally to this work.
 ^{1}Forschungszentrum Jülich, Institute for Energy and Climate Research  Systems Analysis and Technology Evaluation (IEKSTE), 52428 Jülich, Germany
 ^{2}Institute for Theoretical Physics, University of Cologne, 50937 Köln, Germany
 ^{3}German Aerospace Center (DLR), Institute of Networked Energy Systems, Oldenburg, Germany
 ^{4}Department of Computer Science, OsloMet – Oslo Metropolitan University, N0130 Oslo, Norway
 ^{5}Research Domain IV – Complexity Science, Potsdam Institute for Climate Impact Research, Telegrafenberg A31, 14473 Potsdam, Germany
 ^{6}Department of Mathematics and Computer Science, Freie Universität Berlin, Berlin, Germany
 ^{7}Department of Physics, HumboldtUniversität zu Berlin, Newtonstraße 15, 12489 Berlin, Germany
 ^{8}NordSTAR – Nordic Center for Sustainable and Trustworthy AI Research, N0166 Oslo, Norway
 ^{9}Artificial Intelligence Lab, Oslo Metropolitan University, N0166 Oslo, Norway
 ^{10}Technical University of Munich, Germany; School of Engenieering & Design, Earth System Modelling
 ^{11}Global Systems Institute and Department of Mathematics, University of Exeter, United Kingdom
 ^{}These authors contributed equally to this work.
Abstract. Dansgaard–Oeschger (DO) events are sudden climatic shifts from cold to substantially milder conditions in the arctic region that occurred during previous glacial intervals. They can be most clearly identified in paleoclimate records of δ^{18}O and dust concentrations from Greenland ice cores, which serve as proxies for temperature and atmospheric circulation patterns, respectively. The existence of stadial (cold) and interstadial (milder) phases is typically attributed to a bistability of the North Atlantic climate system allowing for rapid transitions from the first to the latter and a more gentle yet still fairly abrupt reverse shift from the latter to the first. However, the underlying physical mechanisms causing these transitions remain debated. Here, we conduct a datadriven analysis of the Greenland temperature and atmospheric circulation proxies under the purview of stochastic processes. Based on the Kramers–Moyal equation we present a onedimensional and twodimensional derivation of the proxies' drift and diffusion terms, which unravels the features of the climate system's stability landscape. Our results show that: (1) in contrast to common assumptions, the δ^{18}O proxy results from a monostable process, and transitions occur in the record only due to the coupling to other variables; (2) conditioned on δ^{18}O the dust concentrations exhibit both mono and bistable states, transitioning between them via a doublefold bifurcation; (3) the δ^{18}O record is discontinuous in nature, and mathematically requires an interpretation beyond the classical Langevin equation. These findings can help understand candidate mechanisms underlying these archetypal examples of abrupt climate changes.
Leonardo Rydin Gorjão et al.
Status: closed (peer review stopped)

RC1: 'Comment on esd202195', Peter Ditlevsen, 04 Jan 2022
In this paper two paleoclimatic ice core records are analyzed. These, the water isotope, d18O, and the dust concentration records are analyzed for the period 5927 kyr BP, which is the glacial period dominated by regular occurrences of DansgaardOeschger events. There are two major points in the paper: Firstly, the data are modelled as a stochastic process using the KramersMoyal equation to investigate the importance of (discontinuous) jumps in the noise. Secondly, the two records are modelled as a twodimensional joined process. It seems to me that the two points are only loosely related, and the authors could consider presenting them in two separate papers.
I enjoyed reading the paper and find it publishable. However, there are a few issues below calling for revisions before publication.
I have two major concerns regarding the two parts, and some minor points:
As to the first point, I have not seen the KramersMoyal (KM) equation applied to these data before, so this is a novel approach. The equation, for which the Taylor expansion of the conditional probability density function is taken to higher order than two, covers the case where the noise term in the governing Langevin equation is not gaussian, but contains jumps. It is stated that in the case of Levy processes the FokkerPlanck equation does not apply. Actually, for the most relevant class of Levy processes, the alphastable Levy processes an extension of the FokkerPlanck equation based on the characteristic function of the alphastable process exists (see: Samorodnitsky and Taqqu (1994) ‘Stable nongaussian processes, Chapman and Hall, NY. or Ditlevsen, PRE, 60, 172179). The challenge in applying the KM equation is the estimate of the higher order coefficients (eq. 6) for the data series: Since the higher order terms are (increasingly) dominated by the extremes in the increments the finite time series very quickly “dilutes”. A main result (eq. 9 and figure 3) includes the sixth moment of the increments. I thus miss an analysis of uncertainties and reliability in these estimates. I find that this is essential for publishing this (nice!) result.
Another point which could be given a little more attention is the fact (as also correctly stated) that the strong timeasymmetry in the data (the sawtooth shape) cannot be captured by the model. How does this influence the relevance in including higher order terms (higher than second order) in eq. 5?
As to the second point, the major results are presented in figure 4. Obviously, when considering a onedimensional record, the drift can always be seen as a result of a potential. This is not the case in two  and higher dimensions, where gradient drift is a nongeneric case. I’m sure that the authors are aware of this, the drift is a twodimensional flow field, as also shown in the small inserts in the subplots of figure 4. I find the construct of pseudoonedimensional potentials (V(x_1x_2)) both confusing and useless. I suggest that the authors consider abandoning this all together (as well as the notion of a potential landscape). The interpretation in figure 4(c) of a double fold bifurcation is obscure, and I believe wrong. What I find interesting is the scatter plot in 4(a), which nicely explain the results in figure 2, namely that the stationary distribution for the dust is bimodal while it is unimodal for the d18O: This corresponds to the marginal distributions in 4(a) (projections onto the axis. The authors could consider analyzing a rotation (linear combination of the two variables) of the data along an axis connecting the two maxima (GS and GI) and a perpendicular direction. In this way one would obtain a “clean” two state dynamics and a “clean” one state perpendicular dynamics. First thing would be to check for independence. Just a suggestion.
Minor points:
Introduction: These data have been analyzed over many years and a lot is known. For a better overview and setting the present work in context, a representative presentation of work done over the years would be useful. There is a strong bias towards very resent publications.
Again: I’m not fan of the “potential landscape” metaphor.
Figure 1: I suggest to plot the dust record upsidedown. This will visualize the strong dependence between the two records, and also make the saw tooth shape in the dust record much more apparent. Make the figure full text width, Ylabel: ln(dust) (no units), d18O (permil). Or normalized w.r.t. std. dev.
L90: A discussion of concentration vs flux of dust could be added
L91: Data is > Data are
L96: This was pointed out by others previously: Rial and Saha (2011), Abrupt Climate Change: Mechanisms, Patterns, and Impacts, Geophysical Monograph Series 193, Mitsui and Crucifix, arXiv:1510.06290, Lohmann and Ditlevsen, 2018, Clim Past 14. ( and probably others).
L110: I do not understand why
L125: Stationarity require certain properties of a(x) (such that the process does not drift to infinity. It is better to denote it “homogenous” or “autonomous”. (same comment on L138).
L126: Langevin equation is continuous > Langevin equation generate realizations which are continuous
Top of Form
L155: There is no Eq 4. (4a and 4b are hardly equations)
L185: The purely continuous process (gaussian process) diffuse proportionally to t^(1/2) not t. That is: (sigma(t)=sqrt(4Dt)). The most natural jump processes in this context are the alphastable processes, they do exhibit similar scaling relations with time (sigma(t)~t^(1/alpha)).
L205: Also referring back to Figure 1: Wouldn’t it be natural to rescale the data before doing the analysis.
L239: I do not understand the statement (which I believe is not correct): This indicates that d18O exhibits faster dynamics than dust. Please explain.
L245: The exotic explanation for monostability through a complex noise structure seems a little out of context here: The reason why the d18O record has a single maximum in the PDF (with a shoulder) is the sawtooth shape of the DO events masking the obvious two state nature of the record.
Section 4.2.1 is obscure to me.
L273 and L288: 4 (d) <> 4(c) .
L456: What happened to the index _t in mu and sigma^2?
I hope my comments are useful.
 AC1: 'Reply on RC1', Leonardo Rydin Gorjão, 25 Feb 2022

RC2: 'Comment on esd202195', Tamas Bodai, 21 Jan 2022
The paper “Changes in stability and jumps in Dansgaard–Oeschger events: a data analysis aided by the Kramers–Moyal equation” analyses d18O and dust data from a Greenland ice core in order to gain further understanding of the famous DansgaardOeschger (DO) events. They preprocess the time series with the aim of establishing a stationary stochastic process. They estimate KramersMoyal (KM) coefficients, which could possibly reveal jumps in the process, outside the framework of the FokkerPlanck equation, corresponding to what can be naively seen as regime transitions. They explore the added value of joint fitting of the d18O and dust data over treating them separately.
I’m not very convinced that the applied methodology is suitable. As far as i see, the authors do not test their nullhypothesis (H0) of a stationary process. I would not think that a stationary process described by the KM equations is consistent with a hypothetical nonstationary process that could not be rejected. Say, we have a nonstationary OrntsetinUhlenbeck process of
dx/dt = a*x + B(t) + c*xi(t), (OU)
where xi is white noise, and B(t) = b*sin(sin(2*t)+t) is some regular nonstationarity. It mimics some regime behaviour with sudden and regular transitions. We can easily see that the pdf of x is bimodal. If we didn't know the underlying process generating eq., and perhaps we somehow overlooked the regularity of the transitions, we might think the underlying model is:
dx/dt = F*x + c*xi(t), (H0)
where F = V(x), V(x) being a doublewell potential function. If we are in the small noise limit, we know that the pdf takes the shape of V(x), so, we could estimate V that way. Furthermore, we can estimate the noise strength ‘c’ in some standard way too. Now the question is if it matters at all that we have an H0 other than the true process OU. It would not matter if in any appreciable way the processes perform the same, i.e., when, loosely speaking, they are consistent or approximately equivalent. For example, we can derive the probability distribution of residence times from H0, and perform a statistical test if our residence time data is consistent with that, or we can reject H0.
We want to perform a socalled crucial experiment (experimentum crucis).
Considering H0 of the authors, another likely feature based on which H0 can be rejected is the saw tooth asymmetry, in particular, that the cold to warm, c2w, transitions are much more rapid than the warm to cold, w2c, ones. Looking at the dust time series of Fig. 1 with much naivity wrt. physics, but with some experience about dynamical systems, i would think that c2w is an attractor crisis, whereas w2c is a noise induced tipping, and there is some slowly drifting control parameter, i.e., a nonstationarity when we exclude that parameter from our state variables.
The authors also make a reference to attractor crisis, in terms of a saddlenode bifurcation, but in some other context. It is the context of slices of a 2D potential function. This does not sound correct.
The paper is very well written in a way, but it doesn’t make for a very pleasant reading journeying through flawed results, starting with the single variable approach, and then — at least as i suspect — even the 2 variable approach.
I attach the pdf of the manuscript with comments saved as annotations. Hopefully the authors find it useful in some way.
Note: I always review nonanonymously, and never make recommendation for or against publication. The recommendation that i make is only to circumvent the rigidity of the submission system, and therefore please consider it void.
Tamas Bodai
 AC2: 'Reply on RC2', Leonardo Rydin Gorjão, 25 Feb 2022
Status: closed (peer review stopped)

RC1: 'Comment on esd202195', Peter Ditlevsen, 04 Jan 2022
In this paper two paleoclimatic ice core records are analyzed. These, the water isotope, d18O, and the dust concentration records are analyzed for the period 5927 kyr BP, which is the glacial period dominated by regular occurrences of DansgaardOeschger events. There are two major points in the paper: Firstly, the data are modelled as a stochastic process using the KramersMoyal equation to investigate the importance of (discontinuous) jumps in the noise. Secondly, the two records are modelled as a twodimensional joined process. It seems to me that the two points are only loosely related, and the authors could consider presenting them in two separate papers.
I enjoyed reading the paper and find it publishable. However, there are a few issues below calling for revisions before publication.
I have two major concerns regarding the two parts, and some minor points:
As to the first point, I have not seen the KramersMoyal (KM) equation applied to these data before, so this is a novel approach. The equation, for which the Taylor expansion of the conditional probability density function is taken to higher order than two, covers the case where the noise term in the governing Langevin equation is not gaussian, but contains jumps. It is stated that in the case of Levy processes the FokkerPlanck equation does not apply. Actually, for the most relevant class of Levy processes, the alphastable Levy processes an extension of the FokkerPlanck equation based on the characteristic function of the alphastable process exists (see: Samorodnitsky and Taqqu (1994) ‘Stable nongaussian processes, Chapman and Hall, NY. or Ditlevsen, PRE, 60, 172179). The challenge in applying the KM equation is the estimate of the higher order coefficients (eq. 6) for the data series: Since the higher order terms are (increasingly) dominated by the extremes in the increments the finite time series very quickly “dilutes”. A main result (eq. 9 and figure 3) includes the sixth moment of the increments. I thus miss an analysis of uncertainties and reliability in these estimates. I find that this is essential for publishing this (nice!) result.
Another point which could be given a little more attention is the fact (as also correctly stated) that the strong timeasymmetry in the data (the sawtooth shape) cannot be captured by the model. How does this influence the relevance in including higher order terms (higher than second order) in eq. 5?
As to the second point, the major results are presented in figure 4. Obviously, when considering a onedimensional record, the drift can always be seen as a result of a potential. This is not the case in two  and higher dimensions, where gradient drift is a nongeneric case. I’m sure that the authors are aware of this, the drift is a twodimensional flow field, as also shown in the small inserts in the subplots of figure 4. I find the construct of pseudoonedimensional potentials (V(x_1x_2)) both confusing and useless. I suggest that the authors consider abandoning this all together (as well as the notion of a potential landscape). The interpretation in figure 4(c) of a double fold bifurcation is obscure, and I believe wrong. What I find interesting is the scatter plot in 4(a), which nicely explain the results in figure 2, namely that the stationary distribution for the dust is bimodal while it is unimodal for the d18O: This corresponds to the marginal distributions in 4(a) (projections onto the axis. The authors could consider analyzing a rotation (linear combination of the two variables) of the data along an axis connecting the two maxima (GS and GI) and a perpendicular direction. In this way one would obtain a “clean” two state dynamics and a “clean” one state perpendicular dynamics. First thing would be to check for independence. Just a suggestion.
Minor points:
Introduction: These data have been analyzed over many years and a lot is known. For a better overview and setting the present work in context, a representative presentation of work done over the years would be useful. There is a strong bias towards very resent publications.
Again: I’m not fan of the “potential landscape” metaphor.
Figure 1: I suggest to plot the dust record upsidedown. This will visualize the strong dependence between the two records, and also make the saw tooth shape in the dust record much more apparent. Make the figure full text width, Ylabel: ln(dust) (no units), d18O (permil). Or normalized w.r.t. std. dev.
L90: A discussion of concentration vs flux of dust could be added
L91: Data is > Data are
L96: This was pointed out by others previously: Rial and Saha (2011), Abrupt Climate Change: Mechanisms, Patterns, and Impacts, Geophysical Monograph Series 193, Mitsui and Crucifix, arXiv:1510.06290, Lohmann and Ditlevsen, 2018, Clim Past 14. ( and probably others).
L110: I do not understand why
L125: Stationarity require certain properties of a(x) (such that the process does not drift to infinity. It is better to denote it “homogenous” or “autonomous”. (same comment on L138).
L126: Langevin equation is continuous > Langevin equation generate realizations which are continuous
Top of Form
L155: There is no Eq 4. (4a and 4b are hardly equations)
L185: The purely continuous process (gaussian process) diffuse proportionally to t^(1/2) not t. That is: (sigma(t)=sqrt(4Dt)). The most natural jump processes in this context are the alphastable processes, they do exhibit similar scaling relations with time (sigma(t)~t^(1/alpha)).
L205: Also referring back to Figure 1: Wouldn’t it be natural to rescale the data before doing the analysis.
L239: I do not understand the statement (which I believe is not correct): This indicates that d18O exhibits faster dynamics than dust. Please explain.
L245: The exotic explanation for monostability through a complex noise structure seems a little out of context here: The reason why the d18O record has a single maximum in the PDF (with a shoulder) is the sawtooth shape of the DO events masking the obvious two state nature of the record.
Section 4.2.1 is obscure to me.
L273 and L288: 4 (d) <> 4(c) .
L456: What happened to the index _t in mu and sigma^2?
I hope my comments are useful.
 AC1: 'Reply on RC1', Leonardo Rydin Gorjão, 25 Feb 2022

RC2: 'Comment on esd202195', Tamas Bodai, 21 Jan 2022
The paper “Changes in stability and jumps in Dansgaard–Oeschger events: a data analysis aided by the Kramers–Moyal equation” analyses d18O and dust data from a Greenland ice core in order to gain further understanding of the famous DansgaardOeschger (DO) events. They preprocess the time series with the aim of establishing a stationary stochastic process. They estimate KramersMoyal (KM) coefficients, which could possibly reveal jumps in the process, outside the framework of the FokkerPlanck equation, corresponding to what can be naively seen as regime transitions. They explore the added value of joint fitting of the d18O and dust data over treating them separately.
I’m not very convinced that the applied methodology is suitable. As far as i see, the authors do not test their nullhypothesis (H0) of a stationary process. I would not think that a stationary process described by the KM equations is consistent with a hypothetical nonstationary process that could not be rejected. Say, we have a nonstationary OrntsetinUhlenbeck process of
dx/dt = a*x + B(t) + c*xi(t), (OU)
where xi is white noise, and B(t) = b*sin(sin(2*t)+t) is some regular nonstationarity. It mimics some regime behaviour with sudden and regular transitions. We can easily see that the pdf of x is bimodal. If we didn't know the underlying process generating eq., and perhaps we somehow overlooked the regularity of the transitions, we might think the underlying model is:
dx/dt = F*x + c*xi(t), (H0)
where F = V(x), V(x) being a doublewell potential function. If we are in the small noise limit, we know that the pdf takes the shape of V(x), so, we could estimate V that way. Furthermore, we can estimate the noise strength ‘c’ in some standard way too. Now the question is if it matters at all that we have an H0 other than the true process OU. It would not matter if in any appreciable way the processes perform the same, i.e., when, loosely speaking, they are consistent or approximately equivalent. For example, we can derive the probability distribution of residence times from H0, and perform a statistical test if our residence time data is consistent with that, or we can reject H0.
We want to perform a socalled crucial experiment (experimentum crucis).
Considering H0 of the authors, another likely feature based on which H0 can be rejected is the saw tooth asymmetry, in particular, that the cold to warm, c2w, transitions are much more rapid than the warm to cold, w2c, ones. Looking at the dust time series of Fig. 1 with much naivity wrt. physics, but with some experience about dynamical systems, i would think that c2w is an attractor crisis, whereas w2c is a noise induced tipping, and there is some slowly drifting control parameter, i.e., a nonstationarity when we exclude that parameter from our state variables.
The authors also make a reference to attractor crisis, in terms of a saddlenode bifurcation, but in some other context. It is the context of slices of a 2D potential function. This does not sound correct.
The paper is very well written in a way, but it doesn’t make for a very pleasant reading journeying through flawed results, starting with the single variable approach, and then — at least as i suspect — even the 2 variable approach.
I attach the pdf of the manuscript with comments saved as annotations. Hopefully the authors find it useful in some way.
Note: I always review nonanonymously, and never make recommendation for or against publication. The recommendation that i make is only to circumvent the rigidity of the submission system, and therefore please consider it void.
Tamas Bodai
 AC2: 'Reply on RC2', Leonardo Rydin Gorjão, 25 Feb 2022
Leonardo Rydin Gorjão et al.
Leonardo Rydin Gorjão et al.
Viewed
HTML  XML  Total  BibTeX  EndNote  

697  153  19  869  11  10 
 HTML: 697
 PDF: 153
 XML: 19
 Total: 869
 BibTeX: 11
 EndNote: 10
Viewed (geographical distribution)
Country  #  Views  % 

Total:  0 
HTML:  0 
PDF:  0 
XML:  0 
 1