巖石圈地幔分層性對克拉通穩(wěn)定性的影響

邱啟佳, 陳林, 張智, 向宵,3, 徐濤, 陳赟,白志明, 梁曉峰, 武澄瀧

1 桂林理工大學(xué)地球科學(xué)學(xué)院, 桂林 541006

2 中國科學(xué)院地質(zhì)與地球物理研究所巖石圈演化國家重點實驗室, 北京 100029

3 中國科學(xué)院大學(xué)地球與行星科學(xué)學(xué)院, 北京 100049

4 中國科學(xué)院地質(zhì)與地球物理研究所礦產(chǎn)資源研究院重點實驗室, 北京 100029

0 引言

克拉通擁有厚度大、密度輕和強(qiáng)度高的巖石圈龍骨,是大陸巖石圈中能夠長期穩(wěn)定存在的古老構(gòu)造單元(Lee et al.,2011;Pearson et al.,2021).大部分克拉通形成于太古代(>2.5 Ga;Pearson and Wittig,2008),具有長期穩(wěn)定性,但仍有些克拉通的巖石圈發(fā)生了明顯減薄,例如華北克拉通和北美Wyoming克拉通(吳福元等,2008;朱日祥等,2011;Dave and Li,2016;圖1a).全球克拉通具有明顯的多樣性特征：(1)不同克拉通之間表面地形高程的差異可達(dá)2 km(Hu et al.,2018);(2)地殼厚度的差異可達(dá)13 km(Laske et al.,2013);(3)巖石圈-軟流圈邊界(LAB)深度的差異可達(dá)200 km(Schaeffer and Lebedev,2013;Afonso et al., 2019;圖1a);(4)地幔根熔融虧損程度的深度分布差異可達(dá)～80 km(Perchuk et al.,2020);(5)巖石圈地幔底部橄欖巖交代層的深度分布差異可達(dá)～50 km(Perchuk et al.,2020).是什么因素決定了克拉通的穩(wěn)定性和多樣性?對此前人提出了多種機(jī)制,如巖石圈拆沉作用(鄧晉福等,1994)、巖石圈地幔水化作用(Kusky et al.,2014)、地幔熱柱侵蝕作用(Gorczyk et al.,2018;向宵等,2023)以及大洋板塊俯沖脫水-交代作用(朱日祥等,2011,2012)等.這些機(jī)制都要求特定的構(gòu)造環(huán)境,然而更為普遍的是古老的厚克拉通通常被年輕的薄活動帶所包圍,兩者之間存在巖石圈厚度差異(圖1a).一方面,活動帶的存在可以通過自身較弱的特性來吸收應(yīng)力,從而對克拉通起保護(hù)作用(Lenardic et al.,2000;Yoshida,2012);另一方面,活動帶與克拉通之間的巖石圈厚度突變會引起邊界驅(qū)動對流(edge-driven convection;King and Anderson,1995,1998),造成板內(nèi)火山活動(Kim and So,2020;Sun and Liu,2023)和克拉通巖石圈減薄(Hieronymus et al.,2007;van Wijk et al.,2010;Hardebol et al.,2012;劉丹紅和陳林,2019),從而降低克拉通的穩(wěn)定性.活動帶對克拉通的作用更多是保護(hù)還是破壞現(xiàn)在仍不清楚.

圖1 全球LAB深度分布及克拉通巖石圈的概念模型(a) 紅色曲線表示太古宙陸核,紫色曲線表示由多個太古宙陸核組成的復(fù)合克拉通,綠色曲線表示由多個復(fù)合克拉通組成的超級克拉通(LAB深度數(shù)據(jù)來自Afonso et al., 2019;克拉通位置據(jù)Pearson et al., 2021); (b) 克拉通巖石圈地幔分為虧損的上層、中間的MLD層及交代的下層(修改自Yuan and Romanowicz,2018),黑色箭頭表示克拉通-活動帶巖石圈臺階誘發(fā)的邊界驅(qū)動對流,紅色箭頭表示軟流圈流場.

針對邊界驅(qū)動對流與克拉通穩(wěn)定性的問題,前人開展了大量研究.大部分學(xué)者主要將克拉通巖石圈地幔視為一個均勻?qū)觼硌芯?發(fā)現(xiàn)較弱的克拉通巖石圈地幔容易被剝離減薄或拉伸變薄(Currie and van Wijk,2016;劉丹紅和陳林,2019),而強(qiáng)流變性質(zhì)和低密度能使克拉通巖石圈保持長期穩(wěn)定(Morency et al.,2002;Lenardic et al.,2003;Hieronymus et al.,2007;Beuchert et al.,2010;van Wijk et al.,2010;Hardebol et al.,2012;Kaislaniemi and van Hunen,2014;Sacek,2017).其中,Wang等(2014)通過對雙層克拉通巖石圈研究也得出了上述結(jié)論,但他們僅考慮了密度分層的影響.另外,研究者們考慮了克拉通巖石圈中弱低速層——中巖石圈不連續(xù)面(mid-lithosphere discontinuity, MLD;Selway et al.,2015)的影響后,發(fā)現(xiàn)克拉通中的弱MLD層與軟流圈的接觸范圍會影響克拉通巖石圈減薄發(fā)生的難易程度,且減薄作用會沿著MLD分布的方向持續(xù)下去(Wang et al.,2018;Shi et al.,2020,2021),若MLD被厚而冷的巖石圈所包圍,則克拉通巖石圈能夠保持長期穩(wěn)定性(Liu et al.,2018).這些研究大都建立在假設(shè)克拉通巖石圈地幔的整體成分和性質(zhì)完全一致的基礎(chǔ)之上,沒有充分考慮克拉通巖石圈地幔物性隨深度的變化.

