Deciphering the Palaeocene–Eocene thermal maximum by Granger causality test

The Palaeocene–Eocene thermal maximum is a global warming period (~ 56 Ma), which is marked by a sharp negative carbon isotope excursion (CIE) that caused by the injection of massive isotopically-light carbon into the ocean-atmosphere. It is often considered that the carbon injection caused global warming. However, several studies have suggested that warming and environmental perturbations precede the onset of the CIE. Here we present Granger test to investigate the detailed mechanisms of this event. We show a shift from climate-warming driving carbon-emission scenario to a scheme in which carbon-injection causing global-warming during the CIE. The initial carbon emission might be from methane hydrates dissociation and/or permafrost thawing, possibly linked with astronomical paced warming. This change of causal direction may result from the warming feedback of the emitted carbon and additional carbon from other sources, such as volcanism, bolide impact, oxidation of marine organic matter, and wildres burning peatlands.


Introduction
The Paleocene-Eocene thermal maximum (PETM, ~ 56 Ma) is a period of rapid and prolonged global warming 1 . Paleothermometer proxies (e.g. δ 18 O (oxygen isotope) and TEX 86 ('tetraether index' of tetraether lipids consisting of 86 carbon atoms)) analyses suggest that global temperature rise approximately 5-8°C during this interval 2,3 . Associated with the climate warming is a dramatic negative carbon isotopic excursion (CIE) of 3-5‰ 4 , which has been interpreted to be caused by the injection of massive 13 C-depleted carbon into the ocean and atmosphere 5 . It is generally considered that the massive carbon injection into the ocean-atmosphere system caused global warming and environmental perturbations 6,7 . However, several studies have found that global warming and environment change precedes the carbon isotope excursion, and proposed that the initial warming may have triggered the release of the isotopically-light carbon sources 8 − 14 . Thus, despite strong correlation exists between the carbon cycle perturbation and the global warming and environment change, the exact driving mechanism between these factors still remains unclear.
In this study, we use the Granger causality test to further investigate the mechanisms between carbon release (inferred from δ 13 C data) and global warming (reconstructed from δ 18 O data). The Granger causality test is rst introduced in econometrics by Nobel Prize winner Clive W. J. Granger 15 . It is used to examine the causality relationship between variables of time series. In brief, if variable Y can be explained by lagged variable X in the statistic model, it is then considered that X granger-cause Y. This method has been used to investigate the causality between greenhouse gas concentrations and global temperature, and has yielded different conclusions [16][17][18][19][20] . Furthermore, in addition to data from 1850 to present day, his method has been applied to climatological series from the more recent past (e.g. time series data recovered from ice cores, covering ~ 800 kyr) [21][22][23][24][25] . In this study, we explore the potential application of this method to deep-time climatological data recovered from core sediments and discuss its limitations. Loading [MathJax]/jax/output/CommonHTML/jax.js Geochemical data generated of samples collected from the lower part of the drill core to the upper part would form a time series dataset. In this study, previously published data from the Millville core (ODP 174AX) are analysed, which is drilled in the expanded New Jersey shelf in the North Atlantic Ocean.
Coupled δ 13 C and δ 18 O data are available across the Paleocene-Eocene interval of the drill core at 2 mm gap ( Fig. 1) 26 . Based on sedimentation rates estimated from the nearby area, the time gap is about 20-25 years 12 . However, carbon cycle-climate modelling has estimated the sampling resolution of ~ 40 years 14 . The exact duration between two adjacent data points is hard to be precisely quanti ed. However, this would not affect the test result, as the sequence of the data is not changed. Oxygen isotope (δ 18 O) varies as a result of temperature uctuations and thus paleoclimate can be reconstructed by analysing δ 18 O of the sediments 27 . Carbon isotope record is used to re ect the addition of δ 13 C-depleted carbon into the ocean and atmosphere, although they do not follow a linear relationship 28 . The data are divided into four intervals, with the rst interval (Q1) covers the data before the onset of the CIE, the second and third intervals (Q2 and Q3) cover the onset of the CIE 26 and the fourth interval (Q4) covers the data after the CIE. The boundary between Q2 and Q3 is drawn based on the carbon isotope pro le, where a pulse of large negative CIE occurred within the onset of the PETM CIE (Fig. 1).

