M1 - Building and characterising a MDC record based on dated methane seepage samples

a) Building the record and testing its robustness

We compiled an original dataset of the documented occurrences of MDC up to 150 Ma by implementing similar compilations available in literature20,90. We counted as a single entry in the dataset any formation hosting MDC in a defined sedimentary basin. When a formation name is not available, MDC falling within an arbitrary 100 km radius within the same basin have been considered as a single entry. The minimum and maximum ages of dating for each carbonate sample have been determined either by radiometric dating (where available) or by considering the age of the hosting rock/sediment. The global MDC record we derived from it is a time series with a sampling time 1 My. Of 190 occurrences, 52 were dated in the last million years. Whilst these represent a sufficient dataset to observe variations in this time span with higher resolution sampling, our focus is on connecting modern seepage emission with those recorded across the longest available time history. We thus focused on the dataset of 138 samples dated between 1 Ma and 150 Ma. To build the record, we binned the ages for each sample in intervals of 1 My. For example, a sample dated between 2.6 and 3.5 Ma will be included in both the 2 Ma and 3 Ma age bins. A second sample dated between 3.5 and 5.2 Ma will be included in the 3 Ma, 4 Ma, and 5 Ma intervals. This binning strategy reduces the influence of samples with higher uncertainties (e.g., dated between 1 Ma and 10 Ma) on our results, so that the main variations will be due to samples dated exactly inside a time interval. The MDC occurrence record across the last 150 My is obtained by counting the occurrences in each interval (Fig. 2).

We estimated the robustness of the record for MDC based on the starting dataset, e.g., how independent it is of the particular selection of data. We tested the statistical robustness of the dataset using a bootstrapping approach (n value = 1000), as it is independent of the underlying statistical distribution of the data. We subtracted 30% of the MDC entries from the original dataset randomly and computed the normalised correlation coefficient between the original and the reduced time series. We repeated the procedure 1000 times and computed the average coefficient and average time-series over all iterations (Fig. 3a). We conclude that the dataset is stable based on the high average correlation coefficient (0.9605) and the minimal variance (twice the standard deviation = 0.0025) of the 30% bootstrap test (Fig. 3a). The correlation coefficient follows a normal distribution with these mean and standard deviation, as shown by a Kolmogorov-Smirnov test using a 5% significance threshold.

b) Assessing instability and removing biasing trends from the record

The MDC record shows important temporal trends that can be modelled (i.e. fitted) using polynomial and sinusoidal regressions. The samples that can be collected are likely more representative of modern than ancient times. Such a difference is likely causing the higher variability of MDC in modern times (Fig. 2). We performed linear, quadratic, cubic and sinusoidal regressions (Supplementary Fig. S3) and plotted the corresponding residuals: here, one or more quasi-periodic trends appear in all residuals, which are not stationary across the last 150 My. This inference is confirmed by a Durbin-Watson hypothesis test performed for each regression: the test outputs zero if the null hypothesis of stationarity in the residual is rejected, i.e., there is absence of autocorrelation among residuals at lag 1. The null hypothesis of absence of autocorrelation among the residuals is rejected at the 5% significance level in all cases (the p-value is 10−5 or lower). We selected a linear (decreasing) regression with time over the other regressions as the residual reduction for regression models of higher grade or complexity does not decrease sufficiently to justify their application91.

c) Spectral analysis, characterization and reconstruction of the MDC record from cyclicities

A strong perturbation may rule the deterministic behaviour of the linear residuals and create cyclicity; such signals can be recognised using matched filters, then calculating the spectrogram of the filtered signal92,93. As we do not have a reference to build a template-wavelet for matched filtering, we used the simplest approach, which is filtering the linear residuals using as filter coefficients the signal time-reversed93. We obtained the spectrogram of the filtered residuals with Hamming windows of 5 My and 60% overlap92 (Supplementary Fig. S4, threshold at 145 dB/My). These parameters allow the reconstruction of peaks while retaining stability. Two signals involve cyclicities with periods shorter than 4 My−1 (dashed white line – this is the limit cyclicity we decided to interpret given a 1 My sampling). They span the 6–21 My and 50–70 My intervals: the two intervals will be discussed in terms of its inherent cyclicities and amplitude variations. We observe that the most important cyclicity retrieved in these two periods is C2 = 1/12 My−1.

In addition to the matched filtering (Fig. S4), we studied cyclicities by computing the spectrograms of the MDC residuals (Supplementary Fig. S5a,b). The spectrograms are computed using Hamming windows of 10 My and 5 My (60% overlap, threshold at 140 dB/My−1) as the variations in windowing allow us to estimate lower and higher frequencies across the signal. For a 10 My window, we estimate a second median cyclicity (red circles – Supplementary Fig. S5b) as C1 = 1/26.66 My−1 ± 45%, where the uncertainty is given by the measurements’ variance. For a 5 My window the spectrograms confirm the onset of a median cyclicity of C2 = 1/12 My−1 ± 20% (Supplementary Fig. S5a) obtained with matched filtering. C1 and C2 stabilize during periods of strong MDC perturbations (50–70 Ma and 6–21 Ma) and become unstable, or disappear, for smaller MDC variations. We note that while C1 is more likely to dominate the entire 0–150 My time span given its lower sensitivity to shorter variations in the MDC controller, C2 is the clearest and most stable cyclicity after 49 My, as already observed using the matched filter. In the following analysis, we discuss the different spectral behaviours in the intervals 1–49 Ma, 50–70 Ma, and 71–150 Ma (Supplementary Fig. S5a, vertical dashed lines).