越來越多的證據(jù)表明克拉通巖石圈地幔的成分和結(jié)構(gòu)在垂向上是非均一的.一方面,在地球物理研究中：聯(lián)合波形和SKS波分裂分析結(jié)合地球化學(xué)資料顯示,整個北美大陸下方巖石圈存在～150 km厚的化學(xué)虧損的巖石圈層和其下的～30-120 km厚的熱巖石圈根(Yuan and Romanowicz,2010);S波接收函數(shù)和表面波各向異性顯示,Kalahari克拉通巖石圈地幔存在～50 km厚的S波各向異性層、～65-115 km厚的虧損層和～60-130 km厚的交代層(Sodoudi et al.,2013);各向異性Ps接收函數(shù)(RF)結(jié)果表明,在北美的Wyoming和Superior克拉通下方發(fā)現(xiàn)多個MLD,且深度和幾何形狀均不相同(Ford et al.,2016);地震照明成像(SDI)結(jié)果顯示,西澳大利亞克拉通巖石圈地幔中存在多個MLD,且代表了巖石圈地幔上下強(qiáng)弱分層的邊界層(Sun et al.,2018);調(diào)和分解的接收函數(shù)分析顯示,澳大利亞的Yilgarn克拉通和北美的Superior克拉通的巖石圈至少分為兩層,在170 km的深度以上存在多個相互獨立的MLD(Chen et al.,2021);大地水準(zhǔn)面的校驗和殘余地形及重力資料顯示,大型克拉通下方的大部分過剩密度均分布在巖石圈地幔的下半部分,且克拉通巖石圈地幔的平均密度要比周圍軟流圈重0.5%～1.2%(Wang et al.,2022a,b).另一方面,在地球化學(xué)研究中：重礦物石榴石的主微量元素分析顯示,加拿大Slave地區(qū)的巖石圈地幔在140～160 km深度分為上下兩層,上方為橫向不均勻的化學(xué)虧損層,下方為性質(zhì)均勻的熔融交代層(Griffin et al.,2004);來自金伯利巖及相關(guān)巖石中的石榴石和鉻鐵礦捕虜體資料顯示,大部分克拉通巖石圈地幔具有巖石較富集的上層(>110 km),巖石虧損的中下層(～110-190 km)和熔體交代的最底層(～170-220 km)(Perchuk et al.,2020).綜上所述,大部分克拉通巖石圈地幔都存在垂向成分和結(jié)構(gòu)的分層現(xiàn)象(圖1b).

巖石圈地幔的成分和結(jié)構(gòu)分層性是如何影響克拉通的穩(wěn)定性?針對這一科學(xué)問題,我們聚焦于克拉通巖石圈地幔的垂向分層特性,采用二維熱-力學(xué)的數(shù)值模擬方法,探討在邊界驅(qū)動對流背景下,克拉通巖石圈地幔的密度、黏度和成分分層結(jié)構(gòu)對克拉通巖石圈穩(wěn)定性的影響.本文先對動力學(xué)方法和模型設(shè)計進(jìn)行介紹,然后系統(tǒng)測試在無MLD、連續(xù)MLD和間隔MLD條件下的克拉通巖石圈地幔密度-強(qiáng)度的分層結(jié)構(gòu),探討在邊界驅(qū)動對流作用下影響克拉通巖石圈長期穩(wěn)定性的因素,最后討論華北克拉通破壞時空范圍的成因機(jī)制.

1 數(shù)值模擬方法

本文利用二維熱-力學(xué)程序(I2VIS;Gerya and Yuen,2003;Perchuk et al.,2020),同時考慮了熱擾動來模擬邊界驅(qū)動對流(劉丹紅和陳林,2019),研究邊界驅(qū)動對流對克拉通巖石圈減薄的動力學(xué)過程.該程序的重點是粒子標(biāo)記(marker-in-cell)技術(shù)和交錯網(wǎng)格有限差分方法,聯(lián)合求解動量、物質(zhì)和能量守恒方程,不僅能考慮到局部和大尺度變形,還可以避免一定程度的數(shù)值發(fā)散.

1.1 控制方程

動量守恒方程：

(1)

(2)

式中,σ′xx、σ′yy和σ′xy為偏應(yīng)力張量;x為水平坐標(biāo),y為垂直坐標(biāo);P為壓力,T為溫度,C為巖石成分,M為部分熔融百分?jǐn)?shù),上述參數(shù)綜合控制著密度ρ;g表示重力加速度.

物質(zhì)守恒方程：

(3)

式中,vx為速度水平分量,vy為速度垂向分量.

能量守恒方程：

(4)

(5)

(6)

式中,Cp表示等壓熱容,文中取值為1000 J·kg-1·K-1;DT/Dt是溫度關(guān)于時間的物質(zhì)導(dǎo)數(shù);qx為熱流水平分量,qy為熱流垂直分量;熱傳導(dǎo)率k受溫度(T)、壓力(P)和巖石成分(C)的綜合影響(Clauser and Huenges,1995);α為熱膨脹系數(shù);Hr、Hs、Ha和HL分別代表放射性、剪切、絕熱和潛熱生熱.

1.2 黏-塑性流變

巖石的流變性質(zhì)受塑性和黏性流變機(jī)制的共同控制,偏應(yīng)力張量和應(yīng)變率張量的關(guān)系如下：

(7)

(8)

塑性變形遵從Drucker-Prager屈服準(zhǔn)則(Ranalli,1995)：

(9)

σyield=C0+Psin(φdry)(1-λ),

(10)

(11)

韌性變形綜合考慮了擴(kuò)散蠕變和位錯蠕變：

