脱气速率对早侏罗世碳循环演变的主导作用

Degassing rate as the main driver of Early Jurassic carbon cycle evolution

摘要

碳循环是联结地球大气圈、水圈、岩石圈和生物圈物质与能量交换的重要过程,深刻影响着地球长期气候环境演变与生命演化。早侏罗世是潘基亚大陆裂解之后的中生代温室气候加剧的重要阶段,全球的碳循环格局在潜在的大陆裂解与大规模岩浆活动的驱动下发生了深度重塑。然而,对于早侏罗世 Hettangian 期至 Pliensbachian 期(约18 Ma)“背景期”碳循环演化特征的研究尚缺乏系统认识,关键难点在于如何通过高分辨率地球化学记录解析其背后的碳循环机制。为此,本研究采用表层地球化学演化模型(CHEES),结合贝叶斯马尔可夫蒙特卡洛(MCMC)反演方法,对来自北欧海道 Mochras 钻孔的高分辨率有机碳同位素记录进行了定量反演分析。研究结果表明:构造过程控制的脱气速率变化可能主导了Hettangian 期至Pliensbachian 期以来的碳循环长期趋势;早 Hettangian 期、Sinemurian-Pliensbachian 转换期及晚 Pliensbachian 期出现的多期次中等幅度的碳同位素负偏中可能存在额外碳释放源贡献。此外,脱气速率以及碳释放速率的反演结果存在405 ka长偏心率的短周期信号以及由偏心率幅度调制产生的约 6 Ma 的长周期信号,揭示早侏罗世碳释放过程可能受控于轨道周期和深部地幔过程或构造活动的节律性共同驱动。本研究丰富了对早侏罗世温室气候形成机制的认识,也为理解该时期全球气候变化与碳循环反馈机制提供了重要参考。

Abstract

The carbon cycle is a key process governing material and energy exchanges among Earth’s atmosphere, hydrosphere, lithosphere, and biosphere, and exerts a profound influence on long-term climate evolution and biological development. The Early Jurassic constitutes a critical stage of intensified Mesozoic greenhouse climate following the breakup of the supercontinent Pangaea. This period witnessed significant restructuring of the global carbon cycle, driven by continental rifting and large-scale magmatic activity. However, the evolution of the carbon cycle during the “background period” from the Hettangian to Pliensbachian stages (~18 Ma) remains poorly understood. A key challenge lies in deciphering the underlying carbon cycle mechanisms from high-resolution geochemical records. To address this, we applied the CHEES (CHemical Evolution of Earth Surface spheres) model coupled with a Bayesian Markov Chain Monte Carlo (MCMC) inversion approach to quantitatively reconstruct carbon cycle dynamics based on high-resolution organic carbon isotope (δ¹³C_org) records from the Mochras borehole in the northern Eurasian Seaway. The results indicate that variations in tectonically controlled degassing rates likely dominated the long-term carbon cycle trends from the Hettangian to Pliensbachian, whereas additional sources of carbon release may have contributed to multiple moderate-magnitude negative δ¹³C_org excursions observed during the Early Hettangian, the Sinemurian-Pliensbachian transition, and the Late Pliensbachian. Moreover, the inversion results clearly reveal a ~6 Ma periodicity modulated by short eccentricity maxima, as well as a 405 ka long eccentricity signal, suggesting that Early Jurassic carbon release processes were likely influenced by deep-seated mantle dynamics and tectonic cycles. This study provides new insights into the formation mechanisms of Early Jurassic greenhouse climates, the role of multi-sphere feedbacks, and the response of the carbon-cycle-climate system to tectonic and orbital forcing.

Graphical abstract

关键词

早侏罗世 / 碳循环 / 气候变化 / 有机碳埋藏

Key words

Early Jurassic / carbon cycle / climate changes / organic carbon burial

0 引言

地球系统由大气圈、水圈、岩石圈、生物圈等多个圈层组成,圈层之间的物质与能量交换循环维系着地球自形成以来的长期演化^[1]。在长地质时间尺度上,多圈层之间的相互作用不仅调控了地球演化过程中碳循环与气候系统的平衡,还驱动了关键生物地球化学过程的演变^[2-5]。例如,深部地幔通过构造活动、火山喷发向大气圈释放CO₂,改变辐射强度从而影响大气CO₂浓度和全球温度^[6];大陆硅酸盐风化消耗大气中的CO₂,通过河流搬运溶解性碳酸盐至海洋,以负反馈作用方式调节着大气与海洋之间的碳平衡^[7-8]。此外,生物圈通过初级生产固定CO₂,生成的有机碳在食物链传递及有机质沉降过程中,部分被埋藏于沉积物中,同样以负反馈机制减缓大气CO₂积累^[9]。这些复杂的圈层相互作用过程,在地质历史上驱动了多次全球性气候事件,系统解析其驱动机制与反馈关系是揭示温室气候形成与演化规律的关键途径。

早侏罗世是中生代温室气候加剧的关键阶段,从中大西洋大火成岩省开启了潘基亚大陆的裂解之后、伴随着潜在的洋中脊扩张加速或深部地幔活动增强,全球碳循环格局经历了深度重塑^[10]。现有研究主要关注于三叠纪-侏罗纪之交(Triassic-Jurassic Boundary event, TJB,约距今201 Ma)^[11-13]和Toarcian期海洋缺氧事件(Toarcian Oceanic Anoxic Event, T-OAE,约183 Ma)^[14-17]两个公认的碳循环极端事件,然而,TJB与T-OAE之间长达约18 Ma(约距今200~185 Ma)的“背景期”却存在明显研究空白。已有证据表明,这一阶段可能并非简单的气候稳定期,可能存在由轨道周期变化引起的缓慢而持续的碳脱气过程,对大气CO₂浓度和全球温度产生了长期影响^[18]。然而,关于其碳排放演变特征、规模及气候反馈效应的系统性研究仍较为缺乏,相关过程尚有待深入揭示。

地球系统箱式模型是研究长时间尺度碳循环变化的重要方法,其通过简化的海洋-大气-沉积物“储库”框架,能够高效整合大气、海洋、生物圈与岩石圈的碳储量与通量关系,定量模拟圈层间的物质交换过程及其反馈机制^[19-24]。相较于计算资源消耗大的耦合地球系统模式^[25],箱式模型凭借计算效率高、参数透明、易于调控等优势,能够有效模拟长时间尺度下碳循环-气候系统的动态演化,特别适用于探索驱动因子(如火山脱气、脱气速率变化)对古气候变化的定量影响。然而,现有模拟方法也面临局限性:一方面,缺乏高分辨率、可靠的古环境地质记录作为约束,导致模拟结果存在较大不确定性;另一方面,模型自身通常采用正演驱动的方式,难以解析脱气演变的时序细节,限制了对高频碳循环变化的解析能力。特别是对于早侏罗世这种长时间尺度、可能包含多阶段脱气脉冲的复杂过程,亟需在现有模型的基础上发展出能与更精细的地球化学记录结合的模拟方法,从而提升模拟的精度和可靠性。

基于此,本文拟围绕以下科学目标展开研究:(1)在箱式模型框架下,结合高分辨率碳同位素和沉积记录,定量约束整个早侏罗世的碳循环过程;(2)厘清早侏罗世碳排放过程的演变特征,揭示其长期变化模式及周期性节律;(3)探讨有机碳埋藏对碳排放变化的响应机制,明确其在碳循环调控中的作用。通过深入探讨上述问题,有望厘清早侏罗世这一长时间尺度以来碳循环演变过程,对于认识温室气候形成机制、圈层反馈效应及气候系统响应具有重要科学价值。

