Stable stadial and interstadial states of the last glacial's climate identified in a combined stable water isotope and dust record from Greenland

During the last glacial interval, the Northern Hemisphere climate was punctuated by a series of abrupt changes between two characteristic climate regimes. The existence of stadial (cold) and interstadial (milder) periods is typically attributed to a hypothesised bistability in the glacial North Atlantic climate system, allowing for rapid transitions from the stadial to the interstadial state – the so-called Dansgaard–Oeschger (DO) events – and more gradual yet still fairly abrupt reverse shifts. The physical mechanisms driving these regime transitions remain debated. DO events are characterised by substantial warming over Greenland and a reorganisation of the Northern Hemisphere atmospheric circulation, which are evident from concomitant shifts in the δ¹⁸O ratios and dust concentration records from Greenland ice cores. Treating the combined δ¹⁸O and dust record obtained by the North Greenland Ice Core Project (NGRIP) as a realisation of a two-dimensional, time-homogeneous, and Markovian stochastic process, we present a reconstruction of its underlying deterministic drift based on the leading-order terms of the Kramers–Moyal equation. The analysis reveals two basins of attraction in the two-dimensional state space that can be identified with the stadial and interstadial regimes. The drift term of the dust exhibits a double-fold bifurcation structure, while – in contrast to prevailing assumptions – the δ¹⁸O component of the drift is clearly mono-stable. This suggests that the last glacial's Greenland temperatures should not be regarded as an intrinsically bistable climate variable. Instead, the two-regime nature of the δ¹⁸O record is apparently inherited from a coupling to another bistable climate process. In contrast, the bistability evidenced in the dust drift points to the presence of two stable circulation regimes of the last glacial's Northern Hemisphere atmosphere.

1 Introduction

Recently, evidence was reported for the destabilisation of climatic subsystems likely caused by continued anthropogenically driven climate change  (e.g. Boers, 2021; Rosier et al., 2021; Boers and Rypdal, 2021). Conceptually, such destabilisation is commonly formulated in terms of bistable dynamical systems that approach a bifurcation in response to the gradual change in a control parameter. This setting offers three mechanisms for the system to transition between two alternative stable states (Ashwin et al., 2012). First, the control parameter may cross a bifurcation, which dissolves the currently attracting state and necessarily entails a transition to the remaining alternative stable state (bifurcation-induced transition). Second, random perturbations may push the system across a basin boundary (noise-induced transition); this is generally more likely the closer the system is to a bifurcation. Third, rapid change in the control parameter may shift the basin boundaries at a rate too high for the system to track the moving domain of its current attractor (rate-induced transition). If global warming – viewed as the control parameter – were to exceed certain thresholds, several elements of the climate system are thought to be at risk to “tip” to alternative stable states (Lenton and Schellnhuber, 2007; Lenton et al., 2008; Armstrong McKay et al., 2022), among them the Greenland ice sheet (Boers and Rypdal, 2021), the Amazon rainforest (Boulton et al., 2022), the Atlantic Meridional Overturning Circulation (AMOC) (Boulton et al., 2014; Boers, 2021), and the West Antarctic ice sheet (Rosier et al., 2021).

The possibility of alternative stable states of the entire climate system or its subsystems (and transitions between these) has been discussed at least since the 1960s (e.g. Ghil, 1975; North, 1975; Stommel, 1961). Empirical evidence, however, that the climate system or its subsystems can indeed abruptly transition between alternative equilibria is available only from proxy records which allow reconstruction of past climatic conditions prior to the instrumental period (e.g. Brovkin et al., 2021; Boers et al., 2022; and references therein). Given that comprehensive Earth system models continue to have problems in simulating abrupt climate changes and especially in reproducing abrupt changes evidenced in proxy records (Valdes, 2011), studying abrupt changes recorded by paleoclimate proxies is key for gaining a better understanding of the physical mechanisms involved and for assessing the risks of future abrupt transitions.

In this context, our study investigates the Dansgaard–Oeschger (DO) events, a series of abrupt warming events over Greenland first evidenced in stable water isotope records from Greenland ice cores (Dansgaard et al., 1982, 1984, 1993; Johnsen et al., 1992; North Greenland Ice Core Projects members, 2004). While locally the temperature increases are estimated to be as large as 16 ^∘C in the annual mean temperature (Kindler et al., 2014), a (weaker) signature of these events can be found in numerous records across the globe (e.g. Voelker, 2002; Menviel et al., 2020; and references therein) indicating changes in other climatic subsystems such as Antarctic average temperatures (e.g. WAIS Divide Project Members, 2015; EPICA Community Members, 2006), the Asian and South American Monsoon system (e.g. Wang et al., 2001; Kanner et al., 2012; Cheng et al., 2013; Li et al., 2017; Corrick et al., 2020), or the Atlantic Meridional Overturning Circulation (AMOC) (e.g. Lynch-Stieglitz, 2017; Henry et al., 2016; Gottschalk et al., 2015). The global puzzle of more or less abrupt shifts in synchrony (within the limits of dating uncertainties) with DO events found in versatile paleoclimate proxy records point to a complex scheme of interactions between climatic subsystems involved in the DO variability that dominated the last glacial period. While multiple lines of evidence indicate a central role of changes in the overturning strength of the AMOC (e.g. Lynch-Stieglitz, 2017; Menviel et al., 2020), to date, there is no consensus about the ultimate trigger of DO events.

