Changes in stability and jumps in Dansgaard–Oeschger events: a data analysis aided by the Kramers–Moyal equation

. 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 (cid:58)(cid:58)(cid:58)(cid:58)(cid:58) During (cid:58)(cid:58)(cid:58) the (cid:58)(cid:58)(cid:58)(cid:58) last (cid:58)(cid:58)(cid:58)(cid:58)(cid:58)

Lastly, the analysis in this work is based on Markovian stochastic processes.To assess whether the data is Markovian, we analyse the auto-covariance function of the increments of the detrended data .The covariance is largely zero everywhere, except for a weak anti-correlation at the shortest lag, supporting the assumption that the data is approximately Markovian (see App. ??) .

Methods
In this section, we present the methods used to study the coupled δ 18 O and dust dynamics, drawing on the theory of stochastic processes and stochastic differential equations.Starting from the well-known Langevin process, we introduce the broader Kramers-Moyal (KM) framework.The Fokker-Planck equation associated with a Langevin process will be expanded to the more general Kramers-Moyal equation in one and two dimensions.While the KM setting is suited for the treatment of complex noise, the interpretation of a deterministic drift component is preserved.This allows us to investigate the stability configuration of the coupled δ 18 O-dust system.We begin by shortly introducing one-dimensional stationary Markovian stochastic processes.