1 CHEES模型及反演方法

1.1 表层地球化学演化模型

本研究用了基于Qiu等^[26]和Zhao等^[24]所建立的6箱(5个海洋箱体+1个大气箱体)表层地球化学演化模型(CHemical Evolution of Earth Surface spheres, CHEES)来进行模拟。该模型一共包括5个海洋箱体,分别是浅大陆边缘海洋库(0~100 m水深,占据海洋表面的10% 和海底面积的5%)、深大陆边缘海洋库(100~1000 m水深,占据海底面积的另外5%)、表层海洋库(0~100 m水深,占据海洋表面的76.5%)、高纬度海洋库(0~250 m水深,占据海洋表面的13.5%)和深层海洋库(100~4 000 m水深,占据海底面积的90%)。模型的框架如图1a所示。海洋模块中考虑了混合和运输两个过程,混合发生在深部水体向浅部水体的扩散过程,运输则是由热盐环流驱动的运输流量,通常在高纬度表层和深层海洋之间或低纬度表层与其他海洋区之间发生。海洋箱体之间的混合流量由以下公式表示:

(1)mix_i,j=(C_i-C_j)*${\mathrm{f}}_{\mathit{m}\mathit{i}{\mathit{x}}_{\mathit{i},\mathit{j}}}$

其中,C_i和C_j是箱体i和j的物质浓度,f_mixi,j是箱体i和j之间的混合强度,即单位时间内由于扩散或湍流等物理过程导致的、从一个箱体传输到另一个箱体的单位体积或单位质量的物质交换速率。

热盐环流驱动的运输流量由以下公式表示:

(2)tran_i,j=C_i*${\mathrm{f}}_{\mathit{c}\mathit{i}\mathit{r}{\mathit{c}}_{\mathit{i},\mathit{j}}}$

其中,C_i是箱体i的物质浓度,f_circi,j是箱体i到j之间的热盐环流速度。

所有的海洋箱体都耦合到一个大气箱体中,并且在模型中,海气交换通量包括二氧化碳和氧气之间的交换,CO₂的海气交换量用以下公式计算:

(3)f_airseai=$\frac{{\mathit{A}}_{\mathit{i}}{\mathit{M}}_{\mathit{a}\mathit{t}\mathit{m}}}{{\mathit{\tau }}_{\mathit{o}\mathit{a}}}$(p(CO_2a)-p(CO_2i))

其中,A_i是每个海洋箱体的表面积,M_atm是现在大气中CO₂的摩尔量,τ_oa是海气交换的时间尺度,p(CO_2a)和p(CO_2i)是分别为大气和表层海洋箱体i中CO₂的分压。O₂的海气交换通量则通过以下公式表示:

(4)f_Oairseai=k_exchi(O_2eqi-O_2i)

其中,k_exchi是氧气的活塞速度,O_2eqi是氧饱和状态(单位为mol/m³),O_2i是表层海洋箱体i中的氧气浓度。

全球平均地表温度(GAST)计算公式^[27]如下:

(5)GAST=288+k_clim$\left(\frac{\mathit{l}\mathit{o}\mathit{g}\left(\frac{\mathit{C}{\mathit{O}}_{2}}{280}\right)}{\mathit{l}\mathit{o}\mathit{g}\left(2\right)}\right)$-7.4$\left(\frac{{\mathit{t}}_{\mathit{g}\mathit{e}\mathit{o}\mathit{l}}}{-570}\right)$

其中,k_clim表示长期气候敏感性,CO₂为大气中二氧化碳浓度(单位:10^-6),t_geol表示以距今时间(百万年)为单位的地质时间并以负值表示。大气中二氧化碳的体积比(10^-6)计算公式为

(6)CO₂=280$\frac{\mathit{C}{\mathit{O}}_{2\mathit{a}}}{\mathit{C}{\mathit{O}}_{2{\mathit{a}}_{0}}}$

此外,考虑到潜在的CO₂浓度增加引起的有机质分馏增加的情况,我们明确地将大气CO₂浓度对有机碳同位素分馏的影响纳入我们的碳循环模型中。具体来说,按照COPSE模型中使用的方法^[28],我们采用了一个将Δ_B与大气CO₂和O₂浓度变化联系起来的函数关系:

(7)Δ_B=${\mathrm{\Delta }}_{\mathit{B}}^{0}$+ε_B(CO₂)+ε(O₂)

其中,${\mathrm{\Delta }}_{\mathit{B}}^{0}$表示规定的碳同位素分馏上限,本研究在这里设置为33。ε_B(CO₂)表示由p(CO₂)变化驱动的有机碳同位素分馏变化。这种分馏反应可由下式经验描述:

(8)ε_B(CO₂)=-$\frac{9}{\sqrt{\mathit{p}\left(\mathit{C}\right.{\mathit{O}}_{2\mathit{a}})}}$

p(CO_2a)表示归一化大气CO₂浓度,以PAL (Pre-Industrial Atmospheric Level)为单位表示。大气氧浓度对碳同位素分馏的影响用Berner等^[29]提出的线性关系表示:

(9)ε(O₂)=J(o-1)

其中,o表示归一化后的大气O₂浓度(以PAL为单位表示),J为曲线拟合参数(在COPSE 定义为5)。

模型中设置了风化和初级生产力两个主要的过程来移除大气-海洋库中的碳。风化模块包括了碳酸盐风化和硅酸盐风化(玄武岩和花岗岩),并且二者的风化强度与全球平均温度直接相关。碳酸盐风化、玄武岩和花岗岩与温度公式用Zhao等^[24]的表达,分别如下:

(10)f_Tcarb=1+0.087(GAST-288)

(11)f_Tbas=e^{0.060 8(GAST-288)}

(12)$\begin{array}{c} (1+0.038(\text { GAST }-288))^{0.65} \\ f_{\text {Tgran }}=\mathrm{e}^{0.0724(\text { GAST }-288)} \\ (1+0.038(\text { GAST }-288))^{0.65} \end{array} $

初级生产力模块模拟了表层海洋箱体中光合作用固定CO₂的过程,并考虑了营养盐(如磷酸盐)和微量元素(如铁)的协同限制作用。初级生产力的计算公式为

(13)Pripi_i=k_ppi*min$\left[{\mathit{C}}_{\mathit{D}{\mathit{P}}_{\mathit{i}}}\mathit{*}\frac{{\mathit{C}}_{\mathit{D}{\mathit{P}}_{\mathit{i}}}}{{\mathit{C}}_{\mathit{D}{\mathit{P}}_{\mathit{i}}}+{\mathit{K}}_{\mathit{p}}}\mathit{*}\frac{{\mathit{C}}_{\mathit{F}\mathit{e}{\mathit{Ⅲ}}_{\mathit{i}}}}{{\mathit{\gamma }}_{\mathit{F}\mathit{e}\mathit{P}}}\mathit{*}\frac{{\mathit{C}}_{\mathit{F}\mathit{e}{\mathit{Ⅲ}}_{\mathit{i}}}}{{\mathit{C}}_{\mathit{F}\mathit{e}{\mathit{Ⅲ}}_{\mathit{i}}}+{\mathit{K}}_{\mathit{F}\mathit{e}}}\right]$