The full power spectrum of our MDC data across the last 150 My is obtained after filtering out all cyclicities of higher frequency than 1/50 My−1 (Supplementary Fig. S6, middle raw). This confirms that C1 and C2 dominate the MDC record. By calculating the power and phase spectra (Supplementary Fig. S6a,b) across the 1–70 Ma and 71–150 Ma time periods we observe that C1 (red vertical line) is effectively a stable cyclicity acting between 1 and 150 Ma, even if for the most recent period only two C1 cycles could be reconstructed. C2 (black line) is, in contrast, absent before 70 Ma. However, after this date, C2 becomes an important controller of MDC cyclicity (amplitude is > 1/2 of that of C1).

We want to see to what extent the entire MDC residual across the last 150 My can be reconstructed from the amplitudes and phases obtained using the spectral analysis only. The analytic signal derived from these parameters (frequencies, amplitudes and phases) reconstructs most of the filtered MDC residual variations (Fig. 4). The fit between C1 and MDC before 70 My is excellent, except for a transient in the period 130–140 Ma (r2 = 0.56). The agreement between the record and model worsens between 0 and 70 My ago (r2 = 0.12) comprising both a period of amplitude decrease (50–70 Ma) and one of strong variations (40–49 Ma). By adding cyclicity C2, the agreement between record and model is partially restored between 0 and 40 Ma (r2 = 0.36).

M2 - Stationary comparison with alternative observations

This work compares the MDC record with its potential geological driver(s) on the various intervals covered by the other investigated time series. We obtained representative time series of: seawater sulfate59, global sea level variation21,22,23, deep-sea temperature24, modelled Organic Carbon Burial (OCB)25, and sediment accumulation (SV)60. We re-sampled all the corresponding time series at 1 My and consider them as a unique ensemble.

We used Principal Component Analysis to discriminate meaningful dependency, i.e. observations that are (or are not) strictly correlated and add to the total variance. The analysis is on data after linear detrending. Three principal components account for most of the variability in the ensemble ( >85% - PC1 = 48%; PC2 = 20%; PC3 = 11%). The most important observations are (Fig. 3b; Supplementary Fig. S2):

(1) PC1 defines three groups of data (Figs. 3b and S2, left): (1) MDC and Sulfates (positively correlated); (2) OCB, SV, and sea level from Müller et al. (2008); (3) the two remaining sea levels and temperature (negatively correlated). More generally, SV only provides a significant contribution to PC3 (Supplementary Fig. S2, left). (2) Different sea level models21,22,23 correspond to different variations for PC1, which accounts for 48% of the variations in our dataset - the corrected Miller et al. 2005 curve shows the highest anticorrelation for both PC1 and PC2; (3) Sea level and deep-sea temperature are positively correlated (r2 = 0.48) and give similar contributions to PC1.

Considering the different contributions of the components to variability (PC3~1/5 PC1 and ~1/2 PC2) we decided to: (1) consider the corrected sea level curve only; (2) discard the sediment volumes from further analyses as they show correlations similar to those of OCB; (3) discard the temperature, which is worse-constrained than other parameters in time and highly correlated to the corrected sea level. While the corrected sea level shows significant anti-correlation with MDC (r2 = 0.34), sulfate and OCB show positive (but statistically insignificant - r2 of 0.12 and 0.09, respectively) correlation with the MDC.

We performed a second Principal Component Analysis to assess the dependence of our inferences on the number of observations included. This time, we only added to MDC the three different observations representing the main variability in the first analysis: MDC, the corrected sea level curve (anticorrelation for PCA), OCB and sulfate (positively correlated for PCA). The first two components account for more than 79% of the variance (Supplementary Fig. S2): PC1 = 55.0%; PC2 = 24.0%; PC3 = 11.87%. The PC1-PC2 plot (right) confirms the strong anticorrelation between MDC and sea level; however, from this graph there is no evidence of a preferential correlation between MDC and OCB or sulfate. This test confirms that the main controllers of MDC in the last 100 My have to be searched between sea level, sulfate and OCB.

The single oscillatory low-frequency signal observed for MDC also underlies the other three records (Fig. S7). These significant (r2 always > 0.4) signals have slightly different periods, (MDC – 74.26 My; sea level – 71.81 My; OCB – 62.5 My; sulfate – 69.93), with MDC having the same phase of OCB and sulfate, and the opposite phase of sea level. As we aim to obtain at least two cycles in the 100 My span, all the time series were thus filtered below 1/50 My−1. If this signal is removed, it is evident that the detrended sulfate residuals (Fig. S7) shows no cycle, a feature already apparent in the sulfate curve (Fig. 2): this low-period cyclicity in sulfate is due to their drastic increase especially between 50 and 40 My, the only non-linear variation in the curve in the last 100 My. The 40–50 My time span corresponds to the time of the anomalous transient in the analytic signal reconstructed from cyclicities only (Fig. 4). It thus highlights the central role of sulfate in triggering major changes in the cyclic behaviour of MDC. While the progressive increase in sulfate between 120 and 50 My might contribute to the strong changes in amplitude variations across the 50–70 My period (Figs S4–S5), the cyclic characteristics of MDC cannot be reconstructed in sulfate.

The power spectral analyses performed on the MDC were thus only applied to the sea level (in the period 0–140 Ma) and OCB (0–100 Ma) records. Supplementary Fig. S8 shows that, before 71 My ago, C1 is the only spike in the sea level spectrum and MDC records. The agreement between MDC and sea level spectra in this time span is particularly strong; we deduce that sea level change is the most significant low-frequency controller of MDC. On the other hand, C1 and C2 are the most relevant cyclicities after 70 Ma only for of MDC and OCB records. Among cyclic observations (thus excluding sulfate) the time correlation analysis described in the main text (Fig. 3c) ranks OCB as the only controllers showing instantaneous correlation. After 70 Ma ago, OCB saturates the MDC record and results a significant instantaneous controller of MDC emissions.