(12)

式中,巖石成分、溫度和壓力等共同控制著擴(kuò)散蠕變黏度ηdiffusion和位錯蠕變黏度ηdislocation(Ranalli,1995)：

(13)

(14)

式中,σcr表示擴(kuò)散蠕變與位錯蠕變發(fā)生轉(zhuǎn)變的臨界應(yīng)力;R表示氣體常數(shù);物質(zhì)常數(shù)(AD)、活化能(Ea)、活化體積(Va)和應(yīng)力指數(shù)(n)是巖石的流變參數(shù).

Peierls變形出現(xiàn)于高壓低溫區(qū)域(Karato et al.,2001)：

(15)

(16)

式中,σpeierls為限制物質(zhì)強(qiáng)度的Peierls應(yīng)力,與式(9)中的屈服應(yīng)力σyield相似;σII表示應(yīng)力張量二階不變量;P、T、R、Ea、Va和n意義與式(14)相同,Apeierls為物質(zhì)常數(shù).

巖石圈的有效黏度ηeff為塑性黏度ηplastic、韌性黏度ηductile和Peierls黏度ηpeierls之間的最小值：

ηeff=min{ηplastic,ηductile,ηpeierls}.

(17)

1.3 部分熔融

巖石的部分熔融百分?jǐn)?shù)M與溫度為簡單的線性關(guān)系(Gerya and Yuen,2003;Burg and Gerya,2005)：

(18)

式中,Ts為巖石固相線的溫度,Tl為巖石液相線的溫度(計算表達(dá)式見表1).

部分熔融巖石的有效密度ρeff：

ρeff=ρs-M(ρs-ρm),

(19)

式中,ρs為巖石固相的密度,ρm為巖石熔融的密度(計算表達(dá)式見表1).

1.4 長波長熱擾動

為了模擬邊界驅(qū)動對流,在軟流圈中加入長波長熱擾動(劉丹紅和陳林,2019)：

(20)

式中,Tp為熱擾動幅度,ylab和Tlab分別表示巖石圈底界面的深度和溫度,y0和x0分別表示模型的垂向深度和水平寬度,y和x分別表示模型中的垂向位置和水平位置.

長波長熱擾動的引入可以在巖石圈地幔下方形成橫向不均一的溫度結(jié)構(gòu),這種溫度異?？梢孕纬纱蟪叨鹊蒯α?在不施加邊界驅(qū)動力的情況,增加必要的熱擾動可以有效地驅(qū)動模型演化(King and Anderson,1998).熱擾動引發(fā)的大尺度地幔對流與克拉通-活動帶巖石圈厚度橫向急劇變化引發(fā)的小尺度邊界驅(qū)動對流相互作用,可以起到抑制小尺度邊界驅(qū)動對流過快發(fā)展的作用(劉丹紅和陳林,2019).經(jīng)過反復(fù)測試后,我們發(fā)現(xiàn)：邊界驅(qū)動對流在Tp小于Tlab的1%時完全占主導(dǎo),大尺度地幔對流在Tp大于Tlab的8%時完全占主導(dǎo).因此,選擇Tp=40 ℃(約Tlab=1350 ℃的3%)來設(shè)計軟流圈溫度場,既能模擬大尺度地幔對流,又有利于形成邊界驅(qū)動對流.

2 模型設(shè)計

基于克拉通的分層結(jié)構(gòu)以及克拉通-活動帶的巖石圈厚度差異,我們設(shè)計了一個初始模型(圖2).模型具有4000 km的水平寬度和1000 km的垂向深度(圖2a),通過由2001×381個節(jié)點組成的非均勻網(wǎng)格系統(tǒng)來劃分模型空間.在水平方向上,模型是2 km的平均網(wǎng)格分辨率.在垂直方向上,模型0～200 km深度的分辨率為1 km,200～300 km深度的分辨率從1 km逐漸降至5 km,300～1000 km的分辨率為5 km.如此,既可以保證克拉通巖石圈地幔處的分辨率,又可以提高計算效率.

圖2 初始模型(a)模型設(shè)計;(b)為(a)中放大的模型區(qū)域,白色實線為等溫線,黑色和紅色實線分別表示0.5 Myr時刻巖石圈垂向剖面的密度和黏度.組分色標(biāo)：0,黏性空氣;5和6,大陸上下地殼;9,克拉通上層地幔巖石圈;10,軟流圈;11,MLD;12,活動帶地幔巖石圈;14,克拉通下層地幔巖石圈.

模型左側(cè)和右側(cè)的塊體分別代表克拉通和活動帶,兩者的地殼性質(zhì)一樣,由均為20 km的上、下地殼組成.地殼下覆由40 km的上層巖石圈地幔、20 km的MLD和60 km的下層巖石圈地幔組成,克拉通巖石圈地?？偣埠?20 km,活動巖石圈地幔性質(zhì)均一,厚度為40 km,克拉通-活動帶厚度差異為80 km(圖2b).我們采用巖石力學(xué)實驗確定的不同礦物的流動律(flow law)來描述模型不同成分層的流變學(xué)性質(zhì)(如, Chen, 2021),以“濕石英”和“斜長石”分別代表上、下地殼的流變學(xué)性質(zhì),以“濕橄欖石”代表MLD的流變學(xué)性質(zhì),以“干橄欖石”代表克拉通和活動帶地幔及軟流圈的流變學(xué)性質(zhì)(具體參數(shù)見表1).克拉通的上地殼和下地殼的地溫梯度分別為15 ℃/km和10 ℃/km,活動帶的上地殼和下地殼的地溫梯度分別為20 ℃/km和12.5 ℃/km,克拉通和活動帶的巖石圈底界面溫度Tlab均為1350 ℃,軟流圈中的絕熱地溫梯度K=0.5 ℃/ km.