The one-dimensional Kramers-Moyal equation
Stochastic processes can be either described in terms of random variables, following a stochastic differential equation as introduced above, or in terms of the evolution of their conditional probability density function p(x, t|x ′ , t ′ ), following a partial differential equation.If a single particle's motion is governed by the Langevin equation, its probability density function p(x, t|x ′ , t ′ ) evolves according to the Fokker-Planck equation, given by Fokker (1913Fokker ( , 1914)); Planck (1917) which we consider in a stationary case as above, i.e.without explicit time dependence of the coefficients D 1 (x) and D 2 (x).
These coefficients directly relate to the drift and diffusion terms given in Eq. (??) by The Fokker-Planck equation can only describe continuous processes (Risken and Frank, 1996;Stemler et al., 2007;Gardiner, 2009;Taba .Giving up the condition of continuity, the temporal evolution of the conditional probability density follows the Kramers-Moyal equation Kramers (1940); Moyal (1949); van Kampen (1961) where D m (x) denotes the mth Kramers-Moyal (KM) coefficient, defined from the corresponding conditional moments M m (x, τ ) of the variable x and a time-lag τ , i.e.

Distinguishing between continuous and discontinuous processes in time series
The striking difference between the Fokker-Planck equation (??) and the Kramers-Moyal equation (??) is the presence of higher-order Kramers-Moyal coefficients D m (x), m > 2, which arise as the direct consequence of discontinuities in a Markovian stochastic process.
Consider the KM coefficients D m (x) estimated from the records of a stochastic process as outlined above.A first metric to discern whether this process is discontinuous is to evaluate the ratio of the fourth KM coefficient to the second one , . This roughly compares the size of the the discontinuous paths with the size of the diffusive effects .Thus, values of the ratio D 4 (x)/D 2 (x) close to zero imply continuous sample paths with no jumps in the data.Values larger than zero, i.e.
The numerical analysis was

Results
This section first discusses the non-parametric estimate of the KM coefficients of the isolated dust and δ 18 O records in a one-dimensional setting.Subsequently, the KM analysis of the two-dimensional coupled system is presented in detail.
Given the high correlation between the dust and the δ 18 O records, the differences in the reconstructed potentials and the ratio between fourth and second KM coefficient are remarkable.At first sight, the monostability of the reconstructed δ 18 O potential contradicts the apparent two regime nature of the time series.There are two possible explanations for this discrepancy: First, regime switching of monostable stochastic process can be achieved through complex noise structures (e.g., Lévy-like noise, generalised Fokker-Planck equations, or fractal motions) (Chechkin et al., 2003(Chechkin et al., , 2004;;Metzler and Klafter, 2004).Secondly, a similar effect can be obtained in a two-dimensional setting if the dynamics of one dynamical variable explicitly depends on the other, which would be impossible to judge from the one-dimensional analysis presented so far.Thus, within the limits of this analysis -that is assuming that the process is Markovian and stationary and that the system under study is fully represented by dust and δ 18 O (no coupling to further hidden variables) -the source of the regime switching must either be endowed by more complex noise processes or by the coupling between the dust and ::::::::: coefficient ::::::: dictates ::: the ::::::::::: deterministic :::::: motion ::: of ::: the :::::: system :::: along ::: the :::: dust :::::::: direction; :::::: therein : the δ 18 O systems.
While we have found the ratio between the fourth and second KM coefficient to be negligible in the case of the dust record, for δ 18 O our analysis yields D 4 (x)/D 2 (x) ∼ 1.We remind the reader that, for a continuous stochastic process x(t), all KM coefficients D m (x) = 0, m > 2, according to Pawula's theorem.When dealing with real world data this is never strictly the case, of course, thus examining the 'ratio of jumps to diffusive motion', that is D 4 (x)/D 2 (x), serves only as a first indication if the process is continuous or not.Still, our results suggest that the dust record can be regarded as a realisation of a continuous stochastic process on the time scale of 5 yr, while the δ 18 O is likely to comprise discontinuities on this time scale.We use the stricter Lehnertz-Tabar Q-ratio to underpin our assessment further.In with respect to τ for the δ 18 O record, suggesting that this stochastic process includes jumps.In contrast, we observe a linear relation between Q(x, τ ) and τ for the dust count, suggesting a purely diffusive process without jumps.We note here that the presence of correlated forms of noise is also sufficient to generate higher-order Kramers-Moyal coefficients (though not affect their scaling or the Q-ratio).However, we exclude this option as the auto-correlation of the increments of the data shows no correlations apart from the shortest increment, as seen in App.??.
The Lehnertz-Tabar Q-ratio of the dust and the δ 18 O concentration, following Eq.(??), in a double-logarithmic scale.For the δ 18 O one observes a constant relation of Q(x, τ ) with τ , indicating that this time series is the realisation of a jumpy (discontinuous) processes.The dust concentration exhibits a linear relation with τ , thus is a purely diffusive process.The state ) is chosen at the maximum of the distribution of the time series.
First we inspect the drift coefficients as before in the one dimensional setting, and reconstruct the conditional potentials according to Eq. (??).This yields two two-dimensional scalar fields whose physical explanatory power is limited to one direction each.
Two-dimensional PDF, potential landscapes, and vector fields.In (a) the PDF.The contour indicates a cutoff > 0.015 of the PDF, the state space we consider henceforth.The dotted elements are the records, separated into stadials (GS) and interstadials (GI).In (b) the effective vector field F eff and in the inset F .In (c) the potential landscape V 1,0 (x 1 |x 2 ) of the dust, conditioned on the δ 18 O.The inset shows D 1,0 (x 1 , x 2 ).In (d) the potential landscape V 0,1 (x 2 |x 1 ) of the δ 18 O, conditioned on the dust.The inset shows D 0,1 (x 1 , x 2 ).For V 0,1 (x 2 |x 1 ), in (d), one finds that the location of the minimum of the landscape changes with the value of the dust conentration and undergoes no bifurcation itself.The system is always mono-stable, yet the minimum is not fixed in state space.In stark contrast, the dust landscape V 1,0 (x 1 |x 2 ), in (c), can show up to three fixed points, depending on the value of δ 18 O.The system exhibits a double-fold bifurcation, transitioning from a single (stable) fixed point for negative values of δ 18 O, bifurcating to three fixed points (two stable), and again returning to a single (stable) fixed point for positive values of δ 18 O This offers a good explanation of the apparent 'regime switching' in the dust record, as the system has two stable fixed points co-existing in some regions of the state space.In (b), the effective vector field shows the direction a conceptual particle follows in this two-dimensional space, telling us how δ 18 O and dust interact and the expected trajectory the coupled system follows.
The study of the conditional potentials provided insight in the underlying static stability configuration of the coupled δ 18 O-dust system.However, to gain a more specific understanding of the dynamics, we investigate the (effective) vector field of motion F (F eff ) in the next step.
The approach was chosen to disclose the dynamical features of this two-dimensional system.In the following, we discuss how our study relates to previously published investigations of the same data and how it contributes to the broad discourse on DO variability.
Although we obtain slightly different and in parts opposing result, our study ties in naturally with previous data-driven analysis of Greenland ice core proxy data.After the work presented by Boers et al. (Boers et al., 2017) it is only the second study to follow a two-dimensional inverse modelling approach with respect to Greenland ice core data.
The non-vanishing fourth KM coefficient in the δ 18 O, which indicates forcing beyond typical Gaussian white noise, could point to an external trigger that directly acts on the Greenland temperatures.A sudden shift in the latter could then entail a regime switch in the atmospheric configuration.However, the interpretation of the fourth KM coefficient is not straightforward and depends on the exact choice of the stochastic process model .The role of discontinuities in the δ 18 O record merits further investigation.Moreover, it should be mentioned that non-Markovian processes, as proposed in Ref. (Boers et al., 2017), can also give rise to higher-order KM coefficients.
The results obtained in our analysis do not give a clear answer to the question for the exact mechanism that triggered DO events.In principle, the revealed double-fold bifurcation would allow for bifurcation-induced transitions and thus for a limit-cycle behaviour.However, the record show that system does not track the stable fixed point branches until the bifurcation points, but tends to transition earlier (not shown) .Also, the structure of the δ 18 O drift is incompatible with a deterministic cyclic motion in the dust-δ 18 O plane.In fact, the specific structure of the double-fold bifurcation leaves room for a weak barrier between stadial and interstadial states in the vicinity of the bifurcation point, thus creating a 'channel'-like passage, through which the system passes.

Conclusion
In this article, we ::: We : have analysed the records of δ 18 O :::: ratios : and dust concentrations from the NGRIP ice core from a data-driven perspective (Ruth et al., 2003;North Greenland Ice Core Projects members, 2004;Gkinis et al., 2014).The central point of our study was to examine the stability configuration of the coupled δ 18 O-dust process by reconstructing its potential landscape.For this aim we utilised the Kramers-Moyal equation which generalises the Fokker-Planck equation in the sense that higher-order Kramers-Moyal coefficients can be related to discontinuities in the stochastic processes.
In a first step, a standard one-dimensional Kramers-Moyal analysis revealed a monostable potential for the isolated δ 18 O record and a bistable one for the dust.This finding calls the prevailing understanding that the Greenland ice core δ 18 O record stem from bistable dynamics into question.The qualitative difference between the reconstructed potentails is remarkable given the high co-variability of the two time series and their synchronous two-regime character, which dominates not only the dust but also the δ 18 O record.Moreover, we found non-vanishing higher-order Kramers-Moyal coefficients for δ 18 O, indicating the presence of discontinuities in the record, assuming the process is truly Markovian.This renders the Langevin equation unsuited to fully describe the underlying process and requires the addition of jumps.In contrast, according to our analysis, the isolated dust record is a continuous process that is described well by the Langevin equation.
, with respect to the parameters a, b GI , and b GS (correspondingly for dust).

Markov property of the data
The Markov property of the data, central in this work, is a necessary property whilst designing Markovian stochastic models to describe any paleo-climatic data, as the name suggests.These include the most commonly used stochastic models, e.g., Langevin processes, based on having independent increments of the data (Risken and Frank, 1996;Friedrich et al., 2011).This is best understood by examining the Chapman-Kolmogorov equation of the joint probability densities p i1,...,in (x 1 , . . ., x n ), with {x i } a collection of random variables ordered with t where the Markov property here allows us to separate each joint probability density as being solely a function of two adjacent segments.One of the most straightforward ways of evaluating the Markov property for data is to examine the auto-correlation function of the increments of the data at the shortest incremental distance.That is, take the data x t , construct the differences ∆x t = x t+1 − x t , and obtain the auto-correlation function ρ(τ ) here µ is the mean and σ 2 the variance of ∆x t .In Fig. ?? we display the auto-correlation functions ρ(τ ) of the two time series, which is only lightly anti-correlated at the shortest lag of τ = 5y.
Here, we include a small note of caution for the interested reader.Pre-processing paleo-climatic data is usually implemented to reduce the noise or remove short or long term trends.A common method to remove long-term trends in the records is to apply a low-pass filter.This will invariantly lead to spurious correlations, thus it should be considered with care, dependent on the data analysis techniques employed.In our case, this would be disastrous.Low-pass filtering the data would create spurious correlations in the incremental time series and Markovianity would be lost.
Autocorrelation ρ(τ ) of the increments ∆x t of δ 18 O and dust records.Both records show a weak anti-correlation at the shortest lag τ = 5y, and no correlation for τ > 5y.We thus consider the data Markovian.
Data availability.All ice core data was obtained from the website of the Niels Bohr Institute of the University of Copenhagen (https://ww w.iceandclimate.nbi.ku.dk/data/); the detailed links are indicated below.The original measurements of δ 18 O ratios and dust concentrations go back to (North Greenland Ice Core Projects members, 2004) and (Ruth et al., 2003), respectively.The 5 cm resolution δ 18 O ratio and dust concentration data together with corresponding GICC05 ages used for this study can be downloaded from https://www.iceandclimate.nbi.In order to carry out the estimation in Eq. ?? we map each data point in the corresponding state space to a kernel density and then take a weighted average over all data points Alike selecting the number of bins in a histogram, when employing kernel-density estimation with an Nadaraya-Watson estimator for the Kramers-Moyal coefficients D m (x), one needs to select both a kernel and a bandwidth (Nadaraya, 1964;Watson, 1964;La .Firstly, the choice of the kernel is the choice of a function K(x) for the estimator f h (x), where h is the bandwidth at a point for a collection {x i } of n random variables.The kernel K(x) is such that K(x) = 1/hK(x/h) and is normalisable ∞ −∞ K(x)dx = 1 (Tabar .The bandwidth is equivalent to the selection of the number of bins, except that binning in a histogram is always "placing numbers into non-overlapping boxes".The optimal kernel is the commonly denoted Epanechnikov kernel (Epanechnikov, 1967) , but Gaussian kernels can be used as well.These nevertheless require a compact support in (−∞, ∞), thus on a computer they require some sort of truncation (even if Fourier space, as the Gaussian shape remains unchanged).
The selection of an appropriate bandwidth h can be aided -unlike the selection of the number of bins -by the Silverman's rule-of-thumb (Silverman, 1998), given by h S = 4σ 5 3n 1 5 , where again σ 2 is the variance of the time series.In In order to correctly retrieve from data the Kramers-Moyal coefficients, we need to evaluate the operation in the Fokker-Planck equation Eq. (??).In fact, we showed that this equation is not sufficient to describe the fast transitions in the δ where M m (x, τ ) is the m-order conditional moment, i.e.
which we introduced in Eq. (??) in a similar notation.If we limit m ≤ 2, we are truly talking about a Langevin process described by the Fokker-Planck equation, with L FP the formal Fokker-Planck operator.We also saw that we can generalise the problem and not truncate the terms at second order, thus including an infinite series of conditional moment M m (x, τ ) would give rise to the Kramers-Moyal equation and the Kramers-Moyal operator L KM .The subsequent second-order correction is showcased here for the Fokker-Planck equation, based on Ref. (Gottschall and Peinke, 2008;Rydin Gorjão et al., 2021).For sake of coherence, we utilise here the second-order corrections to show that the second Kramers-Moyal coefficient -the diffusion strength -can be seen as constant, i.e. not state depended.
In order to solve Eq. (??), one takes the formal step considering an initial conditions δ(x − x ′ ) as a starting point and employing the exponential representation of the operator, which we can decompose it into a power series as From here we consider the first-order and second-order approximation, i.e. truncation of the operator, as Considering only the first-order, ∼ τ , we recover the well-known relation between the conditional moments and the Kramers-Moyal coefficients, given by M m (x, τ ) (m!)τ .
If we now include the second-order approximation, i.e. we consider terms up to ∼ τ 2 , we obtain a corrective term for the second Kramers-Moyal coefficient We employ this correction to our examination to show that the diffusion coefficient, i.e. the amplitude of the fluctuations, is constant in space.In Fig. ??(c) and (g) we display both the first-order and the corrected, second-order diffusion coefficient.In this we can see that utilising solely the first-order correction could lead us to erroneously consider the diffusion term as state dependent (having a parabolic shape), suggesting a multiplicative noise.By implementing the second-order corrective terms a considerable improvement of the estimation is achieved, to what we judge to be simple additive (not state dependent) noise in the records.lag :::::: τ = 5y, ::: and ::: no :::::::: correlation :: for :::::: τ > 5y.::: We :::: thus :::::: consider ::: the ::: data ::::::::: Markovian.
Fig.1(a) for comparison).We find the second KM coefficient to be fairly constant (Fig.3 (c)) and the ratio between fourth Fig. 3 (e)  and (f) display estimates of the drift and the potential landscape for δ 18 O, respectively.Most prominently, the drift has only a single stable fixed point (zero-crossing of the drift), or equivalently the potential function exhibits only a single well.The second KM coefficient is mostly constant, like in the case of dust (Fig.3 (g)).With respect to the normalised units, the first and second KM coefficients of δ 18 O exceed their counterparts for dust by factors of approximately 3-4 and 10, respectively.This indicates that δ 18 O exhibits faster dynamics than dust.Moreover, we find that the fourth KM coefficient D 4 (x) for the δ 18 O is of the same magnitude as the second KM coefficient D 2 (x) ( 4.1 :::::::::: Double-fold :::::::::: bifurcation :: of ::: the :::: dust ::: Fig. ??, we clearly see a constant relation of Q(x, τ ) Fig. :::: 3 (c)).With the position of these stable fixed points depending continuously on δ 18 O ::::: ratios, we find here the characteristic form of a double-fold bifurcation.Fig. ?? (d) suggests

Fig. ? ?
Fig. ?? illustrates the detrending scheme for both time series.Due to the two regime nature of the time series, a simple linear regression overestimates the temperature dependencies (see Fig. ??(c) and (d), dashed line).Instead, we separate the data from Greenland Stadials and Greenland Interstadials and then minimise the expression t i ∈ GS (GI) indicates that a given time t i falls into a Stadial (Interstadial) period.The index i runs over all data points and N denotes the total number of points.The resulting a is used to detrend the original data with respect to the temperature.The detrended data is then normalised by subtraction of its mean and division by the difference between the mean Stadial and mean Interstadial values.Removal of a linear trend in the NGRIP δ 18 O and dust time series(North Greenland Ice Core Projects members, 2004) with respect to a global average surface temperature reconstruction(Snyder, 2016).In panel (a) both original δ 18 O (blue) as well as detrended and normalised (purple) are shown.Idem for the dust record in panel (b) (dark green and light green, respectively).The background temperature given in anomalies to present day climate is shown in both aforementioned panels (red).Panels (c) and (d) show a scatter plot the original δ 18 O and dust data with respect to temporarily corresponding temperature anomalies, respectively.Data from Interstadials (Stadials) is shown in orange (light blue).The black dashed line results from a simple linear fit to the entire data, while the continuous black lines correspond to the fitting scheme that uses a single slope but two different offsets to separately fit the Stadial and Interstadial data.
ku.dk/data/NGRIP_d18O_and_dust_5cm.xls(last access: 11.August 2022).The δ 18 O data shown in Fig.1with 20 yr resolution that cover the period 122-10 kyr b2k are available from https://www.iceandclimate.nbi.ku.dk/data/GICC05modelext_GRIP_and_GISP2_and_res ampled_data_series_Seierstad_et_al._2014_version_10Dec2014-2.xlsx(last access: 08.08.2022) and were published in conjunction with the work byRasmussen et al. (2014) andSeierstad et al. (2014).The corresponding dust data, also shown in Fig.1and covering the period 108-10 kyr b2k, can be retrieved from https://www.iceandclimate.nbi.ku.dk/data/NGRIP_dust_on_GICC05_20y_december2014.txt.The global average surface temperature reconstructions provided bySnyder (2016) and used here for the detrending were retrieved from https: //static-content.springer.com/esm/art%3A10.1038%2Fnature19798/MediaObjects/41586_2016_BFnature19798_MOESM258_ESM.xlsx(last access: 08.08.2022) 7 Nadaraya-Watson estimator of the Kramers-Moyal coefficients and bandwidth selection Fig. ?? three different bandwidths are used to evaluated the various KM coefficient, as given in Fig. ?? The bandwidths are the optimal bandwidth given by the Silverman's rule-of-thumb h S , three times h S , and one-third h S .The effect of the bandwidth selection h S on the KM estimations, in identical fashion to Fig.3.The non-parametric estimates of the first KM coefficient D 1 (x), the associated potential landscape V (x), the second KM coefficient D 2 (x), and the ratio of the fourth to the second KM coefficient D 4 (x)/D 2 (x).Left column for dust, right column for δ 18 O.Three bandwidths used for the Nadaraya-Watson kernel-density estimator: the optimal Silverman's rule-of-thumb h S , three times h S , and one-third h S .The Nadaraya-Watson kernel-density estimator's bandwidths h S for δ 18 O is 0.131 and for dust 0.103.In all cases, the interpretation of the estimator remains the same: bi-stability in the dust, mono-stability in the δ 18 O.In (c) and (g) are included the first-order estimator for the second KM coefficient D 2 (x), i.e. without corrective terms, discussed in App.??.Note that in neither of the examples with different bandwidths we notice a change of the potential shape of the various records.The mono-stability of the potential V (x) of δ 18 O is persistent, as is the bi-stability of the potential V (x) of the dust concentration.