2.1 Rolling window metrics

For the rolling window metrics we follow the methodology first outlined in a climate context by Held and Kleinen (2004) and subsequently used by other studies, including Dakos et al. (2008, 2012a), Lenton et al. (2012a, b), Livina and Lenton (2007), and the Early Warning Signals Toolbox developed based on this work (documented at http://www.early-warning-signals.org (last access: 15 October 2018) and available in as the earlywarnings package in R; R Foundation for Statistical Computing, 2016). After selecting the dataset and for the pre-PETM analysis terminating it just prior to the hypothesised transition to avoid biasing the analysis, the data are first interpolated (using linear interpolation by default with the interp1 function in MATLAB; The MathWorks Inc., 2016) to provide the equidistant data points required for rigorous statistical analysis and assumed by the AR(1) model. However, interpolation itself can introduce statistical artefacts into the analysis as, by definition, the addition of interpolated data points increases self-similarity and thus autocorrelation in the dataset. In palaeorecords this tends to result in an artificial increase in autocorrelation in parts of the dataset with either sparser data points or complete gaps in the data, but in this record there is no marked systematic shift in data time steps (Fig. S5). As a result we also analyse non-interpolated data in order to assess the sensitivity of our results to interpolation. Following this, the data are then detrended by subtracting the smoothed dataset, estimated with a Gaussian kernel smoothing function (using the ksmooth function in R), in order to remove any long-term trends as these are not the focus of the analysis. This makes the dataset stationary – a necessary prerequisite for time series analysis – but this also somewhat reduces the value of lag 1 autocorrelation in the results. Bandwidth is an important consideration in this process and is adjusted heuristically for the datasets in order to best remove long-term trends but leaving short-term fluctuations, in this case giving a Gaussian kernel bandwidth of 0.1. This removes long-term secular trends and orbital cyclicity (> 100 kyr) (Figs. 2, 3a and b), leaving only the short-term noise that reveals the resilience of the underlying longer-term (> 10 kyr) processes.

An autoregressive model of order 1 (AR(1)) is fitted to the data within a rolling window (using the generic_ews function of the earlywarnings package in R). The AR(1) model is of the form $x_{t+\mathrm{1}}={\mathit{\alpha }}_{\mathrm{1}}{x}_t+{\mathit{\epsilon }}_t$, fitted by an ordinary least-squares method with a Gaussian random error and a constant time step. Following previous studies the default window size is set at half the length of the dataset, but as part of our sensitivity testing we also repeat our analyses for window sizes between 25 % and 75 % (Figs. S1–S4). The choice of window length is a trade-off between dataset resolution and the reliability of the estimate of the indicator, with a short window allowing shorter-term changes in indicators to be tracked at the cost of lower estimate reliability and vice versa. On the same rolling window the absolute skewness, kurtosis, and standard deviation of the dataset are also calculated (also using the generic_ews function of the earlywarnings toolbox in R). Detrended fluctuation analysis h value (DFA-h) was also used as an alternative measure to AR(1) for short-term memory and critical slowing down in the dataset bins (performed using the DFA function of the fractal package in R). DFA extracts the fluctuation function over a window s, and if the data are long-term power-law correlated, the fluctuation function F(s) increases as a power law: F(s)∝s^h, where h is the DFA fluctuation exponent (Peng et al., 1994) and reaches value 1.5 at a critical transition (Lenton et al., 2012b; Livina and Lenton, 2007).

Finally, the likelihood of there being a real trend in the results is calculated by estimating the non-parametric Kendall rank-correlation statistic (τ), which measures the strength of an indicator's tendency to increase (≫ 0) or decrease (≪0) against the null hypothesis of randomness (∼0) (also using the generic_ews function of the earlywarnings toolbox in R). However, this statistic is most robust when the trend is consistent over a long period, while increasing but oscillating trends or trends only at the very end of the record can produce weak or even negative values despite a clearly visible trend (Dakos et al., 2012a). We calculate a p value for each metric by bootstrapping our detrended datasets to generate 1000 surrogate records (or ARIMA model-generated for AR(1) and DFA-h) with equivalent mean and variance, re-calculating the metric and Kendall τ value for each, and finding the proportion of Kendall τ values equal to or greater than that of the original palaeorecord (Dakos et al., 2008).

2.2 Binned metrics

As well as performing rolling window time series analysis, we also measure AR(1), SD, skewness, and kurtosis on data (detrended but not interpolated) binned into pre- and post-event bins and excluding the events themselves to provide simple before–after comparisons of changes across the events. To this end the datasets (n=2302) were binned into pre-PETM (n=1331), post-PETM (n=921), PETM to ETM2 (n=593), and post-ETM2 (n=240) bins, excluding data points from during each event so as to avoid biasing by extreme or missing data. A p value is calculated for each metric using a permutation test (i.e. by reshuffling and repartitioning the before–after event data into the same-sized bins 1000 times and comparing the metrics' resultant before–after differences with the observed metric before–after difference), except for AR(1) and DFA-h for which we instead use AR(1) model-derived surrogate data to compare against (i.e. by generating 1000 surrogate datasets with the same AR(1) value, mean, and variance as the before-event bin over the length of the after-event bin and comparing this distribution to the observed after-event AR(1) value).

Table 1Values of binned metrics for both benthic δ¹⁸O (a) and δ¹³C (b). Bins are of all data (detrended but not interpolated) before the PETM, all data after the PETM (including data from within ETM2), all data between the PETM and ETM2, and all data after ETM2. The indicators are AR(1) coefficient, detrended fluctuation analysis h value (DFA-h), standard deviation (SD), skewness (SKEW), and kurtosis (KURT). Each value is followed by a p value (in italics within parentheses and preceded by a ^* if significant) computed using a permutation test, except for AR(1) and DFA-h for which we instead use AR(1) model-derived surrogate data to compare against (see Methods for details).

Download Print Version | Download XLSX

2.3 Non-parametric drift–diffusion–jump model

A model-based alternative to the time series analysis methods (whether rolling window or bin based) above is to fit a general non-parametric drift–diffusion–jump model to the dataset with as a surrogate for an unknown data-generating process (Carpenter and Brock, 2011; Cox and Ross, 1976; Dakos et al., 2012a; Johannes, 2004). In this model functions are estimated for drift, diffusion, and jump processes using non-parametric regression; drift measures the local rate of change, diffusion measures the standard deviation of the relatively small shocks that occur at each time step, and jumps are large intermittent shocks. The conditional variance of the data is also estimated from the non-parametric regression and represents the variance of the data from their conditional expectation estimated using kernel regression. We use the ddjnonparam_ews function in the earlywarnings package on R, with the default options of a bandwidth of 0.6 and 500 points for computing the kernel. We use raw data for this analysis, with no detrending or interpolation and without log transforming the data first. In interpreting the results we focus on the general long-term trends in the estimated terms as many of the shorter-term fluctuations potentially represent model overfitting.

2.4 Limitations

Despite positive results in some palaeoclimate EWS studies (Lenton, 2011), there are several potential issues with searching for resilience indicators in palaeorecords. Palaeorecords suffer from greater dating uncertainties and a less frequent sampling rate than is possible with modern climate data, making robust time series analysis more challenging. Most indicators also do not reveal exact information about the nature of the transition itself, with increasing slowing down and variability detected prior to both catastrophic and non-catastrophic transitions featuring a bifurcation in phase space and even before non-catastrophic transitions without a bifurcation (Kéfi et al., 2013). Concerns have also been raised over the likelihood of producing false positives (when EWSs appear to indicate an impending transition which never occurs) or false negatives (i.e. a “missed alarm” when EWSs may be entirely absent prior to a known critical transition) and the extent to which these methods are dependent on subjective parameter choices (Boettiger et al., 2013; Lenton, 2011). There is a risk that selecting and analysing known or suspected critical transitions in the palaeoclimate record is particularly liable to false positives, as positive indicators at the transitions could potentially have occurred purely by chance rather than due to systemic instability (Boettiger and Hastings, 2012a). However, it has been argued that EWSs can be reliably detected if both increasing autocorrelation and variance are seen prior to the transition rather than one of these indicators alone (Ditlevsen and Johnsen, 2010). Detecting multiple, consistent, and robust signals from the indicators can be indicative of decreasing system resilience even if a catastrophic transition is not reached or is instead triggered by an external perturbation rather than internal processes (Dakos et al., 2015).

3.1 Rolling window metrics

Rolling window metrics prior to the PETM reveal a sudden increase in AR(1) and SD after ∼58.2 Ma in the interpolated benthic δ¹⁸O record associated with a step in the data (and the benthic δ¹³C record peaking), which despite a temporary drop in DFA-h suggests some degree of destabilisation of the slow climate system prior to the PETM (Figs. 2 and S1). However, the subsequent decline of standard deviation after ∼57 Ma (likely to partially be the result of earlier extreme data points from ∼58.5–59.5 Ma leaving the rolling window) does not support a tipping point involving deep ocean temperature at the PETM, which is also indicated by the non-significant bootstrapped p values for the metrics of both interpolated and non-interpolated data. The non-interpolated DFA-h metric does show a significant increase (p=0.015), which would suggest systemic slowing down, but this does not match the non-interpolated AR(1) metric. Alternatively, these results could represent a “missed alarm” as the shorter-term climate dynamics that might be critical to the dynamics of the tipping point are not sufficiently resolved by the available data.

The benthic δ¹³C record shows clearer evidence of declining resilience in the slow components of the carbon cycle, with long-term increases in AR(1), DFA-h, and SD in the run-up to the PETM with steps at ∼58.2 and ∼57.3 Ma which are consistent across the sensitivity analyses (absolute skewness also increases, while kurtosis declines up to the PETM; Fig. S1). Bootstrapped p values indicate that the δ¹³C SD trend is significant (p=0.002) for the interpolated data, while the δ¹³C AR(1) and DFA-h trends are significant (p=0.024 and p=0.021, respectively) for the non-interpolated data. This supports a long-term slowing down in benthic δ¹³C in the late Palaeocene, which may reflect a gradually forced destabilisation of the global carbon cycle prior to the PETM.

Rolling window analysis across the whole of the late Palaeocene–early Eocene (LPEE) interval suggest but cannot prove systemic changes in carbon cycle and climate (in)stability across both the PETM and ETM2 (Figs. 3 and S2). Between the PETM and ETM2 δ¹⁸O AR(1) and DFA-h increase up until ∼200 kyr before ETM2 and SD experiences a small temporary increase followed by a larger decrease. In contrast, all metrics for δ¹³C experience a rapid jump during the PETM and then remain relatively stable until ETM2. Following ETM2, δ¹⁸O AR(1) and DFA-h increase significantly while SD increases slightly, whereas for δ¹³C all metrics (as well as absolute skewness and kurtosis; Fig. S2) consistently increase. However, the bootstrapped p values indicate that none of these trends are significant for the interpolated data, but that for the non-interpolated data the increases in AR(1) and DFA-h for both δ¹⁸O and δ¹³C are highly significant (p=0 to 0.001). This indicates that there is some evidence for slowing down – but not for increased variability – in both the δ¹⁸O and δ¹³C data and therefore the long-term climate system and carbon cycle across the LPEE interval into the Eocene, but this is dependent on not interpolating the data prior to the analysis. It should also be recognised that the abrupt shifts in δ¹³C at the PETM and ETM2 are not fully removed by detrending prior to the analysis; hence, they are at least partly responsible for the upward steps in the indicators at the events.

Figure 4Non-parametric drift–diffusion–jump model functions of benthic δ¹⁸O (a), (c), (e), (g) and δ¹³C (b), (d), (f), (h) across the PETM and ETM2. Panels (a) and (b) illustrate the palaeorecords (black crosses), with the panels below illustrating the model functions for (c, d) conditional variance (blue line), (e, f) diffusion (green line), (g, h) and jump intensity (red line) for each palaeorecord (drift is not plotted). The PETM and ETM2 are marked by the grey bars.

Download

3.2 Binned metrics

To address the issue of large excursions failing to be removed by detrending for the rolling window metrics, as well as the issue of data gaps caused by dissolution at the peak of each event (Littler et al., 2014), we calculate aggregate metrics (i.e. no rolling windows) on the binned data (excluding data from within the events) (Table 1). The binned metrics show significant increases in AR(1) across both the PETM and ETM2 for δ¹⁸O and δ¹³C, with DFA-h also significantly increasing across the PETM and ETM2 for δ¹⁸O and across ETM2 for δ¹³C. SD significantly falls in the interval between the PETM and ETM2 and then significantly increases across ETM2 for both datasets, but the significant increase for δ¹³C for the all data after the PETM is likely to be biased by the inclusion of extreme data points within ETM2. These results are consistent with the rolling window metrics, except for δ¹³C SD, which in the rolling window metrics is higher between PETM and ETM2 than before the PETM and is likely to be the result of excluding the extreme data values during the PETM itself in the binned analysis (with SD beginning to drop just before the ETM2 in the rolling window analysis once the PETM leaves the window). These results indicate both the long-term climate system and carbon cycle slowed down to some extent after both the PETM and ETM2 (but became less variable following the PETM itself until ETM2), providing support for the slow parts of global carbon–climate system being progressively destabilised through the LPEE interval and into the Eocene by both of the hyperthermal events but not for any tipping points.

3.3 Non-parametric drift–diffusion–jump model

Fitting a non-parametric drift–diffusion–jump model to the datasets provides independent model-based metrics to compare to the rolling window metrics, with terms for the conditional variance measuring variance from dataset's conditional mean (estimated by kernel regression), diffusion measuring the standard deviation of regular small shocks at every time step, and jump intensity measuring either irregular large shocks or flickering (Carpenter and Brock, 2011; Dakos et al., 2012a) (Fig. 4). For benthic δ¹⁸O this model reveals an overall increase in conditional variance and diffusion and a reduction in jump intensity ∼2 Myr before the PETM, followed by a slight decrease in conditional variance after the PETM and intermittent spikes in jump intensity and conditional variance during and following ETM2. This suggests the climate shifted to a state with higher variability featuring regular small shocks ∼2 Myr prior to the PETM and became slightly less variable following the PETM, but featured larger irregular shocks during and after ETM2. While this suggests some degree of climate instability in the ∼2 Myr before the PETM and following ETM2, there is no evidence of a critical transition in the climate system at the PETM or ETM2 themselves. In contrast, the benthic δ¹³C model reveals decreasing diffusion and increasing conditional variance and jump intensity in the 1.5 Myr run-up to the PETM, indicating increasing total variability driven by large irregular shocks and consistent with a critical transition being approached in the carbon cycle at the PETM (Dakos et al., 2012a). Conditional variance and jump intensity remain high and diffusion remains low for ∼1 Myr after the PETM, before reversing ∼300 kyr before ETM2 except for brief spikes in both diffusion and jump intensity during and after ETM2. This indicates that variability in the carbon cycle remained high and driven by large shocks for ∼1 Myr after the PETM, but that variability mostly shifted towards smaller regular shocks prior to and after ETM2. This difference in the sources of variance prior to each event (the PETM is preceded by elevated jump intensity and overall conditional variance, whereas ETM2 is preceded by increased diffusion and decreased jump intensity) suggests potentially differing carbon cycle dynamics and drivers prior to each event. The shift in variability before ETM2 also slightly precedes the biotic turnover detected in both marine and terrestrial records in the ∼200 kyr prior to ETM2, despite there being no obvious shift in the palaeorecords that may have driven this turnover (Westerhold et al., 2018). Overall the δ¹³C DDJ results are consistent with elevated carbon cycle instability following the PETM, but suggest that ETM2 was not preceded by the same dynamics as the PETM.