模型的上邊界和左右兩側(cè)邊界為自由滑邊界條件,下邊界為滲透性邊界條件(Burg and Gerya,2005),計算區(qū)域下方500 km處將符合自由滑條件,模型底部(y=1000 km)采用的滲透性邊界條件,可以解決模型底部的部分邊界效應(yīng)問題.模型上邊界溫度為常溫(0 ℃),左右兩側(cè)邊界為絕熱邊界,下邊界為外部溫度固定的邊界條件(Burg and Gerya,2005),計算區(qū)域下方500 km處將符合固定溫度條件,該區(qū)域的熱流和溫度能夠隨著演化而自行調(diào)整.模型演化全程無外力擾動,全由初始狀態(tài)驅(qū)動.模型頂部設(shè)置有20 km厚的偽空氣層,其特點為低密度(1 kg·m-3)和低黏度(1019Pa·s),可模擬減薄過程中的地形演化(Schmeling et al.,2008).

3 模擬結(jié)果

在數(shù)值模擬過程中,我們將活動帶和軟流圈的參考密度和參考黏度保持不變.在改變克拉通巖石圈地幔的密度結(jié)構(gòu)時,通過巖石圈地幔各層厚度來加權(quán)調(diào)節(jié)克拉通上、下層巖石圈地幔的密度,保持整個克拉通巖石圈地幔的參考密度與軟流圈的參考密度基本一致.其中,厚度加權(quán)計算公式為：

ρuclmhuclm+ρMLDhMLD+ρlclmhlclm=ρmlmhmlm+ρamham,

(21)

式中,ρuclm、ρMLD、ρlclm、ρmlm和ρam分別為克拉通上層巖石圈地幔、MLD層、克拉通下層巖石圈地幔、活動帶巖石圈地幔和軟流圈地幔;huclm、hMLD、hlclm、hmlm和ham分別為克拉通上層巖石圈地幔、MLD層、克拉通下層巖石圈地幔、活動帶巖石圈地幔和軟流圈地幔的厚度.

本文系統(tǒng)測試了在邊界驅(qū)動對流的作用下,克拉通巖石圈地幔的垂向密度-黏度分層結(jié)構(gòu)及MLD的分布樣式對克拉通巖石圈穩(wěn)定性的影響.通過對不同的Δρ0、Δη0及MLD樣式A和B進(jìn)行組合,我們以無MLD情況下的Δρ0=40 kg·m-3和Δη0=50的模型Model-11為實驗的參考模型,在此基礎(chǔ)上開展了42組數(shù)值模擬實驗(表2),進(jìn)而研究在邊界驅(qū)動對流的條件下,克拉通巖石圈地幔的不同分層結(jié)構(gòu)對其穩(wěn)定性的影響.

3.1 參考模型(無MLD)

圖3為參考模型Model-11的演化結(jié)果.長波長熱擾動驅(qū)使軟流圈呈現(xiàn)自西向東的流動,流經(jīng)克拉通和活動帶的交界處.由于克拉通和活動帶巖石圈厚度的突變,模型在演化到10 Myr時克拉通-活動帶巖石圈臺階處出現(xiàn)了明顯的邊界驅(qū)動對流(圖3a1,b1).較重的克拉通下層巖石圈地幔受邊界驅(qū)動對流的影響出現(xiàn)對流侵蝕,被拆沉墜入軟流圈并滯留于地幔過渡帶(圖3a2,b2).當(dāng)模型演化至160 Myr,由于邊界驅(qū)動對流的侵蝕,較重的克拉通下層巖石圈地幔自東向西橫向減薄了～200 km,而較輕的上層巖石圈地幔只是輕微變形,巖石圈-活動帶間的過渡帶被剝蝕而變得平緩,地幔對流也隨之減慢(圖3a3,b3,c).無MLD的克拉通巖石圈地幔分層模型結(jié)果顯示：由于克拉通下層巖石圈地幔比軟流圈更重更強(qiáng),克拉通只出現(xiàn)了邊緣減薄現(xiàn)象,其整體巖石圈擁有較長的壽命,具有長期的穩(wěn)定性,總體上與劉丹紅和陳林(2019)呈現(xiàn)的結(jié)果一致.

3.2 MLD連續(xù)分布對克拉通穩(wěn)定性的影響

3.2.1 克拉通巖石圈密度分層結(jié)構(gòu)的影響

這里不改變克拉通下層巖石圈的參考黏度,系統(tǒng)測試克拉通下層巖石圈的參考密度對克拉通穩(wěn)定性的影響.圖4顯示了模型Model-23(Δρ0=40 kg·m-3,Δη0=50)和模型Model-22(Δρ0=20 kg·m-3,Δη0=50)的對比結(jié)果,兩個模型的克拉通下層巖石圈參考黏度均為軟流圈參考黏度的50倍,僅在克拉通上、下層巖石圈地幔的參考密度上存在差異(圖4d).模型Model-23顯示：邊界驅(qū)動對流的下降流會拆掉克拉通下層巖石圈地幔并使其墜入地幔過渡帶,而上升流所附帶的熱物質(zhì)會使MLD發(fā)生部分熔融,克拉通下層巖石圈的拆沉?xí)诳死▋?nèi)部發(fā)育出新的巖石圈臺階,從而使對流減薄自東向西持續(xù),在80 Myr時產(chǎn)生～1450 km的橫向減薄范圍(圖4a1,a2,c).模型Model-22顯示：當(dāng)克拉通下層巖石圈的參考密度減輕20 kg·m-3時,邊界驅(qū)動對流對克拉通的減薄影響減弱,減薄規(guī)模也相對較小,80 Myr的克拉通橫向減薄范圍僅有～150 km(圖4b1,b2,c).在克拉通內(nèi)MLD連續(xù)分布的情況下,克拉通巖石圈地幔內(nèi)部密度分層的差異對比表明：當(dāng)克拉通下層巖石圈更重時,其受到的邊界驅(qū)動對流影響更大.