An important branch of research has assessed the performance of low-dimensional conceptual models in explaining the DO variability in the Greenland ice core records (e.g. Ditlevsen, 1999; Livina et al., 2010; Kwasniok, 2013; Mitsui and Crucifix, 2017; Roberts and Saha, 2017; Boers et al., 2017, 2018; Lohmann and Ditlevsen, 2018a; Vettoretti et al., 2022). Typically, one-dimensional multi- or bistable models (Ditlevsen, 1999; Livina et al., 2010; Kwasniok, 2013; Lohmann and Ditlevsen, 2018a) or two-dimensional relaxation oscillators (Kwasniok, 2013; Mitsui and Crucifix, 2017; Roberts and Saha, 2017; Lohmann and Ditlevsen, 2018a; Vettoretti et al., 2022) have been invoked, forced by either slowly changing climate background variables such as CO₂ or changing orbital parameters, by noise, or by both. To the best of our knowledge, to date Boers et al. (2017) presented the only inverse-modelling approach to simulate a two-dimensional Greenland ice core proxy record – δ¹⁸O and dust – with regards to its DO variability. Likewise in two dimensions, we present here a data-driven investigation of the couplings between Greenland temperatures and the larger-scale Northern Hemisphere state of the atmosphere represented by the North Greenland Ice Core Project (NGRIP) δ¹⁸O ratio and dust concentration records, respectively (North Greenland Ice Core Projects members, 2004; Gkinis et al., 2014; Ruth et al., 2003). Treating the combined δ¹⁸O and dust record as the realisation of a time-homogeneous Markovian stochastic process (Kondrashov et al., 2005, 2015), we reconstruct the corresponding deterministic two-dimensional drift using the Kramers–Moyal equation (Kramers, 1940; Moyal, 1949; Tabar, 2019) and reveal evidence for bistability of the coupled δ¹⁸O–dust “system”. Compared to the previously mentioned studies, this approach has the advantage that the estimation of the drift is non-parametric (i.e. it assumes no a priori functional structure for the drift) and that it assesses the stability configuration of the two-dimensional record as opposed to the numerous studies concerned with one-dimensional proxy records.

In the state space spanned by δ¹⁸O ratios and dust concentrations, based on our results we identify two regions of convergence concentrated around two stable fixed points, which can be associated with Greenland stadials and interstadials. We show that the global bistability is rooted in the dust component of the drift, exhibiting what seems to be a double-fold bifurcation parameterised by δ¹⁸O. This asserts a genuine bistability to the glacial Northern Hemisphere atmosphere. In contrast, the δ¹⁸O drift component is mono-stable across all dust values, suggesting that the two regimes evidenced in past Greenland temperature reconstructions are not the signature of intrinsic bistability but that of coupling another bistable subsystem, which – according to our results – may be the atmospheric large-scale circulation.

This article is structured as follows. We first present the paleoclimate proxies analysed in this study and explain how we pre-processed the data to make them suitable for estimating the two-dimensional drift (Sect. 2). Subsequently, we introduce the two-dimensional Kramers–Moyal equation, which is key for the analysis (Sect. 3). Section 4 provides the reconstruction of the two-dimensional drift, and Sect. 5 discusses the results and how they relate to previous studies. In Sect. 6 we summarise our main findings and detail the research questions that follow from these.

Figure 1 (a) Global-mean surface temperature reconstruction for the last glacial interval as provided by Snyder (2016) and linearly interpolated to a 20-year temporal resolution. The reconstruction is based on a multi-proxy database which comprises over 20 000 sea surface temperature reconstructions from 59 marine sediment cores. The figure shows the anomaly with respect to modern climate (5–0 kyr b2k average; b2k: before 2000 CE). The 20-year mean of δ¹⁸O ratios (b) and accordingly resampled dust concentrations (c) from the NGRIP ice core in Greenland, from 122 kyr and 107 kyr to 10 kyr b2k (Rasmussen et al., 2014; Seierstad et al., 2014; Ruth et al., 2003). The dust data are given as the negative natural logarithm of the actual dust concentrations in order to facilitate visual comparison to the δ¹⁸O data. Panels (d) and (e) show the linear regressions of δ¹⁸O and dust onto the reconstructed global-mean surface temperatures (Snyder, 2016) from (a), carried out separately for Greenland stadials (GSs) and Greenland interstadials (GIs). Panels (f) and (g) show the same proxies as shown in (b) and (c) 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 a shorter period from 59 to 27 kyr b2k. The analysis presented in this study was constrained to this section of the record. The two proxy time series in (f) and (g) have been detrended by removing the slow non-linear change induced by changes in the global background temperatures, based on the regressions from (d) and (e). The data were then binned to equidistant time resolution from the original 5 cm depth resolution. The grey shadings mark the GI intervals according to Rasmussen et al. (2014). Panel (h) shows the histograms of the two time series shown in (f) and (g), respectively. All data are shown on the GICC05 chronology (Vinther et al., 2006; Rasmussen et al., 2006; Andersen et al., 2006; Svensson et al., 2008).

2 Data and pre-processing

The analysis presented here is based primarily on the joint δ¹⁸O ratio and dust concentration time series obtained by the North Greenland Ice Core Project (NGRIP) (North Greenland Ice Core Projects members, 2004; Ruth et al., 2003; Gkinis et al., 2014). From 1404.75–2426.00 m of depth in the NGRIP ice core, data are available for both proxies at a spatially equidistant resolution of 5 cm. This translates into non-equidistant temporal resolution ranging from sub-annual resolution at the end to ∼ 5 years at the beginning of the period, 59944.5–10276.4 years b2k, according to the Greenland Ice Core Chronology 2005 (GICC05), the common age–depth model for both proxies (Vinther et al., 2006; Rasmussen et al., 2006; Andersen et al., 2006; Svensson et al., 2008). Lower-resolution data (20-year means) (Rasmussen et al., 2014; Seierstad et al., 2014) reaching back to the last interglacial period (see Fig. 1) are only used for illustrative purposes but not for the analysis.