初级生产力受到磷酸盐和铁两者中较低可用性的控制,呈现出协同限制的特性。在模型中,表层初级生产力是大气CO₂向生物碳库转化的核心过程,其产物颗粒有机碳部分沉降到深海,为有机碳埋藏和深层营养盐再循环提供基础。

因此,在考虑到有外源碳输入情况下的t时刻的大气CO₂浓度和大气中CO₂的碳同位素、以及5个海洋库中有机碳同位素的通量方程如下:

(14)$\frac{\mathit{d}\left(\mathit{C}{\mathit{O}}_{2\mathit{a}}\right)}{\mathit{d}\mathit{t}}$=-f_airseasm-f_airseas-f_airseah+f_ccdeg+f_ocdeg+f_oxidw-f_locb-f_carbw-2f_silw+f_CO2input

(15)$\frac{\mathit{d}\left({\mathit{\delta }}^{13}{\mathit{C}}_{\mathit{C}{\mathit{O}}_{2\mathit{a}}}\right)}{\mathit{d}\mathit{t}}$=-${\mathrm{f}}_{\mathit{a}\mathit{i}\mathit{r}\mathit{s}\mathit{e}{\mathit{a}}_{\mathit{s}\mathit{m}}}^{13}$-${\mathrm{f}}_{\mathit{a}\mathit{i}\mathit{r}\mathit{s}\mathit{e}{\mathit{a}}_{\mathit{s}}}^{13}$-${\mathrm{f}}_{\mathit{a}\mathit{i}\mathit{r}\mathit{s}\mathit{e}{\mathit{a}}_{\mathit{h}}}^{13}$+$\frac{{\mathit{f}}_{\mathit{c}\mathit{c}\mathit{d}\mathit{e}\mathit{g}}}{1+\frac{1}{{\mathit{R}}_{\mathit{C}}}}$+$\frac{{\mathit{f}}_{\mathit{o}\mathit{c}\mathit{d}\mathit{e}\mathit{g}}}{1+\frac{1}{{\mathit{R}}_{\mathit{G}}}}$+$\frac{{\mathit{f}}_{\mathit{o}\mathit{x}\mathit{i}\mathit{d}\mathit{w}}}{1+\frac{1}{{\mathit{R}}_{\mathit{G}}}}$-${\mathrm{f}}_{\mathit{l}\mathit{o}\mathit{c}\mathit{b}}^{13}$-f_carbw*$\frac{{\mathit{\delta }}^{13}{\mathit{C}}_{\mathit{C}{\mathit{O}}_{2\mathit{a}}}}{\mathit{C}{\mathit{O}}_{2\mathit{a}}}$-2f_silw*$\frac{{\mathit{\delta }}^{13}{\mathit{C}}_{\mathit{C}{\mathit{O}}_{2\mathit{a}}}}{\mathit{C}{\mathit{O}}_{2\mathit{a}}}$+f_CO2input*δ¹³C_CO2input

(16)$\frac{\mathit{d}\left({\mathit{\delta }}^{13}{\mathit{C}}_{\mathit{P}\mathit{O}{\mathit{C}}_{\mathit{s}\mathit{m}}}\right)}{\mathit{d}\mathit{t}}$=Pripr_sm*(δ¹³C_DICsm-Δ_B)-Remin_sm*$\frac{{\mathit{\delta }}^{13}{\mathit{C}}_{\mathit{P}\mathit{O}{\mathit{C}}_{\mathit{s}\mathit{m}}}}{\mathit{P}\mathit{O}{\mathit{C}}_{\mathit{s}\mathit{m}}}$-Tra${\mathrm{n}}_{\mathit{P}\mathit{O}{\mathit{C}}_{\mathit{s}\mathit{m}\_\mathit{s}}}$*$\frac{{\mathit{\delta }}^{13}{\mathit{C}}_{\mathit{P}\mathit{O}{\mathit{C}}_{\mathit{s}\mathit{m}}}}{\mathit{P}\mathit{O}{\mathit{C}}_{\mathit{s}\mathit{m}}}$-Tra${\mathrm{n}}_{\mathit{P}\mathit{O}{\mathit{C}}_{\mathit{s}\mathit{m}\_\mathit{d}\mathit{m}}}$*$\frac{{\mathit{\delta }}^{13}{\mathit{C}}_{\mathit{P}\mathit{O}{\mathit{C}}_{\mathit{s}\mathit{m}}}}{\mathit{P}\mathit{O}{\mathit{C}}_{\mathit{s}\mathit{m}}}$-Sed_sm*$\frac{{\mathit{\delta }}^{13}{\mathit{C}}_{\mathit{P}\mathit{O}{\mathit{C}}_{\mathit{s}\mathit{m}}}}{\mathit{P}\mathit{O}{\mathit{C}}_{\mathit{s}\mathit{m}}}$

(17)$\frac{\mathit{d}\left({\mathit{\delta }}^{13}{\mathit{C}}_{\mathit{P}\mathit{O}{\mathit{C}}_{\mathit{d}\mathit{m}}}\right)}{\mathit{d}\mathit{t}}$=-Remin_dm*$\frac{{\mathit{\delta }}^{13}{\mathit{C}}_{\mathit{P}\mathit{O}{\mathit{C}}_{\mathit{d}\mathit{m}}}}{\mathit{P}\mathit{O}{\mathit{C}}_{\mathit{d}\mathit{m}}}$+Tra${\mathrm{n}}_{\mathit{P}\mathit{O}{\mathit{C}}_{\mathit{d}\_\mathit{d}\mathit{m}}}$*$\frac{{\mathit{\delta }}^{13}{\mathit{C}}_{\mathit{P}\mathit{O}{\mathit{C}}_{\mathit{d}}}}{\mathit{P}\mathit{O}{\mathit{C}}_{\mathit{d}}}$+Tra${\mathrm{n}}_{\mathit{P}\mathit{O}{\mathit{C}}_{\mathit{s}\mathit{m}\_\mathit{d}\mathit{m}}}$*$\frac{{\mathit{\delta }}^{13}{\mathit{C}}_{\mathit{P}\mathit{O}{\mathit{C}}_{\mathit{s}\mathit{m}}}}{\mathit{P}\mathit{O}{\mathit{C}}_{\mathit{s}\mathit{m}}}$-Sed_dm*$\frac{{\mathit{\delta }}^{13}{\mathit{C}}_{\mathit{P}\mathit{O}{\mathit{C}}_{\mathit{d}\mathit{m}}}}{\mathit{P}\mathit{O}{\mathit{C}}_{\mathit{d}\mathit{m}}}$