表2 數(shù)值模擬實驗Table 2 Numerical modelling experiments

3.2.2 克拉通巖石圈強(qiáng)度分層結(jié)構(gòu)的影響

為了研究克拉通巖石圈地幔強(qiáng)度分層結(jié)構(gòu)的影響,克拉通巖石圈地幔的參考密度將保持不變,測試克拉通下層巖石圈地幔的不同參考黏度對克拉通穩(wěn)定性的影響.圖5顯示了模型Model-23(Δρ0=40 kg·m-3,Δη0=50)和模型Model-15(Δρ0=40 kg·m-3,Δη0=1)的對比結(jié)果,兩個模型的克拉通巖石圈地幔都保持相同的密度分層結(jié)構(gòu),僅在克拉通下層巖石圈地幔的黏度上存在差異(圖5d).當(dāng)克拉通下層巖石圈與軟流圈的參考黏度差Δη0=50(Model-23)時,強(qiáng)的克拉通下層巖石圈地幔容易受邊界驅(qū)動對流侵蝕的影響,產(chǎn)生大規(guī)模的塊狀拆沉,75 Myr時的減薄范圍為～1000 km(圖5a1,a2).而當(dāng)克拉通下層巖石圈與軟流圈無參考黏度差(Model-15)時,在75 Myr時克拉通巖石圈自東向西的橫向減薄范圍僅為～300 km,但持續(xù)的對流會使致密下層巖石圈地幔出現(xiàn)重力失穩(wěn)現(xiàn)象,從而在巖石圈底界面繼續(xù)產(chǎn)生向西～900 km的10～15 km厚的滴落狀剝離(圖5b1,b2).在克拉通內(nèi)MLD連續(xù)分布的情況下,克拉通巖石圈地幔內(nèi)部黏度分層的差異對比表明：當(dāng)克拉通下層巖石圈更強(qiáng)時,其受到的邊界驅(qū)動對流影響更大.

3.3 MLD不連續(xù)分布對克拉通穩(wěn)定性的影響

大量的地震觀測表明(Ford et al.,2016;Sun et al.,2018;Chen et al.,2021),MLD可能呈橫向不連續(xù)分布,因此,我們在克拉通巖石圈地幔中設(shè)置了不同間隔距離的不同寬度的多個MLD,以測試MLD的不連續(xù)分布對克拉通穩(wěn)定性的影響.

圖3 參考模型的成分場和溫度場演化白色箭頭和實線分別表示地幔流場的方向和大小.

圖4 克拉通巖石圈密度分層結(jié)構(gòu)對克拉通穩(wěn)定性的影響(Δη0=50)(a1)和(a2)分別為模型Model-23(Δρ0=40 kg·m-3)演化到80 Myr的成分場和密度場;(b1)和(b2)分別為模型Model-22(Δρ0=20 kg·m-3)演化到80 Myr的成分場和密度場;(c)為模型Model-23和Model-22的LAB深度曲線;(d)表示模型Model-23和Model-22在0.4 Myr時X=1800 km處的克拉通巖石圈密度剖面.

圖5 克拉通巖石圈強(qiáng)度分層結(jié)構(gòu)對克拉通穩(wěn)定性的影響(Δρ0=40 kg·m-3)(a1)和(a2)分別為模型Model-23(Δη0=50)演化到75 Myr的成分場和黏度場;(b1)和(b2)分別為模型Model-15(Δη0=1)演化到75 Myr的成分場和黏度場;(c)為模型Model-23和Model-15的LAB深度曲線;(d)表示模型Model-23和Model-15在0.4 Myr時X=1800 km處的克拉通巖石圈黏度剖面.

圖6 MLD的分布樣式對克拉通穩(wěn)定性的影響(Δρ0=40 kg·m-3,Δη0=50,B=10 km)(a1)—(a3)和(b1)—(b3)分別為模型Model-25(A=50 km)和模型Model-27(A=200 km)在5、50和100 Myr的成分場演化;(a4)和(b4)分別為模型Model-25和Model-27的LAB深度曲線.

圖7 克拉通巖石圈密度分層結(jié)構(gòu)對克拉通穩(wěn)定性的影響(A=100 km,B=10 km,Δη0=50)(a1)—(a3)和(b1)—(b3)分別為模型Model-32(Δρ0=20 kg·m-3)和模型Model-26(Δρ0=40 kg·m-3)在10、60、120 Myr的成分場演化;(a4)和(b4)分別為模型Model-32和模型Model-26的LAB深度曲線.

圖8 克拉通巖石圈強(qiáng)度分層結(jié)構(gòu)對克拉通穩(wěn)定性的影響(A=100 km,B=10 km,Δρ0=20 kg·m-3)(a1)—(a3)和(b1)—(b3)分別為模型Model-38(Δη0=20)和Model-32(Δη0=50)在50、100和160 Myr的成分場演化;(a4)和(b4)分別為模型Model-38和模型Model-32的LAB深度曲線.

圖9 巖石圈的不同分層結(jié)構(gòu)對克拉通穩(wěn)定性的影響模式圖不同顏色的圓形代表克拉通巖石圈減薄的程度：綠色表示輕微減薄,藍(lán)色表示橫向減薄范圍不超過200 km,紅色表示橫向減薄范圍超過200 km;不同線段組合的白色圓形代表克拉通巖石圈減薄的減薄方式.