The ratio of stable water isotopes, expressed as δ¹⁸O values in units of permil, is a proxy for the site temperature at the time of precipitation, and hence the abrupt shifts present in the data qualitatively indicate the abrupt warming events over Greenland (Jouzel et al., 1997; Johnsen et al., 2001). The concentration of dust, i.e. the number of particles with a diameter above 1 µm 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). More specifically, it is assumed to be controlled mostly by three factors (Fischer et al., 2007): first, by climatic conditions at the emission source, i.e. the dust storm activity over East Asian deserts preconditioned on generally dry regional climate; second, by the transport efficiency, which is affected by the strength and position of the polar jet stream; and third, the depositional process, which is mostly determined by local precipitation patterns. Correspondingly, the substantial changes in the dust concentrations across DO events are interpreted as large-scale reorganisations of the Northern Hemisphere's atmospheric circulation. Typically, atmospheric changes affecting the dust flux onto the Greenland ice sheet are accompanied by changes in the snow accumulation of opposite sign (e.g. Fischer et al., 2007). This enhances the corresponding change in the recorded dust particle concentration. However, for high-accumulation Greenland ice cores – such as NGRIP – the dust concentration changes still serve as a reliable indicator of atmospheric changes according to Fischer et al. (2007). Since the dust concentrations approximately follow an exponential distribution, we consider the negative natural logarithm of the dust concentration in order to emphasise the similarity to the δ¹⁸O time series. For ease of notation, we always use the term dust (or dust concentrations), although technically we refer to its negative natural logarithm.

In Fig. 1 we show the original low-resolution (b and c) and the pre-processed high-resolution data (f and g) together with corresponding histograms (h), also given in Fig. 3a. Clearly, two regimes can be visually distinguished: Greenland stadials are characterised by low δ¹⁸O ratios and high dust concentrations. Greenland interstadials (grey shading in panels f and g of Fig. 1) in general exhibit the reversed configuration besides a mild trend toward stadial conditions, which can be more or less pronounced during the individual interstadials. In our study, we use the categorisation of the climatic periods as presented by Rasmussen et al. (2014). The two-regime character of the time series translates into a bi-modal histogram of the dust data, as seen in Fig. 1h. In the case of the δ¹⁸O data, the stronger trend during interstadials and the higher relative noise amplitude mask a potential bi-modality, and the histogram appears unimodal. Notice that the somewhat counterintuitive combination of meta-stable distinct dynamical regimes and unimodal distributions of the associated variables has also been discussed in the context of atmospheric dynamics (Majda et al., 2006).

The analysis conducted in this work relies on the following assumptions and technical conditions:

• i.

  The data-generating process is sufficiently time-homogeneous over the considered time period.

• ii.

  The process is Markovian at the sampled temporal resolution.

• iii.

  The data are equidistant in time.

• iv.

  The relevant region of the state space is sampled sufficiently densely by the available data.

With regard to (i), a low-frequency influence of the background climate on the proxy values and on the frequency of DO events is evident (see Fig. 1), with suppressed DO variability during the coldest parts of the glacial and longer interstadials for its warmer parts (e.g. Rial and Saha, 2011; Roberts and Saha, 2017; Mitsui and Crucifix, 2017; Lohmann and Ditlevsen, 2018b; Boers et al., 2017, 2018). We therefore restrict our analysis to the period 59–27 kyr b2k, which is characterised by a fairly stable background climate and persistent co-variability between dust and δ¹⁸O (Boers et al., 2017). To remove the remaining influence of the background climate on the climate proxy records we remove a trend that is non-linear in time from both time series. This trend is obtained by linearly regressing the proxy data against reconstructed global-average surface temperatures (Snyder, 2016). Figure 1 illustrates the detrending scheme; due to the two-regime nature of the time series, a simple linear regression of the proxy variables onto the global-average surface temperatures would overestimate the temperature dependencies. Instead, we separate the data from Greenland stadials and Greenland interstadials and then minimise the quantity