(18)$\frac{\mathit{d}\left({\mathit{\delta }}^{13}{\mathit{C}}_{\mathit{P}\mathit{O}{\mathit{C}}_{\mathit{s}}}\right)}{\mathit{d}\mathit{t}}$=Pripr_s*(δ¹³C_DICs-Δ_B)-Remin_s*$\frac{{\mathit{\delta }}^{13}{\mathit{C}}_{\mathit{P}\mathit{O}{\mathit{C}}_{\mathit{s}}}}{\mathit{P}\mathit{O}{\mathit{C}}_{\mathit{s}}}$+Tran_POCsms*$\frac{{\mathit{\delta }}^{13}{\mathit{C}}_{\mathit{P}\mathit{O}{\mathit{C}}_{\mathit{s}\mathit{m}}}}{\mathit{P}\mathit{O}{\mathit{C}}_{\mathit{s}\mathit{m}}}$-Tran_POCsh*$\frac{{\mathit{\delta }}^{13}{\mathit{C}}_{\mathit{P}\mathit{O}{\mathit{C}}_{\mathit{s}}}}{\mathit{P}\mathit{O}{\mathit{C}}_{\mathit{s}}}$-Sink_POCsd*$\frac{{\mathit{\delta }}^{13}{\mathit{C}}_{\mathit{P}\mathit{O}{\mathit{C}}_{\mathit{s}}}}{\mathit{P}\mathit{O}{\mathit{C}}_{\mathit{s}}}$

(19)$\frac{\mathit{d}\left({\mathit{\delta }}^{13}{\mathit{C}}_{\mathit{P}\mathit{O}{\mathit{C}}_{\mathit{h}}}\right)}{\mathit{d}\mathit{t}}$=Pripr_h*(δ¹³C_DICh-Δ_B)-Remin_h*$\frac{{\mathit{\delta }}^{13}{\mathit{C}}_{\mathit{P}\mathit{O}{\mathit{C}}_{\mathit{h}}}}{\mathit{P}\mathit{O}{\mathit{C}}_{\mathit{h}}}$+Tran_POCsh*$\frac{{\mathit{\delta }}^{13}{\mathit{C}}_{\mathit{P}\mathit{O}{\mathit{C}}_{\mathit{s}}}}{\mathit{P}\mathit{O}{\mathit{C}}_{\mathit{s}}}$-Sink_POChd*$\frac{{\mathit{\delta }}^{13}{\mathit{C}}_{\mathit{P}\mathit{O}{\mathit{C}}_{\mathit{h}}}}{\mathit{P}\mathit{O}{\mathit{C}}_{\mathit{h}}}$

(20)$\frac{\mathit{d}\left({\mathit{\delta }}^{13}{\mathit{C}}_{\mathit{P}\mathit{O}{\mathit{C}}_{\mathit{d}}}\right)}{\mathit{d}\mathit{t}}$=-Remin_POCd*$\frac{{\mathit{\delta }}^{13}{\mathit{C}}_{\mathit{P}\mathit{O}{\mathit{C}}_{\mathit{d}}}}{\mathit{P}\mathit{O}{\mathit{C}}_{\mathit{d}}}$+Sink_POCsd*$\frac{{\mathit{\delta }}^{13}{\mathit{C}}_{\mathit{P}\mathit{O}{\mathit{C}}_{\mathit{s}}}}{\mathit{P}\mathit{O}{\mathit{C}}_{\mathit{s}}}$+Sink_POChd*$\frac{{\mathit{\delta }}^{13}{\mathit{C}}_{\mathit{P}\mathit{O}{\mathit{C}}_{\mathit{h}}}}{\mathit{P}\mathit{O}{\mathit{C}}_{\mathit{h}}}$-Tran_POCddm*$\frac{{\mathit{\delta }}^{13}{\mathit{C}}_{\mathit{P}\mathit{O}{\mathit{C}}_{\mathit{d}}}}{\mathit{P}\mathit{O}{\mathit{C}}_{\mathit{d}}}$-Sed_d*$\frac{{\mathit{\delta }}^{13}{\mathit{C}}_{\mathit{P}\mathit{O}{\mathit{C}}_{\mathit{d}}}}{\mathit{P}\mathit{O}{\mathit{C}}_{\mathit{d}}}$

1.2 反演方法及模拟实验设计

本研究在 Qiu等^[26]提出的马尔可夫蒙特卡洛(MCMC)贝叶斯时序递进式反演模型框架下开展模拟,其方法具备全局优化能力强、对复杂模型适应性高、对初始值不敏感、具有明确概率解释以及动态调整机制等优势,已在早 Toarcian 期多参数联合反演中应用,能够在还原长期趋势的前提下同时有效捕捉高频变化特征。本研究在现有框架下,调整了用于约束模型结果的观测参数,选择仅采用有机碳碳同位素(δ¹³C_org)序列作为约束指标,旨在反演早侏罗世背景期碳循环演化过程及碳排放节奏。反演流程包含6个模块(图1b):“敏感性分析”、“参数更新”、“参数提议”、“模型执行”、“似然计算”和“适应性调整”。在“敏感性分析”模块中,设定了待反演参数的取值范围,采用蒙特卡洛采样输入模型初始化模块。针对本研究目标,重点考虑了气候敏感性参数(设定区间为2.5~6 ℃)。在“参数提议”模块中,初始参数集(ϕ₀)通过时间序列扰动生成新候选参数集(ϕ),模拟退火式递进优化;在“模型执行”模块中,分别运行标准模拟与扰动模拟,获得两组结果。在“似然计算”模块中,针对 δ¹³C 数据计算对数似然度,用以下公式评价模型拟合效果:

(21)l${\mathrm{l}}_{{\mathit{\delta }}^{13}\mathit{C}}$=logpdf(Normal(δ¹³C_observed,δ¹³C_std),δ¹³C_model)

通过迭代更新参数集,提升拟合度,并通过“适应性调整”模块动态优化扰动幅度,平衡全局搜索与局部优化能力。通过该方法,本文基于δ¹³C 观测记录,对早侏罗世背景期碳排放过程及碳循环反馈演化进行了定量反演与约束。

为规避反演过程中陷入局部最优解的风险,我们在本研究中采用了双重策略。首先,在迭代过程中设计了一个由全局到局部的多尺度收敛机制:第一阶段将目标曲线划分为 5 个区段以捕捉全局趋势;第二阶段细分为 15 个区段进行中等尺度优化;第三阶段则基于目标曲线上相邻两点构成的微小区间进行高精度模拟以恢复细节。该机制不仅有效降低了陷入局部最优的可能性,同时加速了模型向全局最优解的逼近。其次,对于所设置的每组反演模拟实验,均执行了超过200次的独立蒙特卡洛模拟,每次均采用随机初始模型参数以确保收敛过程彼此独立,从而显著降低了因单次初始化不良而陷入局部解的概率。如图2所示的模拟实验一中的误差收敛曲线,每条彩色曲线代表一次独立模拟的 Log-Likelihood 收敛过程。结果显示,模型在初始数次迭代后即迅速捕捉到目标曲线的宏观变化趋势(对应于收敛策略的前两个阶段),此后的迭代则集中于细节特征的恢复(对应于第三阶段),表现为 Log-Likelihood 值逐渐趋于稳定并收敛。不同模拟结果之间虽存在一定差异,但整体收敛趋势一致,验证了所提出模型策略的有效性和稳健性。

2 讨论