3.3.1 MLD分布樣式的影響

保持克拉通巖石圈地幔中的密度(Δρ0=40 kg·m-3)和黏度(Δη0=50)分層結(jié)構(gòu)不變,以此測試不同間距的MLD分布對克拉通穩(wěn)定性的影響,模型結(jié)果顯示：MLD分布較為密集的模型Model-25(A=50 km,B=10 km)的橫向減薄規(guī)模小,且減薄速度較為緩慢(圖6a1,a2),在100 Myr時的克拉通巖石圈橫向減薄范圍僅為～200 km(圖6a3),而MLD分布較為稀疏的模型Model-27(A=200 km,B=10 km)的橫向減薄規(guī)模較大,且減薄速度也較快(圖6b1,b2),在100 Myr時的克拉通橫向減薄范圍可達(dá)～850 km(圖6b3).在不連續(xù)MLD分布情況下,更稀疏的MLD分布更有利于邊界驅(qū)動對流對克拉通的減薄破壞(圖6a4,b4).

3.3.2 克拉通巖石圈密度分層結(jié)構(gòu)的影響

以MLD分布較為密集的模型(A=100 km,B=10 km)為基礎(chǔ),保持克拉通下層巖石圈地幔的參考黏度(Δη0=50)不變,測試克拉通巖石圈地幔的不同密度分層結(jié)構(gòu)對克拉通穩(wěn)定性的影響.模型Model-32(Δρ0=20 kg·m-3)同Model-26(Δρ0=40 kg·m-3)比較顯示,更輕的克拉通下層巖石圈的橫向減薄規(guī)模小,且減薄速度慢(圖7a1—a2,b1—b2),在120 Myr時,Model-32的克拉通巖石圈橫向減薄范圍比Model-26少1000 km(圖7a3,b3).在克拉通內(nèi)MLD分布較為密集的情況下,克拉通巖石圈地幔內(nèi)部密度分層的差異對克拉通的穩(wěn)定性影響同連續(xù)MLD分布下的結(jié)果趨勢一致：重的克拉通下層巖石圈受到的邊界驅(qū)動對流影響更大(圖7a4,b4).

3.3.3 克拉通巖石圈強(qiáng)度分層結(jié)構(gòu)的影響

繼續(xù)保持相同的MLD分布樣式(A=100 km,B=10 km),在克拉通下層巖石圈地幔的參考密度(Δρ0=20 kg·m-3)不變的情況下,測試克拉通巖石圈地幔的不同黏度分層結(jié)構(gòu)對克拉通穩(wěn)定性的影響.模型結(jié)果顯示：克拉通下層巖石圈更強(qiáng)的模型Model-32(Δη0=50)同相對較弱的模型Model-38(Δη0=20)對比,巖石圈橫向減薄規(guī)模更大,且減薄速度更快(圖8a1—a2,b1—b2),在160 Myr時Model-38的克拉通巖石圈橫向減薄范圍僅為～500 km(圖8a3),而Model-32能達(dá)～950 km(圖8b3).在克拉通內(nèi)MLD分布較為密集的情況下,克拉通巖石圈地幔內(nèi)部黏度分層的差異對克拉通的穩(wěn)定性影響也同連續(xù)MLD分布下的結(jié)果趨勢一致：更強(qiáng)的克拉通下層巖石圈更容易受到邊界驅(qū)動對流的影響而發(fā)生拆沉(圖8a4,b4).

4 討論

4.1 巖石圈地幔分層性對穩(wěn)定性的影響

克拉通巖石圈地幔分層結(jié)構(gòu)在俯沖構(gòu)造背景下的巖石圈拆沉以及拆沉后再分層的演化過程在模擬實驗中已有研究(如,Peng et al.,2022;Cao and Liu,2023),而圖9顯示了邊界驅(qū)動對流背景下的克拉通巖石圈地幔內(nèi)的垂向密度-黏度分層和MLD分布樣式對其穩(wěn)定性的影響.在不考慮MLD的情況下,模擬結(jié)果和前人研究相似：相對重的下層巖石圈地幔更容易受到外力而發(fā)生巖石圈減薄,其強(qiáng)度增加則有利于保持穩(wěn)定(劉丹紅和陳林,2019;Currie et al.,2023).而考慮MLD時,弱MLD的存在不僅可以導(dǎo)致大范圍的減薄發(fā)生,并且會使克拉通減薄從其邊緣沿著MLD分布的方向進(jìn)行(Wang et al.,2018;Liu et al.,2018).其中,只有當(dāng)克拉通下層巖石圈地幔同軟流圈保持同等密度時,更強(qiáng)的克拉通下層巖石圈地幔才更有利于保持克拉通的穩(wěn)定;而若克拉通下層巖石圈地幔更重時(Δρ0≥20 kg·m-3),邊界驅(qū)動對流會導(dǎo)致致密的克拉通下層巖石圈出現(xiàn)明顯的橫向減薄,并沿著MLD的分布方向持續(xù)進(jìn)行,此時,克拉通下層巖石圈地幔的強(qiáng)度增加則不利于克拉通的穩(wěn)定性.