$$\begin{array}{}\text{(1)}& R={\left(\sum _{i=\mathrm{1}}^{N}x\left(t_{\mathrm{i}}\right)-a\mathrm{\Delta }T\left(t_{\mathrm{i}}\right)-\left\{\begin{array}{l}{b}_{\mathrm{GI}},\mathrm{if}{t}_{\mathrm{i}}\in \mathrm{GI}\\ b_{\mathrm{GS}},\mathrm{if}{t}_{\mathrm{i}}\in \mathrm{GS}\end{array}\right.\right)}^{\mathrm{2}}\end{array}$$

once for x=δ¹⁸O and once for x, taken as dust concentrations. The optimisations yield optimal values for the parameters a, b_GI, and b_GS for dust and δ¹⁸O (see Fig. 1 panels d and e for δ¹⁸O and dust concentration, respectively). For a given time t_i we write t_i∈GS (GI) to indicate that t_i falls into a stadial (interstadial) period. The index i runs over all data points, and N denotes the total number of data points. The resulting slope a is used to detrend the original data with respect to the time-dependent background temperature:

$$\begin{array}{}\text{(2)}& x_{\mathrm{detrended}}\left(t_{\mathrm{i}}\right)=x\left(t_{\mathrm{i}}\right)-a\mathrm{\Delta }T\left(t_{\mathrm{i}}\right).\end{array}$$

Subsequently, the detrended data are normalised by subtracting their respective means and dividing by their respective standard deviations. After the detrending, all stadial (or interstadial) periods exhibit almost the same level of values, which allows the data to be considered as the outcome of a time-homogeneous and long-term stationary process (compare Fig. 1f and g). As a consequence of the regime transitions the process is certainly not stationary on short timescales but only on timescales larger than typical DO cycle periods. Levelling out the differences between the recurring climate periods guarantees a sufficiently dense sampling of the relevant region of the state space (iv) and prevents a blurring of the drift reconstruction (i).

Table 1 Unit root test of the detrended data. ADF refers to the augmented Dickey–Fuller test; ADF-GLS refers to the augmented Dickey–Fuller–GLS test. We reject the presence of a unit root in each of the time series at p<0.05.

Stationarity tests provide further confirmation that the detrended data are free of any slow underlying trends: we have applied two separate tests to assess the stationarity of the detrended data on timescales beyond single DO cycles. These tests are the augmented Dickey–Fuller test (ADF) and the augmented Dickey–Fuller–GLS test (ADF-GLS). Both tests test the possibility of a unit root in the time series (null hypothesis). The alternative hypothesis is that the time series does not have a unit root; i.e. it is stationary. We can safely reject the presence of a unit root in each time series at p<0.05 (see Table 1).

There is a trade-off between conditions (i) and (iv) concerning the choice of the data window. While an even shorter window would assure time homogeneity of the dynamics with higher confidence, the sampling of the state space would become insufficiently sparse. The above choice (59–27 kyr b2k) guarantees a sufficient number of recurrences of the pre-processed two-dimensional trajectory to the relevant state space regions to perform statistical analysis. To obtain time-equidistant records (iii), the data are binned into temporally equidistant increments of 5 years.

Figure 2 Autocorrelation ρ(τ) of the increments Δx_t of δ¹⁸O and dust records. Both records show a weak anti-correlation at the shortest lag τ=5 years and no correlation for τ>5 years. We thus consider the data Markovian.

The question of Markovianity (ii) is the most difficult to answer unambiguously. Here we draw on the following heuristic argument: the autocorrelation functions of the increments of both proxies shown in Fig. 2 exhibit weak anti-correlation at a shift of one time step, while correlations beyond this are negligible. Such a small level of correlation certainly rules out the notion that long-term memory effects played a major role in the emergence of the given time series. Bear in mind that this is a necessary yet not sufficient criterion to consider the data Markovian. For practical reasons, we refrained from further Markovianity tests.

Finally, further preconditions for our endeavour are the fact that the NGRIP record exhibits an exceptionally high resolution (iv) compared to other paleoclimate archives and that the two time series share the same time axis.

3 Methods

In this work, we treat the combined δ¹⁸O and dust record as a trajectory of a two-dimensional, time-homogeneous, and Markovian stochastic process of the form

$$\begin{array}{}\text{(3)}& \mathrm{d}\mathit{x}=\mathit{F}\left(\mathit{x}\right)\mathrm{d}t+\mathrm{d}\mathit{\xi },\end{array}$$

where ξ denotes a general δ-correlated driving noise in the Itô sense. It may be state-dependent – i.e. explicitly depend on x – and contain discontinuous elements. No further specification is needed for the analysis presented here. The reconstruction of the two-dimensional drift F(x) is based on the Kramers–Moyal (KM) equation, which reads

$$\begin{array}{}\text{(4)}& \begin{array}{rl}& {\displaystyle \frac{\partial }{\partial t}}p\left(x_{\mathrm{1}},x_{\mathrm{2}},t\mathrm{|}{x}_{\mathrm{1}}^{\prime },x_{\mathrm{2}}^{\prime },t^{\prime }\right)\\ & =\sum _{i,j=\mathrm{1}}^{\mathrm{\infty }}(-\mathrm{1}{)}^{i+j}\left({\displaystyle \frac{{\partial }^{i+j}}{\partial x_{\mathrm{1}}^i\partial x_{\mathrm{2}}^j}}\right)\\ & D_{i,j}(x_{\mathrm{1}},x_{\mathrm{2}})p\left(x_{\mathrm{1}},x_{\mathrm{2}},t\mathrm{|}{x}_{\mathrm{1}}^{\prime },x_{\mathrm{2}}^{\prime },t^{\prime }\right)\end{array}\end{array}$$

in two dimensions, where $p(x_{\mathrm{1}},x_{\mathrm{2}},t\mathrm{|}{x}_{\mathrm{1}}^{\prime },x_{\mathrm{2}}^{\prime },t^{\prime })$ denotes the probability for the system to assume the state (x₁,x₂) at time t, given that it was in the state $(x_{\mathrm{1}}^{\prime },x_{\mathrm{2}}^{\prime })$ at the time t^′. The coefficients $D_{i,j}(x_{\mathrm{1}},x_{\mathrm{2}})$ of the two-dimensional Kramers–Moyal equation can be estimated – analogously to the one-dimensional coefficients as explained in Tabar (2019) – from a realisation of a two-dimensional stochastic process $\mathit{x}\left(t\right)=\left(x_{\mathrm{1}}\right(t),x_{\mathrm{2}}(t\left)\right)$. The terms D_1,0(x) and D_0,1(x) combine to the deterministic drift that governs the stochastic process:

$$\begin{array}{}\text{(5)}& \mathit{F}(x_{\mathrm{1}},x_{\mathrm{2}})={\left(D_{\mathrm{1},\mathrm{0}}\left(x_{\mathrm{1}},x_{\mathrm{2}}\right),D_{\mathrm{0},\mathrm{1}}\left(x_{\mathrm{1}},x_{\mathrm{2}}\right)\right)}^{\top }.\end{array}$$

In this work, we only consider the first-order KM coefficients. These are sufficient to uncover the deterministic non-linear features behind the stochastic data. Notice that we could formulate our method equally well in terms of the simpler Fokker–Planck equation (Risken and Frank, 1996). However, operating with the Fokker–Planck equation implicitly assumes that the stochastic process under investigation follows a Langevin equation in a strict sense; i.e. the noise term in Eq. (3) would be restricted to the case of Brownian motion. This conflicts with findings from ongoing research which indicate that the description of the driving noise ξ(t) as Brownian motion might not be valid (Rydin Gorjão et al., 2022). The use of the KM instead of the Fokker–Planck equation in this work aims at emphasising that ξ(t) might be more complex than Brownian motion and contain for example discontinuous elements. In principle, for a given stochastic process model, the higher-order KM coefficients can be used to estimate the corresponding noise parameters (see e.g. Anvari et al., 2016; Lehnertz et al., 2018; Rydin Gorjão et al., 2019; Tabar, 2019). However, this is not straightforward in two dimensions, and we deliberately refrain from an upfront selection of a process model in this work. Furthermore, a reliable estimate of higher-order coefficients in two dimensions is prevented by insufficient data density. A general derivation of the Kramers–Moyal equation can be found in Kramers (1940), Moyal (1949), Risken and Frank (1996), Gardiner (2009), and Tabar (2019).

In practice, in order to carry out the estimation of the first-order KM coefficients as defined in Eq. (4) we map each data point in the corresponding state space to a kernel density and then take a weighted average over all data points:

$$\begin{array}{}\text{(6)}& {\displaystyle }\begin{array}{rl}{D}_{\mathrm{1},\mathrm{0}}\left(\mathit{x}\right)& \sim {\displaystyle \frac{\mathrm{1}}{\mathrm{\Delta }t}}\langle \left(x_{\mathrm{1}}\right(t+\mathrm{\Delta }t)-x_{\mathrm{1}}(t\left)\right)\left|\mathit{x}\right(t)=\mathit{x}\rangle \\ & \sim {\displaystyle \frac{\mathrm{1}}{\mathrm{\Delta }t}}{\displaystyle \frac{\mathrm{1}}{N}}\sum _{i=\mathrm{1}}^{N-\mathrm{1}}K\left(\mathit{x}-{\mathit{x}}_{\mathrm{i}}\right)\left(x_{\mathrm{1},i+\mathrm{1}}-x_{\mathrm{1},i}\right)\end{array}\text{(7)}& {\displaystyle }\begin{array}{rl}{D}_{\mathrm{0},\mathrm{1}}\left(\mathit{x}\right)& \sim {\displaystyle \frac{\mathrm{1}}{\mathrm{\Delta }t}}\langle \left(x_{\mathrm{2}}\right(t+\mathrm{\Delta }t)-x_{\mathrm{2}}(t\left)\right)\left|\mathit{x}\right(t)=\mathit{x}\rangle \\ & \sim {\displaystyle \frac{\mathrm{1}}{\mathrm{\Delta }t}}{\displaystyle \frac{\mathrm{1}}{N}}\sum _{i=\mathrm{1}}^{N-\mathrm{1}}K\left(\mathit{x}-{\mathit{x}}_{\mathrm{i}}\right)\left(x_{\mathrm{2},i+\mathrm{1}}-x_{\mathrm{2},i}\right),\end{array}\end{array}$$

with $\mathit{x}=(x_{\mathrm{1}},x_{\mathrm{2}}{)}^{\top }$.

Similar to selecting the number of bins in a histogram, when employing kernel density estimation with a Nadaraya–Watson estimator for the Kramers–Moyal coefficients D_m,n(x), one needs to select both a kernel and a bandwidth (Nadaraya, 1964; Watson, 1964; Lamouroux and Lehnertz, 2009). Firstly, the choice of the kernel is the choice of a function K(x) for the estimator ${\widehat{f}}_h\left(\mathit{x}\right)$, where h is the bandwidth at a point x

$$\begin{array}{}\text{(8)}& {\widehat{f}}_h\left(\mathit{x}\right)={\displaystyle \frac{\mathrm{1}}{nh}}\sum _{i=\mathrm{1}}^{n}K\left({\displaystyle \frac{\mathit{x}-{\mathit{x}}_{\mathrm{i}}}{h}}\right)\end{array}$$

for a collection {x_i} of n random variables. The kernel K(x) is normalised as $\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}K\left(\mathit{x}\right)\mathrm{d}\mathit{x}=\mathrm{1}$ and has a bandwidth h, such that $K\left(\mathit{x}\right)=\mathrm{1}/hK(\mathit{x}/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 non-overlapping boxes”. The optimal kernel is the commonly denoted Epanechnikov kernel (Epanechnikov, 1967), also used here for the analysis of the data:

$$\begin{array}{}\text{(9)}& K\left(\mathit{x}\right)={\displaystyle \frac{\mathrm{3}}{\mathrm{4}}}(\mathrm{1}-{\mathit{x}}^{\mathrm{2}}),\mathrm{with}\mathrm{support}\left|\mathit{x}\right|<\mathrm{1}.\end{array}$$

Gaussian kernels are commonly used as well. Note that these require compact support in $(-\mathrm{\infty },\mathrm{\infty })$; thus on a computer they require some sort of truncation (even in 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 Silverman's rule of thumb (Silverman, 1998), given by

$$\begin{array}{}\text{(10)}& h_{\mathrm{S}}={\left({\displaystyle \frac{\mathrm{4}{\widehat{\mathit{\sigma }}}^{\mathrm{5}}}{\mathrm{3}n}}\right)}^{\frac{\mathrm{1}}{\mathrm{5}}},\end{array}$$

where again σ² is the variance of the time series. We note that the above formula for the ideal bandwidth has been developed for the estimation of the probability density function. As there is currently no consensus on the optimal kernel and bandwidth for the estimation of the KM coefficients, we employ an Epanechnikov kernel with bandwidth h_s throughout our work.

All numerical analyses were performed with Python's NumPy (Harris et al., 2020), SciPy (Virtanen et al., 2020), and pandas (McKinney, 2010). Kramers–Moyal analysis was performed with kramersmoyal (Rydin Gorjão and Meirinhos, 2019). Figures were generated with Matplotlib (Hunter, 2007).

4 Results

We first discuss the two drift components D_1,0 and D_0,1 (see Eq. 5) separately as functions of the two-dimensional space spanned by δ¹⁸O ratios and dust concentrations. In the component-wise analysis, the analysed component takes the role of a dynamical variable, while the respective other assumes the role of a controlling parameter. In this setting, corresponding nullclines can be computed, which reveals the bifurcation and stability structure of the two individual drift components. Intersections of the two components' nullclines yield fixed points of the coupled system, which are stable if both nullclines are stable at the intersection.

Figure 3 Two-dimensional drift reconstruction. (a) PDF of the two-dimensional record, with projections onto both dimensions. Blue and orange dots represent the individual data points from Greenland stadials (GSs) and Greenland interstadials (GIs), respectively. Contour lines are obtained from a kernel density estimate of the data distribution. The dotted contour line indicates a chosen cutoff data density of >0.015 data points per pixel; regions in the state space with lower data density are not considered in the analysis. One pixel has the size of 0.015 × 0.015 in normalised units. (b) The reconstructed vector field F according to Eq. (5). Regions of convergence are apparent and correspond to the GI and GS states of the record. (c) The dust component D_0,1 of the reconstructed drift. The dust's nullcline exhibits an S shape, with two stable branches (orange) and an unstable one in between (red), indicative of a double-fold bifurcation with δ¹⁸O as a control parameter. (d) The δ¹⁸O component D_1,0 of the reconstructed drift. Here, the nullcline is comprised of a single stable branch (orange). The position of the δ¹⁸O fixed point varies with the value of the dust. Fixed points of the coupled system are given by the intersections of the two components' nullclines, marked with an X in panel (b).

4.1 Double-fold bifurcation of the dust

The estimated dust drift $D_{\mathrm{0},\mathrm{1}}(x_{\mathrm{1}}={\mathit{\delta }}^{\mathrm{18}}\mathrm{O},x_{\mathrm{2}}=\mathrm{dust})$ is displayed in Fig. 3c. This coefficient dictates the deterministic motion of the system along the dust direction; therein the δ¹⁸O ratio takes the role of the controlling parameter. We can trace the nullcline's branches, which take a general “S” shape as we vary δ¹⁸O. Hence, depending on the value of δ¹⁸O, there is either one or three fixed points for the motion along the dust direction: for approximately δ¹⁸O $<-\mathrm{1.0}$, there is one stable fixed point; for approximately $-\mathrm{1.0}<{\mathit{\delta }}^{\mathrm{18}}$O <0.9, there are three fixed points, two stable ones and an unstable one between them; for approximately δ¹⁸O >0.9, there is again just one stable fixed point. In fact, the merger of the nullcline's lower stable branch and unstable branch is not fully captured by the reconstruction due to too-low data density in the corresponding region (see Fig. 3c). With the position of these stable fixed points depending continuously on δ¹⁸O ratios, we find here the characteristic form of a double-fold bifurcation, in which δ¹⁸O takes the role of a control parameter.

The dust nullclines' structure supports the possibility for abrupt transitions in two ways: either random fluctuations move the system across the unstable branch (if present, depending on the value of the control parameter), or the control parameter, in this case δ¹⁸O, crosses a bifurcation point, and the currently attracting stable branch merges with the unstable branch. In both cases, the system will transition fairly abruptly to the alternative stable branch. A rate-induced transition seems implausible in this case since the unstable branch is approximately constant with respect to a change in the control parameter (i.e. δ¹⁸O). Thus, a crossing of the unstable branch by means of a rapid shift in δ¹⁸O seems highly unlikely.

4.2 Coupling of the δ¹⁸O drift with the dust

We now focus on the reconstructed drift $D_{\mathrm{1},\mathrm{0}}({\mathit{\delta }}^{\mathrm{18}}\mathrm{O},\mathrm{dust})$ of the δ¹⁸O ratios (Fig. 3d). Its nullcline appears to be an explicit function of the dust; i.e. for each value of dust concentration there is a single stable fixed point along the δ¹⁸O dimension. The position of the fixed point changes with the value for dust in a continuous manner, with a high rate of change for intermediate dust values and a small change for more extreme dust values. These findings suggest that δ¹⁸O follows a mono-stable process whose fixed point is subject to change in response to an “external control” imposed by the dust.

4.3 Combined two-dimensional drift

Figure 3b shows the two-dimensional drift field F(δ¹⁸O,dust) of the coupled system given by Eq. (5). The two fixed points which arise from the intersections of the dust nullcline's stable branches with the δ¹⁸O stable nullcline fall well within the regions of the state space associated with Greenland stadials and interstadials, respectively. The stable regime (δ¹⁸O $\sim -\mathrm{1}$, dust $\sim -\mathrm{1}$) can be identified with Greenland stadials, while the stable regime (δ¹⁸O ∼0.5, dust ∼1) corresponds to Greenland interstadials. Similarly, we can locate an unstable fixed point roughly in between the two observed stabled fixed points of the coupled system. Judging from Fig. 3b the unstable fixed point resembles a saddle with a convergent drift along the δ¹⁸O direction and divergence along the connection line between the two stable fixed points. The system's bistability is inherited from the dust's drift and is not enshrined in the δ¹⁸O ratios. As mentioned previously, Fig. 3b suggests that – starting from a stable fixed point – perturbations along the δ¹⁸O direction will not entail state transitions but instead simply decay until the system reaches the δ¹⁸O nullcline again. In contrast, perturbations along the dust direction may shift the system into the other respective basin of attraction.

Figure 4 Redrawing of Fig. 3 in a rotated state space. The variables p₁ and p₂ represent the rotated time series, onto which the same KM analysis is performed as before. We can observe that even in this rotated setting we cannot disregard the coupling of the two variables. The doubled-fold structure is occluded by the rotation (see drift of the variable p₁ in panel c). The drift of p₂ remains dependent on p₁ (see panel d). We can thus conclude that the observed coupling is not an artefact of the initial state space used and is an intrinsic characteristic of the two proxies.

4.4 Rotation of the state space and the presence of a non-negligible interplay of the dust and δ¹⁸O

Above we argue for the existence of a double-fold bifurcation in the dust variable. In order to show that the coupling of the dust and δ¹⁸O is not a spurious result of the initial state space, we conduct an analogue analysis using a rotated state space. To rotate the state space we employ principal component analysis and obtain a new set of variables $\mathit{p}=(p_{\mathrm{1}},p_{\mathrm{2}})$, with $\mathit{p}=\mathbf{U}({\mathit{\delta }}^{\mathrm{18}}\mathrm{O},\mathrm{dust}{)}^{\top }$, where U is given by

$$\begin{array}{}\text{(11)}& \mathbf{U}=\left[\begin{array}{cc}-\mathrm{0.707}& -\mathrm{0.707}\\ -\mathrm{0.707}& \mathrm{0.707}\end{array}\right].\end{array}$$

In Fig. 4 we redraw Fig. 3 in the rotated state space; we observe that (i) the nullcline of p₁ is now almost independent of p₂, and (ii) the p₂ nullcline is still strongly dependent on p₁, while none of the rotated variables show any bifurcation. Overall, the dynamics of the dust–δ¹⁸O can be explained as we introduce in Sect. 4.2, with two basins of attraction being separated by a saddle. In particular, the assessment of the drift in the rotated state space shows that the data cannot be described by a simple two-dimensional double-well potential with two axes of symmetry and decoupled dynamics along them.

5 Discussion

We use the two-dimensional Kramers–Moyal equation to investigate the deterministic drift of the combined dust and δ¹⁸O record from the NGRIP ice core for the time interval 59–27 kyr b2k, which exhibits pronounced DO variability. The reconstructed stability structure with two basins of attraction and a separating saddle is consistent with the regime switches observed simultaneously in both components of the record: in the δ¹⁸O–dust plane, the basins of attraction are located such that a transition from one to the other entails a change in both components. However, the analysis of the vector field (Fig. 3b) does not indicate any clear paths the system takes in order to transition between stadial and interstadial states. The shape of the nullclines can, in principle, allow for a situation where a perturbation along the δ¹⁸O direction pushes the dust across its bifurcation point, triggering a transition of the dust, which in turn stabilises the δ¹⁸O perturbation. The combined drift F(x₁,x₂), however, exhibits strong restoring forces along the δ¹⁸O direction, which render this mechanism rather implausible. Viewed from either stable fixed point, perturbations along the dust direction could in contrast push the system across the basin boundary relatively easily. Certainly, a combination of noise along both directions may also be able to drive the system across the region of weak divergence that separates the two attractors. We note that the mild relaxation that is typical for Greenland interstadials cannot be explained by the results of this analysis alone. As mentioned previously, an estimation of the KM coefficients for the individual, univariate δ¹⁸O and dust time series indicates that at least the δ¹⁸O noise comprises non-Gaussian and potentially discontinuous components, which could play a central role with respect to the transition between the two identified stable states of the drift (Rydin Gorjão et al., 2022). However, there remain discrepancies to be reconciled in the analysis of the higher-order KM coefficients of the individual δ¹⁸O and dust time series such that arguments about the role of non-Gaussian noise in the regime transitions remain speculative at this point. Ideally, higher-order KM coefficients should be computed for the two-dimensional record; however, this is prevented by the low data resolution.

In the following, we discuss how the results presented here relate to the findings of previous studies. An important branch of research around DO events draws on low-dimensional conceptional modelling and, related to that, inverse modelling approaches with model equations being fitted to ice core data. Many of these studies build on stochastic differential equations and in particular on Langevin-type equations. Our study follows the same key paradigm, regarding the paleoclimate record as the realisation of a stochastic process. However, as far as we know, it is the first study to assess the two-dimensional drift non-parametrically in the δ¹⁸O–dust plane.

For the period investigated here Livina et al. (2010) attested bistability to both the δ¹⁸O and the dust component individually by fitting a Langevin process to a 20-year-mean version of the NGRIP record. Later Kwasniok (2013) and Lohmann and Ditlevsen (2018a) showed – using techniques from Bayesian model inference – that a two-dimensional relaxation oscillator model outperforms a simple double-well potential in terms of simulating the NGRIP δ¹⁸O record. Such a relaxation oscillation still relies on a fundamental bistability in the variable that is identified with δ¹⁸O ratios. A physical interpretation for a FitzHugh–Nagumo-type DO model is provided by Vettoretti et al. (2022).

Our results contradict the interpretation that δ¹⁸O ratios and therewith Greenland temperatures bear an intrinsic bistability. In the two-dimensional setting, the apparent two-regime nature of the δ¹⁸O record can be explained by the control that the dust exerts on the δ¹⁸O fixed points and the corresponding location of the two stable fixed points in the two-dimensional drift. Since we find the bistability of the reconstructed coupled system rooted in the dust, our analysis suggests that the atmosphere may have played a more active role in stabilising the two regimes that dominated the last glacial's Northern Hemisphere climate than many AMOC-based explanations of the DO variability suggest (Ganopolski and Rahmstorf, 2002; Clark et al., 2002; Vettoretti and Peltier, 2018; Li and Born, 2019; Menviel et al., 2020). Similarly, the observation that dust perturbations may induce state transitions may be seen as a hint that random perturbation of the atmospheric circulation can trigger DO events as proposed by, for example, Kleppin et al. (2015). In this regard, it should be noted that a multistability of the latitudinal jet stream position has been suggested – although in a somewhat different setting and sense – based on an investigation of reanalysis data of modern climate (Woollings et al., 2010). In contrast, one would not expect two distinct stable Greenland temperature regimes with all controlling factors kept fixed. This is in line with the bistability of the dust drift and the mono-stable δ¹⁸O drift revealed in our analysis.

Clearly, the state space spanned by δ¹⁸O and dust is a very particular one. On the one hand, the interpretation of the two proxies as indicators of Greenland temperatures and the hemispheric circulation state of the atmosphere bears qualitative uncertainties and should certainly not be considered a one-to-one mapping. On the other hand, other climate subsystems not directly represented in the data analysed here, like the AMOC for example, are likely to have played an important role in the physics of DO variability as well. Even if δ¹⁸O ratios and dust concentrations were to exclusively represent Greenland temperatures and the atmospheric circulation state, the recorded climate variables were certainly highly entangled with other climate variables such as the AMOC strength, the Nordic Sea's and North Atlantic's sea ice cover, or potentially North American ice sheet height (e.g. Menviel et al., 2020; Li and Born, 2019; Boers et al., 2018; Zhang et al., 2014; Dokken et al., 2013). In our analysis, such couplings are subsumed in the δ-correlated noise term ξ – an approach which may rightfully be criticised to be overly simplistic. However, given the lack of climate proxy records that jointly represent more DO-relevant components of the climate system on the same chronology, the chosen method reasonably complements existing data-driven investigations of DO variability. For example Boers et al. (2017) similarly examined the dynamical features of the combined δ¹⁸O–dust record. They proposed a third-order polynomial two-dimensional drift in combination with a non-Markovian term and Gaussian white noise to model the coupled dynamics. While our approach is limited to a Markovian setting, it allows for more general forms of drift (and noise). Being non-parametric, it does not rely on prior model assumptions in this regard. It is not per se clear how the couplings to “hidden” climate variables (i.e. those not represented by the analysed proxy record) influence the presented drift reconstruction, and there is certainly a risk of missing a relevant part of the dynamics.

6 Conclusions

We have analysed the records of δ¹⁸O ratios and dust concentrations from the NGRIP ice core from a data-driven perspective. The central point of our study was to examine the stability configuration of the coupled δ¹⁸O–dust process by reconstructing its two-dimensional drift. Our findings indicate a mono-stable δ¹⁸O drift whose fixed point's position is an explicit function of the dust. The dust variable, in contrast, seems to undergo a double-fold bifurcation parameterised by δ¹⁸O, with a change from a single (stable) fixed point to three fixed points (two stable, one unstable) and again to a single (stable) fixed point, from small to large values of the δ¹⁸O ratio. Together, the drift components yield two stable fixed points in the coupled system surrounded by convergent regions in the δ¹⁸O–dust state space, in agreement with the two-regime nature of the coupled record. Judging from the reconstructed drift, perturbations along the dust dimension are more likely to trigger a state transition, which points to an active role of atmospheric circulation in DO variability.

Importantly, our findings question the prevailing interpretation of the two regimes observed in the isolated δ¹⁸O record as the direct signature of an intrinsic bistability. Such an intrinsic bistability can be confirmed only for the dust variable. Regarding δ¹⁸O ratios as a direct measure of the local temperature, it seems plausible that not the temperature itself is bistable but rather that the bistability is enshrined in another climate variable – or at least a regional-scale climate process or a combination of processes – that drives Greenland temperatures. The apparent two-regime nature of the δ¹⁸O record would thus only be inherited from the actual bistability of other processes. This may be the atmospheric circulation as represented by the dust proxy or another external driver not directly represented by the analysed data.

Similar investigations to ours should be applied to other pairs of Greenland proxies to investigate the corresponding two-dimensional drift. Finally, our study underlines the need for higher-resolution data, as the scarcity of data points is a limiting factor for the quality of non-parametric estimates of the KM coefficients.