碳同位素记录被普遍视为记录大气-海洋系统碳循环扰动的重要地球化学指标,能够敏感反映外源碳输入事件及碳库变化过程。本研究选用2020年Storm等^[18]报道的北欧海道 Mochras 钻孔(Mochras borehole)中自晚三叠世 Rhaetian 期至早侏罗世 Toarcian 期、跨度约 18 Ma (约距今201~185.5 Ma)的高分辨率有机碳同位素(δ¹³C_org)数据作为反演目标曲线。该数据集目前被认为是该时段内相对可靠、分辨率最高的有机碳同位素序列。Mochras 剖面中观察到的脂屑镜质组呈现出高度降解的形态以及颗粒的荧光特性指示了有机质来源于海洋藻类。此外,Storm等^[18]将 Mochras 钻孔的有机碳同位素曲线与 Paris 盆地 Sancerre-Couy 钻孔^[30](上Hettangian统至Pliensbachian阶)的无机碳同位素记录及 Bristol Channel 盆地^[31](下Hettangian统 至下 Sinemurian统)无机碳同位素记录进行对比,结果显示 Mochras 有机碳同位素的波动具有全球一致性,反映全球碳循环扰动信号并非区域效应。这为本研究基于箱式模型从全球尺度探讨碳循环演变提供了坚实的数据支持。为探讨该δ¹³C_org曲线背后的环境驱动机制,本研究设计了四组模拟实验,针对不同碳释放情景进行反演分析:模拟一设定外源碳输入量为0,允许脱气速率在每个时间步动态变化。通过拟合模拟δ¹³C_org曲线至目标曲线,反演得到整个模拟时段内脱气速率的演变过程。模拟二固定脱气速率,设定额外碳释放端员同位素组成为 -5‰,允许碳释放量在各时间步变化,旨在反演模拟来自火山碳源的释放通量变化特征。模拟三与模拟四与模拟二设计一致,固定脱气速率并反演外源碳释放总量,差异在于假设不同碳源端员同位素组成;并且模拟三中设为 -25‰,模拟热成因甲烷释放过程;模拟四中设为-60‰,模拟生物成因甲烷释放场景。通过上述多组情景模拟,结合贝叶斯反演框架,旨在定量约束早侏罗世从距今200 Ma开始18 Ma以来的δ¹³C_org曲线所记录的碳释放过程及其对碳循环动态演化的驱动效应。

四组反演实验结果均较好地再现了目标δ¹³C_org曲线(图2a,3a)中观测到的主要变化趋势及高频波动特征。模拟 δ¹³C_org曲线成功捕捉了 Hettangian 至中 Sinemurian 期间约 4‰ 的正偏移,中 Sinemurian 至中 Pliensbachian 期间持续约 4 Ma、幅度约 5‰ 的负偏移及其后恢复阶段,以及晚 Pliensbachian 期快速负偏移事件(幅度约 4‰)。值得注意的是,在设定外源碳输入量为零、仅允许脱气速率变化的模拟一中,模型在早 Sinemurian 期 δ¹³C_org正偏移阶段的拟合效果较差;而在引入外源碳输入的模拟二、三、四中则不存在该问题,整体拟合效果明显提升。这一差异主要源于所选用的δ¹³C_org数据受有机碳同位素分馏效应的影响,而该分馏过程在模型中受到大气CO₂浓度变化的调控。

具体而言,大气 CO₂浓度及温度变化可能会影响碳在大气-海洋系统中不同碳库(如$\mathrm{HCO}_{3}^{-}$、$\mathrm{CO}_{3}^{2-}$)与沉积相(如碳酸盐、有机分子)之间的同位素分馏行为^[32]。类似于现代植物或浮游植物(如球石藻)对碳源的偏好性,其¹²C的优先利用程度可能会随CO₂浓度和温度的变化而动态调整^[33]。在模拟一中,由于岩浆脱气碳源的同位素组成较重(约 0‰),为重现 Hettangian 期内多次 δ¹³C_org负偏移事件,模型需模拟出较大的碳释放量,导致大气CO₂浓度在正偏阶段之前已升高至相对较高水平(接近 4 000×10^-6)。进入正偏阶段后,尽管脱气速率已大幅降低,但因先前累积的高CO₂浓度难以快速降低,导致 δ¹³C_org模拟值难以达到观测曲线中的高正值。此外,为维持模型平衡,脱气速率在该阶段已接近设定下限,进一步限制了拟合能力。相比之下,模拟二、三、四中所释放碳端员同位素值更负(分别为 -5‰、-25‰ 和 -60‰),使得正偏移前期大气CO₂浓度处于较低水平,正偏阶段碳释放量趋近于零。在碳移除机制作用下,CO₂浓度进一步降低,有利于完整重现正偏阶段δ¹³C_org曲线。此外,碳端员同位素值越负,正偏阶段δ¹³C_org拟合效果越优。另一个显著特征表现于 Sinemurian-Pliensbachian 界限的负偏移及晚 Pliensbachian 期负偏移阶段。尽管模拟一亦能模拟出这两个负偏阶段,但对应的有机碳埋藏速率和初级生产力已接近模型设定的上限;而模拟二、三、四中上述指标保持在合理范围,且碳端员同位素越负,对应的有机碳埋藏和初级生产力响应越低。这是因为,碳端员越负,模拟δ¹³C_org所需的碳释放量越小,相应引起的大气CO₂浓度波动及增温幅度亦更有限,进而对磷风化输入、海洋初级生产力提升及有机碳埋藏的促进作用亦相对减弱。