不同構(gòu)造背景(俯沖、地幔柱)下的克拉通巖石圈中的MLD位置分布與外擾動產(chǎn)生的變形會影響克拉通巖石圈的穩(wěn)定性(Shi et al.,2021,2022),而且在邊界驅(qū)動對流的構(gòu)造背景下,克拉通巖石圈中MLD的存在會使其減薄破壞呈現(xiàn)一個前期慢速減薄和后期快速減薄的非均勻過程(圖10).當(dāng)克拉通下層巖石圈地幔更重(Δρ0≥20 kg·m-3)以及更強(qiáng)(Δη0≥20)時,在克拉通巖石圈減薄過程中,其減薄方式會發(fā)生轉(zhuǎn)變,具體體現(xiàn)在長期(～60-100 Myr)的慢速的小規(guī)模絲狀減薄向短期(～10-20 Myr)的快速的大規(guī)模塊狀拆沉的過渡,且克拉通巖石圈前期的橫向減薄范圍少(～200-400 km),后期的橫向減薄范圍多(～600-800 km)(如, Model-23),此時,若增強(qiáng)下層巖石圈地幔的強(qiáng)度,則會使減薄方式轉(zhuǎn)變的時間提前,即克拉通巖石圈的大規(guī)模塊狀拆沉提前發(fā)生,甚至克拉通巖石圈的減薄方式一開始就是塊狀拆沉(如, Model-24).而當(dāng)克拉通巖石圈的MLD呈不連續(xù)分布(B=10 km)時,與MLD連續(xù)分布的情況下的克拉通巖石圈地幔密度和黏度分層結(jié)構(gòu)的演化趨勢一致,會有一個減薄方式轉(zhuǎn)變的過渡,但由于MLD的橫向不連續(xù),減薄方式的轉(zhuǎn)變時間往后延長了～30 Myr,并且會降低克拉通的橫向減薄速度,從而延長克拉通的壽命(圖10),且更密集的MLD分布樣式有利于克拉通巖石圈抵抗邊界驅(qū)動對流的侵蝕(圖6).

圖10 不同巖石圈地幔分層的克拉通巖石圈減薄隨時間的水平范圍虛線和實線分別表示在無MLD(黑色)、連續(xù)MLD(藍(lán)色)和不連續(xù)MLD(紅色)情況下,克拉通下層巖石圈地幔與軟流圈密度差值Δρ0為20 kg·m-3和40 kg·m-3的模型在150 Myr時的克拉通減薄的水平范圍(黏度差值Δη0固定為50).

4.2 對華北克拉通巖石圈減薄過程的啟示

華北克拉通是研究克拉通破壞的典型實例之一(朱日祥等,2011;Wu et al.,2019).通過古生代的金伯利巖和新生代的玄武巖地幔包體約束表明,古生代(～480 Ma)華北克拉通的巖石圈厚約200 km,而新生代則不足80 km(Menzies et al.,1993;Griffin et al.,1998);地震接收函數(shù)顯示,現(xiàn)今華北克拉通西部的巖石圈尚存200 km,而東部的巖石圈僅有60～80 km(Chen et al.,2008).巖漿巖的性質(zhì)和同位素示蹤研究顯示華北克拉通的減薄可能始于晚石炭紀(jì)-晚三疊世(300～220 Ma; Yang et al.,2008;徐義剛等,2009),在早白堊世(約130～120 Ma)達(dá)到減薄峰期(姜耀輝等,2005;Hou et al.,2007),但并沒有減薄到現(xiàn)今巖石圈的厚度(徐義剛,2004,2006).這表明華北克拉通東部巖石圈的厚度自古生代以來至少減薄100 km,而西部巖石圈尚保持完整.目前,主流觀點認(rèn)為華北克拉通東部的破壞是由太平洋板塊俯沖引起的巖石圈地幔熔/流體交代作用及其伴隨的海溝后撤引發(fā)的強(qiáng)烈伸展作用共同作用的結(jié)果(吳福元等,2000;Zheng et al.,2006,2007;Xu,2007;朱日祥等,2011,2012).目前認(rèn)為華北克拉通的破壞呈現(xiàn)時空不均勻性的特點(徐義剛,2004;Zheng et al.,2006,2007),且存在不止一期的巖石圈減薄破壞,華北克拉通的減薄很可能是一個循序漸進(jìn)的漫長過程(徐義剛等,2009;劉丹紅和陳林,2019).而最新的華北克拉通鎂鐵質(zhì)巖漿的平衡壓力研究顯示華北克拉通破壞的時間集中在早白堊世(～120 Ma)期間,華北克拉通的部分東部巖石圈在極短的時間內(nèi)(10 Ma)從200 km減薄到35 km(Chen et al.,2023),這種短時間內(nèi)的大規(guī)模減薄難以用太平洋板塊俯沖模型來解釋.巖石地球化學(xué)證據(jù)表明：華北克拉通破壞前的巖石圈地幔根部的黏度是軟流圈的2～8倍(Wang et al.,2022),在華北克拉通中部的巖石圈地幔中有上部的殘留古老地幔和下部的改造地幔的雙層特征(湯艷杰等,2021);地球物理證據(jù)表明：在現(xiàn)今的華北克拉通西部陸塊和中部跨華北造山帶的巖石圈地幔中分別觀測到MLD的存在,且MLD存在的深度與華北克拉通東部陸塊的LAB深度一致(Chen et al.,2014).綜上所述,華北克拉通在破壞前存在明顯的巖石圈地幔成分和結(jié)構(gòu)的分層,并且破壞前的東部巖石圈地幔中可能存在MLD,但與中部巖石圈地幔中的MLD間隔較遠(yuǎn).

