Discontinuous stochastic forcing in Greenland ice core data

ORIGINAL ARTICLE Climate Dynamics (2025) 63:465 https://doi.org/10.1007/s00382-025-07880-9 Keno Riechers, Andreas Morr and Leonardo Rydin Gorjão have contributed equally to this work. Andreas Morr [email protected] 1 Climate Physics, Max Planck Institute for Meteorology, Hamburg, Germany 2 Earth System Modelling, School of Engineering and Design, Technical University of Munich, Munich, Germany 3 Research Domain IV – Complexity Science, Potsdam Institute for Climate Impact Research, 14473 Potsdam, Germany 4 Department of Epileptology, University of Bonn Medical Center, 53105 Bonn, Germany 5 Helmholtz Institute for Radiation and Nuclear Physics, University of Bonn, 53115 Bonn, Germany 6 Interdisciplinary Centre for Complex Systems, University of Bonn, 53175 Bonn, Germany 7 Department of Computer Science, OsloMet – Oslo Metropolitan University, 0130 Oslo, Norway 8 School of Economics, Innovation and Technology, Kristiania University of Applied Sciences, Kirkegata 24-26, 0153 Oslo, Norway 9 Global Systems Institute and Department of Mathematics, University of Exeter, Exeter, UK 10 Institute of Energy and Climate Research (IEK-10), Forschungszentrum Jülich, 52428 Jülich, Germany 11 Institute for Theoretical Physics, University of Cologne, 50937, Köln, Germany 12 Department of Environmental Sciences Faculty of Science, Open University, Heerlen, The Netherlands 13 Faculty of Science and Technology, Norwegian University of Life Sciences, 1432 Ås, Norway Abstract Paleoclimate proxy records from Greenland ice cores, archiving e.g. δ18 O as a proxy for surface temperature, show that sudden climatic shifts called Dansgaard–Oeschger events (DO) occurred repeatedly during the last glacial interval. They comprised substantial warming of the Arctic region from cold to milder conditions. Concomitant abrupt changes in the dust concentrations of the same ice cores suggest that sudden reorganisations of the hemispheric-scale atmospheric circulation have accompanied the warming events. Genuine bistability of the North Atlantic climate system is commonly hypothesised to explain the existence of stadial (cold) and interstadial (milder) periods in Greenland. However, the physical mechanisms that drove abrupt transitions from the stadial to the interstadial state, and more gradual yet still abrupt reverse transitions, remain debated. Here, we conduct a one-dimensional data-driven analysis of the Greenland temperature and atmospheric circulation proxies under the purview of stochastic processes. We take the Kramers–Moyal equation to estimate each proxy’s drift and diffusion terms within a Markovian model framework. We then assess noise contributions beyond Gaussian white noise. The resulting stochastic differential equation (SDE) models feature a monostable drift for the Greenland temperature proxy and a bistable one for the atmospheric circulation proxy. Indicators of discontinuity in stochastic processes suggest to include higher-order terms of the Kramers–Moyal equation when modelling the Greenland temperature proxy’s evolution. This constitutes a qualitative difference in the characteristics of the two time series, which should be further investigated from the standpoint of climate dynamics. Received: 22 February 2025 / Accepted: 15 September 2025 © The Author(s) 2025 Discontinuous stochastic forcing in Greenland ice core data KenoRiechers1· AndreasMorr2,3 · KlausLehnertz4,5,6 · Pedro G.Lind7,8 · NiklasBoers2,3,9 · DirkWitthaut10,11 · Leonardo RydinGorjão12,13 1 3 K. Riechers et al. 1 Introduction Paleoclimate proxy records provide evidence for past abrupt climate shifts from regional to at least hemispheric scale (e.g. Menviel et al. 2020; Brovkin et al. 2021; Boers et al. 2022). Long-term climate simulations suggest that anthropogenic global warming could trigger structurally similar transitions in several Earth system components in the future, i.e., that these components could ‘tip’ to a qualitatively different state (e.g. Lenton et al. 2008, 2019; Boers 2021; Armstrong et al. 2022; Boulton et al. 2022; Wang et al. 2023). Such catastrophic shifts would have severe consequences on societies and ecosystems and may even unleash feedbacks, further increasing the global mean temperature. However, the assessment of potentially upcoming tipping points is challenging as the capability of modern complex climate models to simulate climate tipping dynamics is still limited (Valdes 2011; Liu et al. 2017; Wang et al. 2023). In light of this, the study of past abrupt climate shifts may provide insights into the processes involved in climate tipping events. Furthermore, past events may serve as benchmarks for the performance of fully coupled models in simulating the non-linear and high-dimensional dynamics that could lead to tipping events. In this context, we reassess here two proxy time series from the NGRIP ice core (North Greenland Ice Core Projects members 2004), which feature pronounced imprints of abrupt climatic transitions, by means of the Kramers–Moyal equation. Agnostic time series models, i.e., models whose dynamics appear to reproduce nature but are not entirely based on physical mechanisms, have played a major role in furthering the debate on climate tipping phenomena (e.g. Riechers et al. 2023a; Boers et al. 2017; Mitsui and Crucifix 2017; Kwasniok 2013; Lohmann and Ditlevsen 2018; Dakos et al. 2008; Bochow and Boers 2023). The ability to produce quantitatively similar dynamical behaviour building only on heuristic physical assumptions facilitates the statistical analysis of tipping phenomena, employing methods of stochastic analysis (Lenton et al. 2012; Morr and Boers 2024; Morr et al. 2024). The common concept of a climate tipping element is that of a dynamical system whose current stable equilibrium state is prone to annihilation in a dynamic bifurcation (Scheffer et al. 2009; Ashwin et al. 2012; Boers et al. 2022). This typically involves the reduction of complex, high-dimensional dynamics to just a few (if not one) summary observables that may be modelled in terms of stochastic differential equations (SDEs), i.e., as random dynamical systems. Therein, the noise term reflects the action of the unresolved dynamics on the summary observable (Hasselmann 1976). A common choice is to force the resolved variables with Gaussian white noise, but this approach may be overly simplistic in many situations. In particular, in the context of climate tipping points, a deviation from Gaussian white noise has important implications for the detection of early warning signals and for the probability of premature noise-induced tipping (Ditlevsen 1999; Lucarini et al. 2022; Benson et al. 2024; Kuehn ey al. 2022; Morr and Boers 2024). Here, we investigate the famous heavy-oxygen δ18 O record from the NGRIP ice core (North Greenland Ice Core Projects members 2004). The data shows that repeated decadal-scale warming events of regionally up to 16 ◦ C in amplitude, known as Dansgaard–Oeschger events, punctuated the North Atlantic climate throughout the last glacial interval (Dansgaard et al. 1984; Broecker et al. 1985; Johnsen et al. 1992; Dansgaard et al. 1993; Kindler et al. 2014). The sudden temperature increases were followed by a phase of moderate cooling before the temperatures ultimately relaxed back to colder levels in a second phase of more abrupt cooling. The two distinct cold and mild regimes are termed stadials and interstadials, respectively. In line with the SDE approach outlined above, we regard the δ18 O and dust concentration records as realisations of one-dimensional Markov processes and estimate the corresponding KM coefficients (Tabar 2019). The two records exhibit concomitant shifts, which are interpreted as sudden adjustments of global mean temperature and reorganisations of the atmospheric circulation of at least hemispheric scale (Fuhrer et al. 1999; Ruth et al. 2003, 2007; Schüpbach et al. 2018). Previous studies have motivated through statistical means the employment of a Markovian framework for these dynamics Riechers et al. (2023b); Kwasniok (2013). There have also been conceptual arguments of time-scale separation that lend this framework credence Gottwald (2021); Riechers et al. (2024). Under this modelling assumption, the KM coefficients are closely related to the Fokker–Planck equation of time-evolving diffusive systems. We estimate from the data and subsequently compare the two Kramers–Moyal expansions with respect to their implied stochastic model structure. Specifically, we investigate whether the dynamics can each be represented by a canonical Langevin approach or whether a discontinuous noise component, such as Poisson jump diffusion, is needed. This article is structured as follows: In Sec. 2 we briefly introduce the two paleo-climatic proxies that we examine. Subsequently, in Sec. 3, we detail the Kramers–Moyal expansion in one dimension as the prime method to construct time series models including noise and possibly discontinuous elements. Section 4 presents the results of this analysis: Herein, we show the monoand bistability of the obtained models of the two records and discuss the need to choose a noise model different from Gaussian white noise. 1 3 465 Page 2 of 14 Discontinuous stochastic forcing in Greenland ice core data In Sec. 5 we discuss our findings and relate them to previous work. Sec. 6 summarises our key findings and draws conclusions. 2 Data and pre-processing This work relies on the δ18 O and dust concentration records obtained by the North Greenland Ice Core Project (NGRIP) (Ruth et al. 2003; North Greenland Ice Core Projects members 2004; Gkinis et al. 2014). From 1404.75 m to 2426.00 m of depth the joint record is provided at 5 cm equidistant resolution. This translates to the time span from 59945 yr to 10276 yr b2k (before 2000 CE) with ∼5 yr resolution for the oldest and sub-annual resolution for the most recent part of the record (Fig. 1a and b). For the analysis, the data was rescaled, binned to an equidistant time axis of 5-year resolution, detrended, and normalised (see Appendix A for details). The concentration of dust, i.e., the number of particles with a diameter larger than 1 µ m per ml, is commonly interpreted as a proxy for the state of the hemispheric atmospheric circulation (e.g. Fischer et al. 2007; Ruth et al. 2007; Schüpbach et al. 2018; Erhardt et al. 2019). In particular, the dust storm activity and dryness over East Asian desserts, the strength and position of the polar jet, and local precipitation patterns govern the emission, transport, and deposition of the dust, respectively (Fischer et al. 2007; Erhardt et al. 2019). Correspondingly, the substantial changes in the dust concentrations at DO events are interpreted as large-scale reorganisations of the Northern Hemisphere’s atmospheric circulation. In agreement with a widespread convention, we rescale the dust record by taking the net negative logarithm (e.g. Ditlevsen 1999; Mitsui and Crucifix 2017; Boers et al. 2017; Riechers et al. 2023a). In this form, the dust record exhibits a high degree of correlation with the δ18 O record (Boers et al. 2017). In order to reduce the influence of slow changes in the background climate, we restricted the analysis to the period 59–27 kyr b2k and applied further detrending with respect to a Northern Hemisphere temperature reconstruction provided by Snyder (2016) (see Fig. 1c and d and App. Appendix A). The concentration of stable water isotopes expressed as δ18 O values in units of permil is a proxy for the site temperature at the time of precipitation (Jouzel et al. 1997; Gkinis et al. 2014). Fig. 1 Trajectories of the 20-year mean of δ18 O (a) and accordingly resampled dust concentrations (b) from the NGRIP ice core in Greenland, from 122 kyr and 107 kyr to 10 kyr before 2000 CE (b2k), respectively (Ruth et al. 2003; Rasmussen et al. 2014; Seierstad et al. 2014). The dust data is given as the negative natural logarithm of the actual dust concentrations, in order to facilitate comparison to the δ18 O data. Panels (c) and (d) show the same proxies but at a higher resolution of 5 years (North Greenland Ice Core Projects members 2004; Gkinis et al. 2014; Ruth et al. 2003) and over the shorter period from 59 to 27 kyr b2k. The analysis presented in this study was constrained to this segment of the records. The two proxy time series in (c) and (d) have been detrended by linearly regressing the data against reconstructed global mean surface temperatures (Snyder 2016) and removing the apparent background-temperature-driven slow change. The grey shadings mark the Greenland interstadial (GI) intervals according to (Rasmussen et al. 2014). All data are shown on the GICC05 chronology (Vinther et al. 2006; Rasmussen et al. 2006; Andersen 2006; Svensson et al. 2008). The data were binned to equidistant time resolution from its original 5 cm depth resolution (see App. Appendix A for further details on the data processing (Riechers et al. 2023a)) 1 3 Page 3 of 14 465 K. Riechers et al. 1974), a Markov stochastic process of the form Eq. (1) generally features discontinuous paths with non-zero probability. Path-wise continuity is only one of many notions of continuity in stochastic processes. Another is the continuity criterion for Markov processes provided by Gardiner (2009), which requires for a process to be continuous that C (x, t, δ) = lim τ→0 1 τP(|xt+τ−xt|>δ) = lim τ→0 1 τˆ | x ′− x | >δ p(x′,t+τ|x, t)dx′! =0, (5) for all δ , x, and t. In words, this means that the probability of a particle deviating from a reference position more than δ in a time interval τ decreases faster than linearly with τ . The presence of higher-order KM coefficients in the corresponding Kramers–Moyal equation is a necessary, yet not sufficient criterion for a given process to be discontinuous under this latter notion. 3.1 Estimating Kramers–Moyal coefficients The central entry point for this work is Eq. (3). It provides a means to estimate the KM coefficients Dm(x) directly from data, i.e., from a recorded realisation of a stochastic process, provided that the following assumptions are fulfilled (to a reasonable degree): i) The observed process is a Markov process, ii) the process is time-homogeneous, i.e., the dynamics did not change over time, iii) the state space is sampled sufficiently densely, iv) and the sampling time is short compared to the characteristic time scale of the dynamics. Under these conditions, the evaluation of the conditional statistical moments M(x, τ) at the shortest available time lag ∆t given by the sampling rate yields a good estimate for the KM coefficients: ˆ D m(x)= 1 m! 1 ∆t⟨ (xt+∆t − xt)m | xt=x ⟩≈ Dm(x) , (6) wherein the ensemble average in Eq. (3) is replaced by the average over the available data ⟨·⟩ . Our numerical implementation of Eq. (6) is based on the Nadaraya–Watson estimator which is detailed in App. Appendix B. 3 Methods Our starting point is a (time-homogeneous) Markov stochastic process xt of the form dxt=f(xt)dt+σ(xt)dξt, (1) where dξt denotes an arbitrary uncorrelated stochastic force. The temporal evolution of the associated conditional probability function p(x, t+τ|x′ ,t) then follows the Kramers–Moyal equation (Kramers 1940; Moyal 1949; Kampen 1961; Gardiner 2009; Risken and Frank 1996; Tabar 2019): ∂ ∂τ p(x, t+τ|x′ ,t)= ∞ ∑ m=1 ( −∂ ∂x )m Dm(x)p(x, t+τ|x′ ,t) . (2) The Kramers–Moyal (KM) coefficients Dm(x) are related to the conditional moments Mm(x, τ) of order m of the stochastic variable x at a time-lag τ by D m(x)= 1 m!lim τ→0 1 τMm(x, τ) =1 m! lim τ→0 1 τˆ (x′ − xt)mp(x′,t+τ | x, t)dx′ . (3) In the special case that the stochastic force in Eq. (1) is given by Gaussian white noise (i.e., it can be expressed by the increments of a Wiener process Wt ), only the first two terms on the right of Eq. (2) contribute and the Kramers–Moyal equation reduces to the better-known Fokker–Planck equation (Fokker 1913, 1914; Planck 1917). With dξt=dWt , Eq. (1) becomes the Langevin equation and the resulting process is then referred to as a Langevin process1. For Langevin processes the relation D 1(x)=f(x)and D2(x)= 1 2 σ2(x) , (4) between the KM coefficients, the drift f(x) and the diffusion σ(x) , holds in general. The other way around, if higher-order moments contribute to the Kramers–Moyal equation, the underlying process cannot be a standard Langevin process. In that case, ξt does not correspond to a Wiener process but has instead a more complex form. However, the first two KM coefficients would still be dominated by the process’ drift and diffusion. While a Langevin process consists, with probability 1, of continuous sample paths (e.g. Theorem 5.1.1 in Arnold 1 There is no agreement on the use of the term Langevin process. Some authors consider Lévy-driven equations as such Langevin equations, others prefer to refer to Langevin processes as those that are solely driven by Gaussian/Brownian noise. 1 3 465 Page 4 of 14 Discontinuous stochastic forcing in Greenland ice core data Langevin from jump-diffusion processes, namely the Θ -ratio Θ( x, τ)= 3 M2 ( x, τ )2 M4(x, τ)∼{1 , Langevin , 1 τ ,jump-diffusion, (9) and the Q-ratio (Lehnertz et al. 2018) Q (x, τ)= M 6 (x, τ) 5 M 4( x, τ )∼{ τ, Langevin , constant,jump-diffusion. (10) For details on the derivation of these relationships, we refer the interested reader to (Tabar 2019). Observing either of the scalings given in Eqs. (9) and (10), respectively, can aid in deciding between employing a Langevin or jump-diffusion model. These relationships are specifically derived for the jumpdiffusion model. Compared to the ratio of D4 and D2 discussed above, the results here are more prone to invalidity due to unjustified modelling assumptions on the real data. For different noise models than the ones introduced above, different scaling behaviours of these ratios with respect to τ will arise. Data from, e.g., a non-Markovian system may, on the other hand, exhibit the described behaviour while actually harbouring entirely different internal dynamics. In this work, we focus on distinguishing between the Langevin and Poisson jump-diffusion models as two archetypical (dis-) continuous stochastic models. Observing any other scaling in Q or Θ may hint at a third model being more appropriate to reproduce the time series dynamics. However, in the context of continuous versus discontinuous stochastic models, considering the two discussed models yields essential information. 4 Results Figure 2 shows the first and second KM coefficients, and the ratio of the second to the fourth KM coefficients, as estimated from the dust and δ18 O time series according to Eq. (6). The corresponding Θ and Q ratios are presented in Fig. 3. a Dust record: For the dust, the constructed drift D1(x) in Fig. 2b exhibits two separate stable states that match the maxima of the probability density function in Fig. 2a. The second KM coefficient D2(x) in Fig. 2c is approximately constant. The ratio between the fourth and the second KM coefficients in Fig. 2d is smaller than 0.1 on the entire state space probed by the time series. For large portions of the dust’s state space, we find in Fig. 3 a decrease of the Θ(x, τ) ratio with increasing τ , similar to a 1/τ behaviour. This applies, in particular, at the stable equilibria of the drift, 3.2 Estimators of discontinuous motion Once the KM coefficients are estimated from the data, one can draw inference on the most fitting choice of the noise model dξt . Vanishing higher-order moments ( m>2 ) classify the model as a Langevin process. In contrast, demonstrable contributions of these moments suggest that the process is best modelled by including noise beyond a Wiener process (see e.g. Kampen 1961; Van Kampen 2007; Gardiner 2009; Tabar 2019; Lin 2023). The finite sampling time step ∆t introduces a bias for the estimators ˆ Dm( x ) (Kurth et al. 2021). As a consequence, even for a Langevin process the expected values for the higher-order KM estimators differ from zero. A first pragmatic metric to discern whether a studied process is a Langevin process or not is to evaluate the ratio between the fourth KM coefficient and the second, i.e., D4(x)/D2(x) . This gauges the distributional tail of all immediate disturbances originating from x. It therefore offers a non-parametric insight into whether a fat tail of disturbances is needed to recreate the dynamics at the considered sampling rate. Such conclusions would be largely model-independent and do not explicitly rely on the Markovianity of the data. Small values ≲0.1 are typically regarded as a justification for a Langevin description. Values D4(x)/D2(x)≳0.1 point to non-diffusive motion (i.e., forcing beyond Gaussian white noise). This metric offers a first insight into whether a discontinuous noise term ξt is needed to model the process (Gao et al. 2016; Lu and Duan 2020; Lucarini et al. 2022). When the Langevin process model is contrasted with a jump-diffusion model of the form (Tabar 2019; Lin 2023) d x t= f ( x t)d t + σ ( x t)d W t+ η ( x t)d J (λ) t, (7) the assessment can be further refined. Here, J( λ ) t denotes a Poissonian jump process characterised by the rate λ . The jump amplitude is determined by the Gaussian stochastic variable η(x) . For this specific process model, the KM coefficients read (Tabar 2019) D 1 (x)=f(x), D2(x)=1 2σ(x)2+1 2λ(x)⟨η(x)2⟩, D m(x)= 1 m! λ(x) ⟨ η(x)m ⟩ ,for m>2 , (8) where ⟨·⟩ expresses the expected value. Similarly, the bias of the KM estimators defined by Eq. (6), when applied to a jump-diffusion process sampled at finite time step ∆t , can be derived analytically. These considerations offer two additional metrics to distinguish 1 3 Page 5 of 14 465 K. Riechers et al. variable. Finally, the ratio D4(x)/D2(x)≳0.3 is 10 times larger for δ18 O than for the dust. The δ18 O record exhibits a mostly constant Θ(x, τ) -ratio with respect to τ , as seen in Fig. 3. It is slightly below but still close to 1 for large parts of the state space. The corresponding Q(x, τ) -ratio is likewise constant ( ≈1 ) with respect to τ , with variations in both directions. 5 Discussion The assessment of the KM coefficients and the scaling of Θ and Q ratios from the dust and the δ18 O records provides some insight into how to best model the proxy time series within the framework of one-dimensional stochastic processes. For the dust, we find bistability of the estimated model’s drift. The small D4(x)/D2(x) ratio and the linear increase of the Q(x, τ) with increasing τ indicate that this process can, in fact, be modelled as a Langevin process. Only the assessment of the Θ(x, τ) ratio calls this conclusion into question. For a Langevin process, this ratio should be equal to one, but we observe a 1/τ -like scaling for small values of τ , in line with an underlying jump-diffusion process. where the data availability is the best and our estimation is most robust. The dust Θ(x, τ) -ratio is close to 1 only in a region of its state space where its probability density has a local minimum ( −0.3≲dust ≲0.3 ). The corresponding Q(x, τ) -ratio shows quite a distinct linear increase with increasing τ – at least for small values of τ . For larger values of τ , Q(x, τ) is constant. b δ18O record: In the case of δ18 O, the drift has only one zero-crossing. This seems to explain the unimodal distribution of the data, though this broader distribution could also be caused by larger observational noise in the record. The mono-stability of the drift would not be affected by timeand state-independent observational noise and can therefore be seen as a more direct insight into the potential underlying dynamics. We note that NGRIP data products that provide δ18 O concentrations at a lower time resolution of 20or 50-year time steps exhibit a bimodal distribution. For the purposes of our analysis, however, only the highest available sampling rate of time series data should be used so as to curtail the biases incurred in the KM estimations. With respect to the normalised units, the first and second KM coefficients of δ18 O exceed their counterparts for dust by factors of approximately 4 and 10, respectively. This indicates that δ18 O was subjected to stronger noise while simultaneously stronger deterministic forces acted on the Fig. 2 The probability density function (PDF) of (a) the dust and (e) δ18 O. The non-parametric estimates of the (b, f) first KM coefficient D1(x) and (c, g) the second KM coefficient D2(x) . The ratio between the fourth and the second KM coefficient D4(x)/D2(x) (d and h). All KM coefficients are evaluated at the shortest available time step ∆t=5 yr of the time series. The estimated dust drift is bistable, while that of δ18 O is monostable. The second KM coefficient D2(x) is relatively constant for both records. The ratio D4(x)/D2(x) is small ( ≲0.1 ) for the dust record. Yet, it is non-negligible for δ18 O ( ≳0.3 ) in large parts of the state space, suggesting that the driving noise in a stochastic model for these time series should not be exclusively Gaussian white noise. Details on the choice of kernel and bandwidth used for the KM coefficient estimation, as well as an analysis of the influence of the kernel bandwidth, can be found in App. B 1 3 465 Page 6 of 14 Discontinuous stochastic forcing in Greenland ice core data with two regimes in the presence of a single equilibrium and time asymmetry (Chechkin et al. 2003, 2004; Metzler and Klafter 2004; Yang et al. 2020). Given the high degree of visual similarity between the dust and the δ18 O records, the differences in the reconstructed potentials and the ratio between the fourth and the second KM coefficient are remarkable. This accentuates the need for careful statistical analysis when devising time series models for non-linear systems with abrupt transitions. Adopting a generalised Langevin equation with a bistable drift term, Ditlevsen (Ditlevsen 1999) showed that the noise in the calcium concentration record from the GRIP ice core can be modelled with an α -stable component. Calcium concentrations are typically considered equivalent to dust concentrations (cf. (Fuhrer et al. 1993; Ruth et al. 2002, 2003; Fischer et al. 2007)). We cannot directly assess the presence of α -stable noise in the NGRIP dust record. This is because noise models with infinite statistical moments, which can We note that a simple Langevin model with a bistable drift and purely diffusive noise can produce the regime shifts observed in the dust record. However, such a model is unlikely to reproduce the asymmetric shape of the interstadial phases evident in the record. For the δ18 O record, the results are exactly the opposite. The constructed drift function exhibits only a single stable equilibrium. The observed quantities D4(x)/D2(x) and Q(x, τ) provide evidence for relevant contributions from higher-order KM coefficients. The Θ(x, τ) ratio, however, is close to one in agreement with a Langevin model. A Langevin model together with the evidenced single equilibrium of the drift function clearly fails to explain the two regimes of the δ18 O record, and the apparent time asymmetry. Taken together, we conclude that the evidence speaks in favour of introducing discontinuities to the driving noise model rather than against it. Complex noise, i.e., noise beyond a Wiener process, could indeed be a way to reproduce time series Fig. 3 The Θ(x, τ) -ratio and Q(x, τ) -ratio of the dust and the δ18 O concentration. (a) The Θ(x, τ) -ratio of the dust is not close to 1 for a large range of τ and x, particularly around the peaks of the bimodal dust distribution. (b) From τ=5 to roughly τ= 50 the Q(x, τ) -ratio increases linearly with τ , consistent with a continuous process, yet for τ>50 the Q(x, τ) -ratio is nearly constant, consistent with a discontinuous process. (c) On the one hand, the Θ(x, τ) -ratio of the δ18 O points to being different from 1 along the peak of the distribution, yet not sufficiently conclusive to ascertain if the δ18 O is discontinuous. (d) On the other hand, the Q(x, τ) -ratio is arguably constant over τ , consistent with a discontinuous δ18 O. For visualisation purposes, Q(x, τ) of δ18 O is multiplied by 0.6 to match the scale of Q(x, τ) of the dust (Harris et al. 2020; Virtanen et al. 2020; Rydin Gorjão and Meirinhos 2019; Rydin Gorjão et al. 2023; Hunter 2007) 1 3 Page 7 of 14 465 K. Riechers et al. stability of the actual physical processes. It is, however, notable that these findings are inconsistent with the hypothesis that past Greenland temperatures were governed by intrinsically bistable dynamics (Livina et al. 2010; Kwasniok 2013). For the atmospheric circulation, there is no such inconsistency. We stress that this inconsistency with respect to previous studies arises under the application of differing modelling assumptions and data pre-processing procedures. The potential influence of non-Markovianity in the dynamics or complex measurement noise cannot be quantified. Disentangling these confounding effects is possible but demands large amounts of data Böttcher et al. (2006). We nonetheless maintain that this novel perspective is a valuable data point for further conceptual and physical considerations of DO events. We found that stochastic forcing should include terms beyond Gaussian white noise when modelling the δ18 O record. This renders the Langevin approach insufficient to accurately reproduce the time series characteristics, drawing attention towards including discontinuous elements. For the dust record, similar indications could be found, though these have not been as convincing. In physical terms, complex noise could have played a central role in the emergence of DO events. Our analysis does not provide direct evidence for a causal relation between discontinuous driving noise and the regime switches of the North Atlantic region’s climate during the last glacial. Yet, it motivates further exploration of this issue along the lines of Gottwald (2020) and Riechers et al. (2023a). The possibility that climate tipping elements are subject to non-Gaussian noise in today’s climate should receive greater consideration. The corresponding implications on the stability of these elements and the ability to detect early warning signals should be investigated. Appendix A: Data detrending As mentioned in Sec. 2, this study focuses on the period 59–27 kyr b2k. Detrending of the data is necessary to ensure that the time series are time-homogeneous stationary processes, which is an underlying assumption for the Kramers–Moyal analysis performed in our investigation. To compensate for the influence of the background climate on the climate proxy records of dust and δ18 O, we remove a linear drift with respect to reconstructed global average surface temperatures (Snyder 2016) from both time series. Figure 4 illustrates the detrending scheme. Due to the two-regime nature of the time series, a simple linear regression would overestimate the temperature dependencies. Instead, we separate the data from Greenland stadials (GS) and Greenland interstadials (GI) and then minimise the expression be found in α -stable distributions, are inherently incompatible with the Kramers–Moyal framework. Yet, our results corroborate the notion that Greenland ice core records bear the signature of non-Gaussian noise, though in our analysis this arises primarily for the δ18 O record. Related to this, Gottwald (2020) recently formulated a conceptual model of DO events wherein α -stable noise plays a central role as an event trigger, later extended by Riechers et al. (2023a). From the perspective of theoretical stochastic modelling, it is worth noting that the α -stable noise model leads to a pathwise continuous process, in contrast to the Poisson jumpdiffusion model discussed in this work. Employing the continuity notion of Eq. (5), however, both of these models would be considered discontinuous. We have to state that the interpretation of higher-order KM coefficients is not straightforward and depends on the exact choice of the stochastic model. A direct causal relation between the DO events and discontinuous noise cannot be inferred without further ado within this study, but the role of discontinuities in the proxy records merits further investigation. It has been observed in complex model simulations that (stochastic) atmospheric anomalies can indeed drive regime changes in the North Atlantic region (Drijfhout et al. 2013; Kleppin et al. 2015). Together with the apparent aptitude of non-Gaussian noise models for Greenland temperature and Northern Hemisphere atmospheric circulation proxies, this motivates further research on the effect that non-Gaussian noise could have on climate tipping elements in present-day climate. If both Greenland temperatures and the state of the Northern Hemisphere atmospheric circulation were subject to non-Gaussian noise, and if indeed pulses of this noise triggered transitions between stadial and interstadial regimes, this would have important implications for our conception of stability of certain climate tipping elements. The possibility that climate tipping elements are nowadays likewise subject to non-Gaussian stochastic forcing warrants more attention. 6 Conclusion In this work, we presented a data-driven analysis of the δ18 O and dust concentration records from the NGRIP ice core, based on the Kramers–Moyal equation. This equation generalises the Fokker–Planck equation by allowing for arbitrarily complex uncorrelated driving noise dξt . In particular, such noise may result in a discontinuous process. The estimation of the KM coefficients yielded a monostable drift for the isolated δ18 O record and a bistable one for the dust. The analysis of the resulting agnostic time series models does not allow for conclusions about the dynamical 1 3 465 Page 8 of 14 Discontinuous stochastic forcing in Greenland ice core data D m(x)∼ 1 m! 1 ∆t⟨(xt+∆t−xt)m|xt=x⟩ ∼1 m! 1 ∆t 1 N N−1 ∑ i=1 K(x−xi)(xi+1 −xi)m . (B1) Similarly to selecting the number of bins in a histogram, for a kernel-density estimation, we select both a kernel and a bandwidth (Nadaraya 1964; Watson 1964; Lamouroux and Lehnertz 2009). The kernel is a function K(x) for the estimator  fh(x) , where h is the bandwidth at a point x, following  f h(x)= 1 nh n  i=1 K  x−xi h  (B2) for a collection {xi} of n random variables. The kernel K(x) is normalisable ´K(x)dx=1 and has a bandwidth h, such that K(x)=K(x/h)/h (Rydin Gorjão et al. 2019; Tabar 2019; Davis and Buffett 2022). The bandwidth h is equivalent to the selection of the number of bins, except that binning in a histogram is always ‘placing numbers into nonoverlapping boxes’. We use an Epanechnikov kernel R 2= N ∑ i=1 ( Xti−aX∆T(ti)− { bGI,if ti∈GI bGS,if ti ∈ GS )2 , (A1) with X either dust or δ18 O and with respect to the parameters aX , bGI , and bGS . For a given time ti∈GS (GI) indicates that ti falls into a stadial (interstadial) period. The resulting aX is used to detrend the original data with respect to the temperature. The detrended data are subsequently normalised by subtraction of their mean and division by their standard deviation. Appendix B: Nadaraya–Watson estimator of the KM coefficients and bandwidth selection In order to carry out the estimation in Eq. (6) we map each data point in the corresponding state space to a kernel density and then take a weighted average over all data points Fig. 4 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 (light brown) as well as detrended and normalised (purple) time series are shown. Likewise for the dust record in panel (b) (light brown and green, respectively). The background temperature is given in anomalies with respect to the mean over the investigated period (blue). Panels (c) and (d) show scatter plots of 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 lines correspond to the fitting scheme that uses a single slope but two different offsets to separately fit the stadial and interstadial data 1 3 Page 9 of 14 465