厘清地质历史演化过程中碳释放的来源,有助于深入理解古气候变化的驱动机制及碳循环-气候系统的耦合过程,但这是一个关键、复杂且具有挑战性的问题。由于在早侏罗世约 18 Ma 的记录中尚缺乏其他可靠指标(如温度、大气 CO₂浓度)可用于反演约束,因而本研究不同于Qiu等^[26]研究中基于多指标联合反演CO₂释放量、甲烷释放量及来源的方法^[26],本次反演只能采用一个“已知数”(有机碳同位素)来求解“一个未知数”(碳来源)的方法。虽然这一方法无法精确解析早侏罗世各阶段碳释放量及其具体来源,但通过设计的四组模拟实验,仍可对碳释放演化趋势提供一定的定量约束。在长周期尺度上,构造过程主控的脱气速率变化仍是早侏罗世碳释放的主要来源,模拟结果表明脱气速率变化可解释δ¹³C_org曲线中大部分长期趋势。然而,在若干关键时期,如早 Hettangian 期、Sinemurian-Pliensbachian 转换期及晚 Pliensbachian 期,δ¹³C_org的显著负偏移特征表明极可能存在额外碳释放源的贡献。对于早 Hettangian 期,δ¹³C_org曲线在整体正偏趋势中叠加有脉冲式负偏移,可能反映出晚三叠生物大灭绝后不稳定恢复阶段的碳循环扰动。此类快速脉冲负偏移特征可由少量生物成因甲烷释放驱动,由于甲烷同位素端员极负(约 -60‰),即使少量释放亦可引起显著δ¹³C_org扰动。值得强调的是,尽管中大西洋大火成岩省(CAMP)活动和潘吉亚裂解过程本身持续时间较长(百万年尺度^[10],但相关构造活动可能在早 Hettangian 期阶段性增强,诱发短时效的环境扰动。例如,沿大陆边缘的断裂和喷发可能在特定时间窗内集中释放大量甲烷水合物,从而驱动脉冲性碳同位素负偏移。也就是说,CAMP 和潘吉亚裂解提供了长期背景条件,而早 Hettangian 的快速负偏移则反映了其中某一特定时段的强烈响应。Sinemurian-Pliensbachian 转换期δ¹³C_org曲线表现出持续约4 Ma的长时间尺度负偏移,该特征不太可能由生物成因甲烷长期释放驱动,因为缺乏合理的长期甲烷形成、储存与释放机制。虽然之前有种解释认为在气候冷却趋势下,伴随¹²C 富集的深层水体上涌或循环,是作为驱动Sinemurian-Pliensbachian 转换期的另一种可能机制^[34-35],但更可能的解释是构造脱气速率增强叠加火山CO₂释放,共同驱动δ¹³C_org的缓慢负偏移过程。晚 Pliensbachian 期的正偏移通常被认为与高海平面时期和温暖气候条件下广泛沉积同位素偏轻的有机质相关,目前尚未发现其与大规模火山活动存在明确的时间对应关系^[36]。随后发生的气候冷却、海平面下降以及水体混合恢复,被认为通过沉积物再搬运、氧化和异养再矿化作用释放了积累的同位素偏轻的碳,进而导致晚 Pliensbachian 期和早Toarcian期界限出现的负偏移^[34,37-38]。此外,来自Hg元素富集的证据表示,在 Pliensbachian -Toarcian期界限出现的负偏移与Hg浓度和Hg/TOC含量比值的增加对应良好,表明了可能也与Karoo LIP 的启动有关^[39]。尽管有学者提出Karoo LIP 岩浆侵入Karoo盆地二叠系煤层,可能通过热解作用释放热成因甲烷^[40-41],但已有煤层样品分析表明,距岩浆侵入体不同距离处煤层δ¹³C_org无显著变化,且镜质组反射率和脂质组含量均偏低,表明热成因甲烷释放量有限^[42]。综上所示,尽管受限于单指标反演,模拟结果仍揭示出早侏罗世碳释放过程存在明显的阶段性差异,构造脱气、火山CO₂和生物成因甲烷可能在不同阶段占据不同主导地位。

4组模拟结果显示,有机碳埋藏速率的高值主要对应于δ¹³C_org负偏移阶段,而正偏移阶段埋藏速率相对较低(图3和图4)。这一结果表明,早侏罗世δ¹³C_org变化并非简单由有机碳埋藏驱动,相反的,有机碳埋藏是作为响应脱气速率变化的一种反馈调节机制。具体而言,负偏阶段受控于大规模脱气和外源碳释放,导致大气CO₂浓度上升与增温效应,进而通过加速磷输入和促进海洋初级生产力显著增强有机碳埋藏。然而,由于外源输入通量远大于生物泵反馈作用,δ¹³C_org仍表现为负偏移。随着碳释放终止,脱气速率降低,大气CO₂浓度下降,δ¹³C_org转为正偏,此时有机碳埋藏虽仍维持一定水平,但相较于负偏阶段已显著降低。此外,需要特别指出的是,本研究所指的δ¹³C_org正偏,实质上更可能是负偏之后的恢复过程,即δ¹³C_org值向基值回升,而非超过基值的绝对正偏。因此,这一正偏与地质历史上由有机碳埋藏增多引起的正偏不同,主要反映脱气速率降低带来的外源碳输入减少所致的碳同位素恢复过程。综上,早侏罗世δ¹³C_org的主控机制为外源碳输入速率的变化,有机碳埋藏对其演变趋势具有重要但次要的缓冲调节作用,而非决定性主导。这一认识与 “埋藏主控假说”存在区别,强调了构造因素导致的脱气速率变化在地质历史时期碳循环中的主导地位。

我们对脱气速率以及火山CO₂的释放速率进行了频谱分析,结果显示除了稳定且显著的405 ka长偏心率周期外,还存在显著的约6 Ma周期信号(图5),其来源尚未明晰。应用最广泛的理论天文解如La2004中记录的最长尺度周期为4.5 Ma偏心率调制周期。尽管4.5 Ma周期在奥陶纪—志留纪碳和氧同位素序列中被识别^[43],但其理论信号强度远远弱于2.4 Ma周期,其在地层中的保存机制仍然存在疑问。超过4.5 Ma的长尺度周期未被记录在理论天文解中,这些周期的来源解释及驱动机制是近年来地球科学领域的研究热点之一^[44]。一个可能的解释是多圈层相互作用,即轨道周期作用与深部过程的耦合作用放大了理论上较弱的长尺度周期,并产生了一系列超长尺度调幅周期^[45]。约36 Ma周期和约9 Ma周期得到了广泛的研究^[46-52]。已有研究表明由地幔对流与俯冲板片相互作用驱动的36 Ma地幔-岩石圈深水循环周期调控了约36 Ma构造-海平面波动,进而对约36 Ma生物多样性波动产生驱动作用^[48,50-51]。巧合的是,太阳-行星系统垂向振荡过程中,地球穿越银道面的周期约为约(32±3) Ma,这意味着约36 Ma构造-气候-生物协同变化可能受地球在太阳系及银河系中的运动调控^[47,50-51]。类似案例也出现于约9 Ma周期尺度。约9 Ma周期最初被认为源于构造作用,但中-新生代La2004偏心率振幅调制包络线频谱分析显示其存在天文信号特征,可能为长尺度偏心率周期^[49,52];而中-新生代碳同位素记录与中生代生物硅序列也呈现约9 Ma准周期,表明约9 Ma周期对气候波动具有调制作用^[52-53]。相比于上述长尺度周期以及熟知的约2.4 Ma偏心率周期和约1.2 Ma斜率周期,约6 Ma周期研究较少。本研究给出了更直观的证据,揭示了受控于深部地质活动的脱气速率变化以及火山CO₂的释放速率存在显著的约6 Ma周期,其可能与轨道周期存在因果/协同关系,从而进一步调控了表生环境变化。这说明天文-构造-气候-生物协同变化可能从约36 Ma周期尺度拓展至本次研究提出的约6 Ma周期的更小尺度。综上,在小于36 Ma的时间尺度上,约6 Ma与约9 Ma周期均可能存在天文-构造-表生环境的耦合关系。尽管这些较小尺度周期是否与银河系动力学运动相关仍需进一步验证,但本研究证实:6 Ma周期尺度上存在天文-构造-气候的多圈层相互作用,该作用通过调控碳循环波动,促进有机碳的周期性埋藏与亏损。

尽管本研究基于高分辨率δ¹³C_org数据初步揭示了早侏罗世碳释放过程的阶段性特征,但当前反演框架仍存在一定局限,主要体现在参数约束单一,尚不足以全面还原碳循环-气候系统中多圈层之间的复杂反馈机制。未来研究可从以下几个方面进一步拓展与深化:一方面,结合温度指标、大气CO₂浓度重建结果、金属元素及其稳定同位素等多种地球化学记录,构建多参数联合贝叶斯反演框架,有望提升反演结果的精度与可靠性,并进一步解析不同碳源(如岩浆碳、生物成因甲烷)与碳汇(特别是有机碳埋藏)在碳循环演化中的相对贡献;另一方面,模型层面可引入更加复杂的圈层耦合过程,例如考虑岩石圈-大气圈-海洋圈之间的物质通量、深部碳储库输入(如板块俯冲、洋中脊扩张)与表层碳输出之间的动态平衡过程,提升模型对地质时期多圈层联动机制的模拟能力。此外,有机碳埋藏是地球碳循环中最关键的长期碳汇之一,连接生物圈的初级生产与岩石圈的碳储存过程,在维持大气 CO₂ 浓度稳态、调节地表温度、促进氧气积累及支撑生态系统演替中发挥了不可替代的作用。作为圈层耦合的重要纽带,有机碳埋藏同时受到气候变迁、海洋生产力、沉积环境与构造活动等多重因子的影响,其变化直接影响地球系统对碳扰动事件的响应与缓冲能力。未来研究需进一步聚焦不同地质时期,特别是在温室-冰室转变、大洋缺氧事件及生物大灭绝等关键节点上,有机碳埋藏通量的时空演变规律与驱动机制,借助高分辨率地球化学记录与定量耦合模型,系统解析其反馈路径及相对贡献。这不仅有助于厘清显生宙以来碳汇动态如何塑造全球气候-生态格局,也为预测未来多圈层响应路径提供关键的古气候视角与理论支持。

3 结论

(1)本研究基于高分辨率 Mochras 钻孔δ¹³C_org记录,采用贝叶斯马尔可夫蒙特卡洛反演方法,结合CHEES箱式模型对早侏罗世约 18 Ma 时段内的碳循环演化过程进行了定量反演,通过设计四组不同碳释放情景的模拟实验,系统探讨了构造脱气、火山CO₂、生物成因甲烷和热成因甲烷释放在δ¹³C_org波动中的可能贡献机制。

(2)早侏罗世Hettangian 期以来至Pliensbachian 期长达约18 Ma由构造过程主控的脱气速率变化可解释δ¹³C_org曲线中的大部分长期趋势。而在早 Hettangian 期、Sinemurian-Pliensbachian 转换期及晚 Pliensbachian 期,额外碳释放源的贡献尤为关键。

(3)首次在受控于深部地质活动的脱气速率变化以及火山CO₂的释放速率序列中识别出显著的约6 Ma周期,其可能来源于天文轨道和地球深部过程的协同作用。可能存在约6 Ma尺度上的天文-构造-气候的多圈层相互作用,调控了表生环境和碳循环波动,促进了气候的周期性变化和有机碳的周期性埋藏与亏损。

衷心感谢期刊编辑团队在稿件处理过程中展现出的高效与专业。同时,我们向参与审稿的二位专家致以诚挚的谢意,他们富有洞察力的意见和建议对提升本论文的质量起到了至关重要的作用。

参考文献

[1] JIN Z, WANG X, WANG H, et al. Organic carbon cycling and black shale deposition: an Earth System Science perspective[J]. National Science Review, 2023, 10(11).

[2] MÜLLER R D, DUTKIEWICZ A. Oceanic crustal carbon cycle drives 26-million-year atmospheric carbon dioxide periodicities[J]. Science Advances, 2018, 4: eaaq0500.

[3] 刘勇胜, 陈春飞, 何德涛, 等. 古亚洲洋碳酸盐俯冲再循环及其对华北克拉通岩石圈组成的影响[J]. 中国科学:地球科学, 2021, 51(10): 1753-1772.

[4] 赵亮, 郭正堂, 袁怀玉, 等. 深部碳过程动力学模拟: 最新进展与工作流程设计[J]. 中国科学:地球科学, 2023, 53(3): 444-460.

[5] BLACK B A, KARLSTROM L, MILLS B J W, et al. Cryptic degassing and protracted greenhouse climates after flood basalt events[J]. Nature Geoscience, 2024, 17(11): 1162-1168.

[6] 宗克清, 何德涛, 陈春飞, 等. 深部碳循环的环境气候效应[J]. 岩石学报, 2022, 38(5): 1389-1398.

[7] BRUNE S, WILLIAMS S E, MÜLLER R D. Potential links between continental rifting, CO₂ degassing and climate change through time[J]. Nature Geoscience, 2017, 10(12): 941-946.

[8] GERNON T M, HINCKS T K, MERDITH A S, et al. Global chemical weathering dominated by continental arcs since the mid-Palaeozoic[J]. Nature Geoscience, 2021, 14(9): 690-696.

[9] HÜLSE D, ARNDT S, WILSON J D, et al. Understanding the causes and consequences of past marine carbon cycling variability through models[J]. Earth-Science Reviews, 2017, 171: 349-382.

[10] KENT D V, OLSEN P E, WANG H, et al. Correlation of sub-centennial-scale pulses of initial Central Atlantic Magmatic Province lavas and the end-Triassic extinctions[J]. Proceedings of the National Academy of Sciences of the United States of America, 2024, 121(46): e2415486121.

[11] HEIMDAL T H, JONES M T, SVENSEN H H. Thermogenic carbon release from the Central Atlantic magmatic province caused major end-Triassic carbon cycle perturbations[J]. Proceedings of the National Academy of Sciences of the United States of America, 2020, 117(22): 11968-11974.

[12] BOND A D, DICKSON A J, RUHL M, et al. Globally limited but severe shallow-shelf euxinia during the end-Triassic extinction[J]. Nature Geoscience, 2023, 16(12): 1181-1187.

[13] BOS R, ZHENG W, LINDSTROM S, et al. Climate-forced Hg-remobilization associated with fern mutagenesis in the aftermath of the end-Triassic extinction[J]. Nature Communications, 2024, 15(1): 3596.

[14] HESSELBO S P, GRÖCKE D R, JENKYNS H C, et al. Massive dissociation of gas hydrate during a Jurassic oceanic anoxic event[J]. Nature, 2000, 406: 392-395.

[15] HESSELBO S P, JENKYNS H C, DUARTE L V, et al. Carbon-isotope record of the Early Jurassic (Toarcian) Oceanic Anoxic Event from fossil wood and marine carbonate (Lusitanian Basin, Portugal)[J]. Earth and Planetary Science Letters, 2007, 253(3/4): 455-470.

[16] HESSELBO S P, PIEŃKOWSKI G. Stepwise atmospheric carbon-isotope excursion during the Toarcian Oceanic Anoxic Event (Early Jurassic, Polish Basin)[J]. Earth and Planetary Science Letters, 2011, 301(1/2): 365-372.

[17] REMÍREZ M N, GILLEAUDEAU G J, GAN T, et al. Carbonate uranium isotopes record global expansion of marine anoxia during the Toarcian Oceanic Anoxic Event[J]. Proceedings of the National Academy of Sciences of the United States of America, 2024, 121(27): e2406032121.

[18] STORM M S, HESSELBO S P, JENKYNS H C, et al. Orbital pacing and secular evolution of the Early Jurassic carbon cycle[J]. Proceedings of the National Academy of Sciences of the United States of America, 2020, 117(8): 3974-3982.

[19] LENTON T M, DAINES S J, MILLS B J W. COPSE reloaded: an improved model of biogeochemical cycling over Phanerozoic time[J]. Earth-Science Reviews, 2018, 178: 1-28.

[20] BERNER R A. GEOCARB-II: a revised model of atmospheric CO₂ over Phanerozoic time[J]. American Journal of Science, 1994, 294(1): 56-91.

[21] BERNER R A, KOTHAVALA Z. GEOCARB III: a revised model of atmospheric CO₂ over Phanerozoic time[J]. American Journal of Science, 2001, 301(2): 182-204.

[22] BERNER R A. Geocarbsulf: a combined model for Phanerozoic atmospheric O₂ and CO₂[J]. Geochimica et Cosmochimica Acta, 2006, 70(23): 5653-5664.

[23] ZEEBE R E. LOSCAR: long-term ocean-atmosphere-sediment carbon cycle reservoir model v2.0.4[J]. Geoscientific Model Development, 2012, 5(1): 149-166.

[24] ZHAO M, MILLS B J W, HOMOKY W B, et al. Oxygenation of the Earth aided by mineral-organic carbon preservation[J]. Nature Geoscience, 2023, 16: 262-267.

[25] LI X, HU Y, GUO J, et al. A high-resolution climate simulation dataset for the past 540 million years[J]. Scientific Data, 2022, 9(1): 371.

[26] QIU R, YU Z, MILLS B J W, et al. Pulsed biogenic methane emissions coupled with episodic warming during the Toarcian Oceanic Anoxic Event[J]. Proceedings of the National Academy of Sciences of the United States of America, 2025, 122(36): e2423598122.

[27] DAL CORSO J, MILLS B J W, CHU D, et al. Permo-Triassic boundary carbon and mercury cycling linked to terrestrial ecosystem collapse[J]. Nature Communications, 2020, 11(1): 2962.

[28] BERGMAN N M, LETON T M, WATSON A J. COPSE: a new model of biogeochemical cycling over Phanerozoic time[J]. American Journal of Science, 2000, 304: 397-437.

[29] BERNER R A, PETSCH S T, LAKE J A, et al. Isotope fractionation and atmospheric oxygen: implications for Phanerozoic O₂ evolution[J]. Science, 2000, 287(5458): 1630-1633.

[30] PETI L, THIBAULT N, CLÉMENCE M E, et al. Sinemurian-pliensbachian calcareous nannofossil biostratigraphy and organic carbon isotope stratigraphy in the Paris Basin: calibration to the ammonite biozonation of NW Europe[J]. Palaeogeography, Palaeoclimatology, Palaeoecology, 2017, 468: 142-161.

[31] XU W, RUHL M, HESSELBO S P, et al. Orbital pacing of the Early Jurassic carbon cycle, black-shale formation and seabed methane seepage[J]. Sedimentology, 2016, 64(1): 127-149.

[32] SCHUBERT B A, JAHREN A H. Global increase in plant carbon isotope fractionation following the Last Glacial Maximum caused by increase in atmospheric pCO₂[J]. Geology, 2015, 43(5): 435-438.

[33] WITKOWSKI C R, WEIJERS J W H, BLAIS B, et al. Molecular fossils from phytoplankton reveal secular pCO₂ trend over the Phanerozoic[J]. Science Advances, 2018, 4: eaat4556.

[34] ROSALES I, QUESADA S, ROBLES S. Geochemical arguments for identifying second-order sea-level changes in hemipelagic carbonate ramp deposits[J]. Terra Nova, 2006, 18(4): 233-240.

[35] SILVA R L, DUARTE L V. Organic matter production and preservation in the Lusitanian Basin (Portugal) and Pliensbachian climatic hot snaps[J]. Global and Planetary Change, 2015, 131: 24-34.

[36] DE LENA L F, TAYLOR D, GUEX J, et al. The driving mechanisms of the carbon cycle perturbations in the late Pliensbachian (Early Jurassic)[J]. Scientific Reports, 2019, 9(1): 18430.

[37] KORTE C, HESSELBO S P, ULLMANN C V, et al. Jurassic climate mode governed by ocean gateway[J]. Nature Communications, 2015, 6: 10015.

[38] KORTE C, HESSELBO S P. Shallow marine carbon and oxygen isotope and elemental records indicate icehouse-greenhouse cycles during the Early Jurassic[J]. Paleoceanography, 2011, 26(4): PA4219.

[39] FENDLEY I M, FRIELING J, MATHER T A, et al. Early Jurassic large igneous province carbon emissions constrained by sedimentary mercury[J]. Nature Geoscience, 2024, 17: 241-248.

[40] MCELWAIN J C, WADE-MURPHY J, HESSELBO S P. Changes in carbon dioxide during an oceanic anoxic event linked to intrusion into Gondwana coals[J]. Nature, 2005, 435(7041): 479-482.

[41] HEIMDAL T H, GODDERIS Y, JONES M T, et al. Assessing the importance of thermogenic degassing from the Karoo Large Igneous Province (LIP) in driving Toarcian carbon cycle perturbations[J]. Nature Communications, 2021, 12(1): 6221.

[42] GRÖCKE D R, RIMMER S M, YOKSOULIAN L E, et al. No evidence for thermogenic methane release in coal from the Karoo-Ferrar large igneous province[J]. Earth and Planetary Science Letters, 2009, 277(1/2): 204-212.

[43] SPROSON A D. Pacing of the latest Ordovician and Silurian carbon cycle by a similar to 4.5 Ma orbital cycle[J]. Palaeogeography Palaeoclimatology Palaeoecology, 2020, 540: 109543.

[44] LASKAR J, ROBUTEL P, JOUTEL F, et al. A long-term numerical solution for the insolation quantities of the Earth[J]. Astronomy & Astrophysics, 2004, 428(1): 261-85.

[45] ZHANG R, JIN Z, LI M, et al. Long-term periodicity of sedimentary basins in response to astronomical forcing: review and perspective[J]. Earth-Science Reviews, 2023, 244: 104533.

[46] BOULILA S, HAQ B U, HARA N, et al. Potential encoding of coupling between Milankovitch forcing and Earth’s interior processes in the Phanerozoic eustatic sea-level record[J]. Earth-Science Reviews, 2021, 220: 103727.

[47] BOULILA S, LASKAR J, HAQ B U, et al. Long-term cyclicities in Phanerozoic sea-level sedimentary record and their potential drivers[J]. Global and Planetary Change, 2018, 165: 128-136.

[48] BOULILA S, PETERS S E, MÜLLER R D, et al. Earth’s interior dynamics drive marine fossil diversity cycles of tens of millions of years[J]. Proceedings of the National Academy of Sciences of the United States of America, 2023, 120(29): e2221149120.

[49] BOULILA S, GALBRUN B, LASKAR J, et al. A -9 Myr cycle in Cenozoic δ¹³C record and long-term orbital eccentricity modulation: is there a link?[J]. Earth Planetary Science Letters, 2012, 317: 273-281.

[50] RAMPINO M R, CALDEIRA K. Correlation and cyclicity of stratigraphic sequence boundaries and chronostratigraphic stage boundaries of the last 253 My: synchrony of tectonism, sea level, climate and biotic change[J]. Earth-Science Reviews, 2025, 265: 105100.

[51] RAMPINO M R, CALDEIRA K, RODRIGUEZ S. Cycles of -32.5 My and -26.2 My in correlated episodes of continental flood basalts (CFBs), hyper-thermal climate pulses, anoxic oceans, and mass extinctions over the last 260 My: connections between geological and astronomical cycles[J]. Earth-Science Reviews, 2023, 246: 104548.

[52] MARTINEZ M, DERA G. Orbital pacing of carbon fluxes by a -9-My eccentricity cycle during the Mesozoic[J]. Proceedings of the National Academy of Sciences of the United States of America, 2015, 112(41): 12604-12609.

[53] IKEDA M, TADA R. A 70 million year astronomical time scale for the deep-sea bedded chert sequence (Inuyama, Japan): implications for Triassic-Jurassic geochronology[J]. Earth and Planetary Science Letters, 2014, 399: 30-43.

基金资助

地球多圈层相互作用的油气富集理论项目(THEMSIE04010101)

国家自然科学基金项目(42202162)

中国科学院地质与地球物理研究所重点部署项目(IGGCAS-202201)