MLD不連續(xù)分布的巖石圈分層模型結(jié)果顯示(圖9),在開始的很長一段時間內(nèi)(90～130 Myr),克拉通下層巖石圈地幔只有其邊緣向西200～400 km的局部區(qū)域發(fā)生減薄,之后僅約30 Myr就能完成繼續(xù)向西約700 km的大規(guī)模減薄,這與上述的破壞時間集中的特征一致,這表明華北克拉通巖石圈地幔中很可能存在橫向MLD不連續(xù)分布的垂向密度-黏度分層結(jié)構(gòu).因此,我們提出一個可能的華北克拉通破壞過程的新模型：(1)130 Ma之前減薄開始：華北克拉通的東部巖石圈由于受到周圍的年輕大洋板塊或造山帶引起的邊界驅(qū)動對流,其邊緣的MLD剛開始受到的外力擾動影響不大,巖石圈產(chǎn)生規(guī)模較小、速度較慢且持續(xù)時間很長的絲狀減薄(圖11a);(2)130 Ma到120 Ma期間減薄進(jìn)入峰期：當(dāng)華北克拉通東部巖石圈中的橫向?qū)挾容^窄的單個MLD被完全拆掉后,持續(xù)的對流開始剝離相鄰且橫向?qū)挾容^長的MLD,巖石圈減薄由絲狀剝離轉(zhuǎn)變?yōu)閴K狀拆沉,導(dǎo)致華北克拉通巖石圈在短時間內(nèi)產(chǎn)生大規(guī)模的減薄(圖11b);(3)120 Ma至現(xiàn)今減薄減緩或結(jié)束：當(dāng)巖石圈減薄到達(dá)相鄰MLD間隔距離很大的華北克拉通中部時,對流減薄作用很難影響到華北克拉通西部內(nèi)被厚巖石圈包圍的MLD,巖石圈減薄方式又轉(zhuǎn)變?yōu)榻z狀剝離,且減薄速度減緩,直至形成現(xiàn)今位于華北克拉通中部的平緩巖石圈臺階(圖11c).由此,我們認(rèn)為克拉通巖石圈地幔的垂向密度-強(qiáng)度分層結(jié)構(gòu)和MLD的分布樣式在華北克拉通破壞過程中發(fā)揮了重要作用.

5 結(jié)論

本文采用二維熱-力學(xué)數(shù)值模擬方法探究在邊界驅(qū)動對流條件下,克拉通巖石圈地幔的垂向密度-強(qiáng)度分層結(jié)構(gòu)及其內(nèi)部MLD的分布樣式對克拉通穩(wěn)定性的影響,獲得的主要結(jié)論如下：

(1)當(dāng)克拉通巖石圈內(nèi)部不存在MLD時,克拉通巖石圈減薄的方式主要為絲狀剝離,在下層巖石圈地幔較輕(Δρ0≤40 kg·m-3)且較強(qiáng)(Δη0≥1)的情況下,其減薄水平范圍小于100 km,且主要發(fā)生在克拉通邊緣處.

(2)當(dāng)克拉通巖石圈內(nèi)部存在一層連續(xù)的MLD時,在克拉通下層巖石圈地幔較重(Δρ0>0 kg·m-3)的情況下,克拉通巖石圈減薄的水平范圍隨著克拉通下層巖石圈地幔密度的增加而增大,若克拉通下層巖石圈較弱(Δη0=1),巖石圈地幔的減薄受對流侵蝕和重力失穩(wěn)的共同影響;若克拉通下層巖石圈較強(qiáng)(Δη0≥1),巖石圈地幔的減薄僅受對流侵蝕影響,其減薄方式隨著其強(qiáng)度的增加從緩慢絲狀剝離轉(zhuǎn)變?yōu)榭焖賶K狀拆沉,甚至從開始階段就越過絲狀剝離而直接變成塊狀拆沉.

(3)當(dāng)克拉通巖石圈內(nèi)部存在一層不連續(xù)分布的MLD時,在克拉通下層巖石圈地幔的密度和黏度不變的條件下,單個MLD的寬度A=100 km且相鄰兩個MLD之間的間隙B=10 km的不連續(xù)分布樣式能夠?qū)p薄方式改變的時間延后～30 Myr,并降低減薄的速度,有效延緩巖石圈地幔的減薄.間隙越稀疏巖石圈減薄越快,反之亦然,減薄方式早期表現(xiàn)為緩慢絲狀剝離,中后期表現(xiàn)為快速塊狀拆沉.

(4)華北克拉通破壞具有持續(xù)時間長(>100 Ma)、峰期集中(130～120 Ma)的特點,并且其破壞的前緣位置與中部造山帶下方MLD的東端對應(yīng).根據(jù)模擬結(jié)果,我們認(rèn)為華北克拉通破壞的時空范圍呈現(xiàn)的不均一現(xiàn)象,可能主要是由巖石圈地幔的垂向密度-強(qiáng)度分層結(jié)構(gòu)和MLD不連續(xù)分布引發(fā)的巖石圈地幔逐段拆沉所導(dǎo)致.

圖11 華北克拉通分期破壞過程的示意圖(a) 130 Ma之前,華北克拉通巖石圈開始減薄,減薄規(guī)模小、速度慢; (b) 130 Ma到120 Ma期間,華北克拉通巖石圈進(jìn)入減薄峰期,減薄規(guī)模大、速度快; (c) 120 Ma至現(xiàn)今,華北克拉通巖石圈減薄緩慢,甚至趨近于停止.

致謝謹(jǐn)以此文祝賀滕吉文先生90華誕暨從事地球物理工作70年.感謝瑞士蘇黎世聯(lián)邦理工學(xué)院Taras Gerya教授提供的I2VIS程序.感謝審稿人對本文提出的建設(shè)性意見.感謝中國科學(xué)院地質(zhì)與地球物理研究所田小波研究員、顏智勇博士、崔曉娜博士、司潤港博士對本文的建議和幫助.本文的模擬計算是在國家超級計算天津中心和北京超級云計算中心完成.