Granger causality test
Granger causality was rst introduced in econometrics to test potential causality relationships between variables of time series data 15 . It is based on incremental forecasting power. Here we consider two nested models, the unrestricted regression model where µ 1 and µ 2 are constants; a, b, and c are the coe cients of the models; ε t and v t are univariate white noise; i is the lag of the models and k is the highest lag. For example, Y t is the current sample, and Y t-1 is the sample before Y t . Y t-1 is one sampling gap stratigraphically lower than sample Y t . The Granger causality relationship between variable Y and variable X is determined by comparing the estimation accuracies for this climatic factor by the unrestricted and the restricted models. The null hypothesis of noncausality corresponds to: A signi cant statistics (P-value < 0.05) implies that the null hypothesis is rejected. The increase in the explanatory power of the unrestricted model (including variable X) compared to the restricted model (not including variable X) suggests causality from X to Y. Before we perform a Granger causality test, it is necessary to determine if the time series is stationary, as nonstationary of time series may cause spurious causality results 29 . This can usually be done with the unit root test process 30 . If the series has a unit root, then the time series is not stationary. Nonstationary time series data can achieve stationary by taking the rst difference. For nonstationary time series data, the Toda-Yamamoto method can be applied 31 . This procedure would require the maximum order of integration (d) of the investigated series.
The optimal lag number is selected by using Akaike information criteria (AIC) 32 .

Unit root test
In this study, the augmented Dicky-Fuller (ADF) test method is used for unit root test 30 . The null hypothesis is that the series has a unit root and contains a stochastic trend. The model of the augmented Dickey-Fuller (ADF) test is speci ed as follow: where y is the investigated variable and w t is a random error term. The lagged rst differences of the dependent variables provide correction for possible serial correlation. The optimum number of lags (p), in this paper, was selected based on log likelihood 33 , Hannan-Quinn criterion 34 , and Akaike information Criterion 32 . The null hypothesis is given by γ = 0. The alternative hypothesis is that the series is stationary. Unit root tests are also performed for cases 1) without trend with constant and 2) without trend or constant and the decision is made based on three-step ADF test 35 . Unit root test was performed on the data for the whole period and the four quarters individually with results listed in Table 1. Notes: The numbers in the parentheses indicate the optimal lag length suggested by the Akaike Information Criterion (AIC). t + c means the estimation with trend and constant. c means estimation with constants only.
In order to provide some robustness in deciding the stationarity of the series, additional tests like

Phillips-Perron (PP), Kwiatkowski-Phillips-Schmidt-Shin (KPSS) and Dickey-Fuller Generalized Least
Squares (DF-GLS) are conducted [36][37][38][39] . The PP and DF-GLS unit root tests have a null hypothesis of a unit root process, while the KPSS test have a null hypothesis of a stationary series. These results are listed in Table 2 Table 3. quarters are stationary, despite the whole period δ 13 C data are not stationary (Table 2). Likewise, stationarity has also been suggested for recent glaciation and interglaciation transitions that are tend to be assumed nonstationary 22,40 . The difference between stationary and nonstationary conclusions drawn by different authors has been attributed to the differing time-spans of data analysed 22 .

Causality test results
The rst quarter (Q1) covers the majority of the data before the CIE (Fig. 1). In theory, the δ 13 C and δ 18 O should behave independently in steady-state. However, a causality relationship is detected from δ 13 C to δ 18 O (Fig. 1, Table 3). Generally, the δ 13 C driving δ 18  The second quarter (Q2) covers the initial excursions for both the δ 13 C and δ 18 O data (Fig. 1). Previously, the theory that global warming driving the carbon release mainly comes from observed geological records that the proxies for environmental perturbations start to shift before the CIE 9-13 . It usually needs expanded sections with high sedimentary rate to clearly discern the sequence of different datasets. For example, at the Wilson Lake section (New Jersey shelf), climate warming precedes the onset of the CIE as inferred from the TEX 86 thermometer record 12 . This is also supported by the biotic change that the onset of the Apectodinium (a subtropical genus) acme is detected ~ 40 cm below the CIE 12 . The nearby Bass River section also shows a similar pattern, with half of the total warming occurred below the CIE 12 . The Millville site is close to the Wilson Lake and Bass River sections, however, δ 13 C and δ 18 O data from this section start to shift simultaneously (Fig. 1). It is not feasible to spot their sequence and causality relationship based on the isotope stratigraphy pro le alone. Even so, our Granger causality test result is able to show a signi cant (at 5%) Granger causality from the δ 18 O to δ 13 C in this period ( Fig. 1; Table 3 Climate warming causing the carbon injection that responsible for the CIE is consistent with the hypothesis of methane release from gas hydrate system 43,44 . Gas hydrates are crystalline compounds that formed under high-pressure how temperature conditions, mainly containing water and CH 4 , with δ 13 C typically < − 60‰ 45 . Large volumes of gas hydrates are stored in continental slope and its size is sensitive to sea oor temperature change. In the event of global warming, enormous gas hydrates were converted to free methane gas. The massive amount of free 13 C-depleted methane would be able to drive the observed CIE. Permafrost thawing is another alternative/additional source of carbon 46 . A substantial amount of soil organic carbon was stored in the Antarctica and Circum-Arctic, and this terrestrial carbon is readily discharged when the temperature reached a certain threshold. The initial climate change might be caused by Earth's astronomical con gurations that favourable inducing warming, which is the original proposed mechanism of the methane hydrate dissociation hypothesis 44,47,48 . A weak causality from δ 13 C to δ 18 O (at 10%) has also been found in this interval. This can be interpreted by the positive feedback from the carbon release. Once a signi cant amount of methane was released, the increased concentration of CH 4 and CO 2 (by oxidization of CH 4 ) in the atmosphere would be able to further warm the ocean.
In the third quarter (Q3), only the causality relationship from δ 13 C to δ 18 O is signi cant (at 10%, Table 3).
This indicates that the carbon emission from submarine gas hydrates and permafrost may have reached its maximum. The geological process has shifted from a warming-driven carbon releasing dominated scenario to the scheme that carbon release is driving the global warming and environmental changes. This shift of driving mechanism might be caused by the positive feedback of the released carbon from gas hydrates destabilization and/or permafrost thawing, together with additional carbon injection from sources such as volcanic outgassing and thermogenic volatile and methane during the emplacement of the NAIP 49,50 , wild res burning peatlands 51 and oxidation of organic matter in shallow continental seaways 52 .
When coming into the fourth quarter (Q4), there is no causal relationship between the δ 13 C and δ 18 O data, despite that this interval is still at the 'body' of the CIE (Fig. 1, Table 3). It is possible that each pro le of δ 13 C and δ 18 O is progressing independently within their own cycle, or only minor causality exists that is not detected by this method. Over the whole investigated period, no causality exists for the δ 13 C and δ 18 O data (Fig. 1). This might result from the fact that there are multiple directions of causality within the whole period. We have noted several large gaps in the data due to the CaCO 3 dissolution 53 . However, the effect of the missing data is hard to evaluate.
The above discussion is based on assumption of stationary Q2 and Q3 interval series data. We acknowledge that some unit root test results do not suggest stationarity of the series, although the majority of unit root tests yielded stationarity (Tables 1 and 2). In the case of nonstationarity, these data would achieve stationary by taking the rst difference -I(1). The Toda-Yamashita method was then applied to investigate the causality for data in Q2 and Q3. Results suggest that the signi cant causality Loading [MathJax]/jax/output/CommonHTML/jax.js previously found on the level data no longer exist (Table 3). Nonetheless, a relatively small P-value is found from δ 18 O to δ 13 C for Q2 (0.1305) and from δ 13 C to δ 18 O for Q3 (0.1176) that have the same direction as previous Granger test results for the level data ( Table 2, 3). In addition, from Table 2's results, we found supporting evidence that δ 18 O and δ 13 C are stationary time series by using different unit-root testing methods (PP, DF-GLS and KPSS).

Quarters' boundaries and data size
Granger causality result suggests that the whole time series data of δ 13 C and δ 18 O show no signi cant causal relations, and that signi cant causality only emerges when pieces of the data are separated from each other into "quarters" (Fig. 1). The choice of these quarters is clearly in uenced by the authors' interpretation. We chose the quarters based on the behaviour of the carbon and oxygen isotope trends, i.e. before and after the PETM excursions (Q1 and Q4; Fig. 1). The PETM excursion is further divided into two parts (Q2 and Q3) by another obviously dramatic carbon isotope shift. However, even with the reference to the same criteria, location of the boundaries would be different by each observer and thus the result of the test would inevitably be affected.
In addition to the uncertainty regarding the quarters' boundaries, division of the dataset would also reduce the data size. We acknowledge that the data size is not large (281 for each series), compared with other subjects where this method is applied, like economic studies (e.g., research about the stock market). The fact that signi cance was only found in separated quarters rather than the whole series further indicates the notion of "P-hacking", i.e., if enough datasets are tested, a low P-value is inevitably found 54 . However, as previously stated, we divide the dataset based on background geological information about the PETM 26 and observable jumps in the scatter plot ( Fig. 1

Linearity and potential omitted variables
To apply the granger test, an underlying assumption is that the δ 13 C and δ 18  are probably not appropriately characterized as having a direct cause and effect relationship. If they are causally involved, there are certainly other variables that are omitted in the bivariate system, such as pCO 2 59 . To put this in a simpli ed PETM scenario, the δ 13 C shift was driven by addition of isotopicallylight carbon ( 13 C depleted) into the ocean and atmosphere -a process would unavoidably increase the pCO 2 level and the temperature. The carbon emission could also be caused by temperature rise as discussed previously. The oxygen isotope fractionation of the seawater (and thus the carbonates) is linked with temperature, and paleotemperature can be reconstructed with oxygen isotope data using some sort of equations (e.g., T (°C) = 16.5-4.3(δ 18 O CaCO3 -δ18O w-AMW ) + 0.14(δ 18 O CaCO3 -δ 18 O w-AMW ) 2 ) 60 . There are methods that can be used to estimate ancient pCO 2 levels, however, that would generally require other geochemical proxies (e.g. boron isotopes), and may involve large error with the estimation. Assuming linear relationship is not contradicting the true general co-movements between δ 13 C and δ 18 O, after controlling for structure breaks of the data. However, we would like to emphasise that the real geological scenario would be much more complicated than that as explained above.

Age-model
Granger causality is typically applied to time series whose timing is well known. Unlike recent monitored climate series data, climate series data recovered from drill cores involve some sort of uncertainty regarding the timing. A constant sedimentation rate (thus a linear depth-time translation) is obviously di cult to meet across this interval. According to nearby area sedimentation rates, the time gap is calculated to be about 20-25 years 12 . However, Zeebe et al., (2016) 14 estimated a sampling resolution of 40 years based on carbon cycle-climate modelling. Indeed, numerous studies have suggested highly variable sedimentation rates during the PETM as a response to hydrological perturbations 61 . However, our unequal spaced data is not affecting our causality's conclusion under linear's hypothesis, since our data has neither serious missing observations' problem nor high density of observations in the same period.

Conclusions
In summary, the Granger causality test suggests that the initial carbon release of the CIE might be caused by astronomical paced climate warming 8,46,47