岩石圈地幔分层性对克拉通稳定性的影响

邱启佳, 陈林, 张智, 向宵,3, 徐涛, 陈赟,白志明, 梁晓峰, 武澄泷

1 桂林理工大学地球科学学院, 桂林 541006

2 中国科学院地质与地球物理研究所岩石圈演化国家重点实验室, 北京 100029

3 中国科学院大学地球与行星科学学院, 北京 100049

4 中国科学院地质与地球物理研究所矿产资源研究院重点实验室, 北京 100029

0 引言

克拉通拥有厚度大、密度轻和强度高的岩石圈龙骨,是大陆岩石圈中能够长期稳定存在的古老构造单元(Lee et al.,2011;Pearson et al.,2021).大部分克拉通形成于太古代(>2.5 Ga;Pearson and Wittig,2008),具有长期稳定性,但仍有些克拉通的岩石圈发生了明显减薄,例如华北克拉通和北美Wyoming克拉通(吴福元等,2008;朱日祥等,2011;Dave and Li,2016;图1a).全球克拉通具有明显的多样性特征：(1)不同克拉通之间表面地形高程的差异可达2 km(Hu et al.,2018);(2)地壳厚度的差异可达13 km(Laske et al.,2013);(3)岩石圈-软流圈边界(LAB)深度的差异可达200 km(Schaeffer and Lebedev,2013;Afonso et al., 2019;图1a);(4)地幔根熔融亏损程度的深度分布差异可达～80 km(Perchuk et al.,2020);(5)岩石圈地幔底部橄榄岩交代层的深度分布差异可达～50 km(Perchuk et al.,2020).是什么因素决定了克拉通的稳定性和多样性?对此前人提出了多种机制,如岩石圈拆沉作用(邓晋福等,1994)、岩石圈地幔水化作用(Kusky et al.,2014)、地幔热柱侵蚀作用(Gorczyk et al.,2018;向宵等,2023)以及大洋板块俯冲脱水-交代作用(朱日祥等,2011,2012)等.这些机制都要求特定的构造环境,然而更为普遍的是古老的厚克拉通通常被年轻的薄活动带所包围,两者之间存在岩石圈厚度差异(图1a).一方面,活动带的存在可以通过自身较弱的特性来吸收应力,从而对克拉通起保护作用(Lenardic et al.,2000;Yoshida,2012);另一方面,活动带与克拉通之间的岩石圈厚度突变会引起边界驱动对流(edge-driven convection;King and Anderson,1995,1998),造成板内火山活动(Kim and So,2020;Sun and Liu,2023)和克拉通岩石圈减薄(Hieronymus et al.,2007;van Wijk et al.,2010;Hardebol et al.,2012;刘丹红和陈林,2019),从而降低克拉通的稳定性.活动带对克拉通的作用更多是保护还是破坏现在仍不清楚.

图1 全球LAB深度分布及克拉通岩石圈的概念模型(a) 红色曲线表示太古宙陆核,紫色曲线表示由多个太古宙陆核组成的复合克拉通,绿色曲线表示由多个复合克拉通组成的超级克拉通(LAB深度数据来自Afonso et al., 2019;克拉通位置据Pearson et al., 2021); (b) 克拉通岩石圈地幔分为亏损的上层、中间的MLD层及交代的下层(修改自Yuan and Romanowicz,2018),黑色箭头表示克拉通-活动带岩石圈台阶诱发的边界驱动对流,红色箭头表示软流圈流场.

针对边界驱动对流与克拉通稳定性的问题,前人开展了大量研究.大部分学者主要将克拉通岩石圈地幔视为一个均匀层来研究,发现较弱的克拉通岩石圈地幔容易被剥离减薄或拉伸变薄(Currie and van Wijk,2016;刘丹红和陈林,2019),而强流变性质和低密度能使克拉通岩石圈保持长期稳定(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)通过对双层克拉通岩石圈研究也得出了上述结论,但他们仅考虑了密度分层的影响.另外,研究者们考虑了克拉通岩石圈中弱低速层——中岩石圈不连续面(mid-lithosphere discontinuity, MLD;Selway et al.,2015)的影响后,发现克拉通中的弱MLD层与软流圈的接触范围会影响克拉通岩石圈减薄发生的难易程度,且减薄作用会沿着MLD分布的方向持续下去(Wang et al.,2018;Shi et al.,2020,2021),若MLD被厚而冷的岩石圈所包围,则克拉通岩石圈能够保持长期稳定性(Liu et al.,2018).这些研究大都建立在假设克拉通岩石圈地幔的整体成分和性质完全一致的基础之上,没有充分考虑克拉通岩石圈地幔物性随深度的变化.

越来越多的证据表明克拉通岩石圈地幔的成分和结构在垂向上是非均一的.一方面,在地球物理研究中：联合波形和SKS波分裂分析结合地球化学资料显示,整个北美大陆下方岩石圈存在～150 km厚的化学亏损的岩石圈层和其下的～30-120 km厚的热岩石圈根(Yuan and Romanowicz,2010);S波接收函数和表面波各向异性显示,Kalahari克拉通岩石圈地幔存在～50 km厚的S波各向异性层、～65-115 km厚的亏损层和～60-130 km厚的交代层(Sodoudi et al.,2013);各向异性Ps接收函数(RF)结果表明,在北美的Wyoming和Superior克拉通下方发现多个MLD,且深度和几何形状均不相同(Ford et al.,2016);地震照明成像(SDI)结果显示,西澳大利亚克拉通岩石圈地幔中存在多个MLD,且代表了岩石圈地幔上下强弱分层的边界层(Sun et al.,2018);调和分解的接收函数分析显示,澳大利亚的Yilgarn克拉通和北美的Superior克拉通的岩石圈至少分为两层,在170 km的深度以上存在多个相互独立的MLD(Chen et al.,2021);大地水准面的校验和残余地形及重力资料显示,大型克拉通下方的大部分过剩密度均分布在岩石圈地幔的下半部分,且克拉通岩石圈地幔的平均密度要比周围软流圈重0.5%～1.2%(Wang et al.,2022a,b).另一方面,在地球化学研究中：重矿物石榴石的主微量元素分析显示,加拿大Slave地区的岩石圈地幔在140～160 km深度分为上下两层,上方为横向不均匀的化学亏损层,下方为性质均匀的熔融交代层(Griffin et al.,2004);来自金伯利岩及相关岩石中的石榴石和铬铁矿捕虏体资料显示,大部分克拉通岩石圈地幔具有岩石较富集的上层(>110 km),岩石亏损的中下层(～110-190 km)和熔体交代的最底层(～170-220 km)(Perchuk et al.,2020).综上所述,大部分克拉通岩石圈地幔都存在垂向成分和结构的分层现象(图1b).

岩石圈地幔的成分和结构分层性是如何影响克拉通的稳定性?针对这一科学问题,我们聚焦于克拉通岩石圈地幔的垂向分层特性,采用二维热-力学的数值模拟方法,探讨在边界驱动对流背景下,克拉通岩石圈地幔的密度、黏度和成分分层结构对克拉通岩石圈稳定性的影响.本文先对动力学方法和模型设计进行介绍,然后系统测试在无MLD、连续MLD和间隔MLD条件下的克拉通岩石圈地幔密度-强度的分层结构,探讨在边界驱动对流作用下影响克拉通岩石圈长期稳定性的因素,最后讨论华北克拉通破坏时空范围的成因机制.

1 数值模拟方法

本文利用二维热-力学程序(I2VIS;Gerya and Yuen,2003;Perchuk et al.,2020),同时考虑了热扰动来模拟边界驱动对流(刘丹红和陈林,2019),研究边界驱动对流对克拉通岩石圈减薄的动力学过程.该程序的重点是粒子标记(marker-in-cell)技术和交错网格有限差分方法,联合求解动量、物质和能量守恒方程,不仅能考虑到局部和大尺度变形,还可以避免一定程度的数值发散.

1.1 控制方程

动量守恒方程：

(1)

(2)

式中,σ′xx、σ′yy和σ′xy为偏应力张量;x为水平坐标,y为垂直坐标;P为压力,T为温度,C为岩石成分,M为部分熔融百分数,上述参数综合控制着密度ρ;g表示重力加速度.

物质守恒方程：

(3)

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

能量守恒方程：

(4)

(5)

(6)

式中,Cp表示等压热容,文中取值为1000 J·kg-1·K-1;DT/Dt是温度关于时间的物质导数;qx为热流水平分量,qy为热流垂直分量;热传导率k受温度(T)、压力(P)和岩石成分(C)的综合影响(Clauser and Huenges,1995);α为热膨胀系数;Hr、Hs、Ha和HL分别代表放射性、剪切、绝热和潜热生热.

1.2 黏-塑性流变

岩石的流变性质受塑性和黏性流变机制的共同控制,偏应力张量和应变率张量的关系如下：

(7)

(8)

塑性变形遵从Drucker-Prager屈服准则(Ranalli,1995)：

(9)

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

(10)

(11)

韧性变形综合考虑了扩散蠕变和位错蠕变：

(12)

式中,岩石成分、温度和压力等共同控制着扩散蠕变黏度ηdiffusion和位错蠕变黏度ηdislocation(Ranalli,1995)：

(13)

(14)

式中,σcr表示扩散蠕变与位错蠕变发生转变的临界应力;R表示气体常数;物质常数(AD)、活化能(Ea)、活化体积(Va)和应力指数(n)是岩石的流变参数.

Peierls变形出现于高压低温区域(Karato et al.,2001)：

(15)

(16)

式中,σpeierls为限制物质强度的Peierls应力,与式(9)中的屈服应力σyield相似;σII表示应力张量二阶不变量;P、T、R、Ea、Va和n意义与式(14)相同,Apeierls为物质常数.

岩石圈的有效黏度ηeff为塑性黏度ηplastic、韧性黏度ηductile和Peierls黏度ηpeierls之间的最小值：

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

(17)

1.3 部分熔融

岩石的部分熔融百分数M与温度为简单的线性关系(Gerya and Yuen,2003;Burg and Gerya,2005)：

(18)

式中,Ts为岩石固相线的温度,Tl为岩石液相线的温度(计算表达式见表1).

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

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

(19)

式中,ρs为岩石固相的密度,ρm为岩石熔融的密度(计算表达式见表1).

1.4 长波长热扰动

为了模拟边界驱动对流,在软流圈中加入长波长热扰动(刘丹红和陈林,2019)：

(20)

式中,Tp为热扰动幅度,ylab和Tlab分别表示岩石圈底界面的深度和温度,y0和x0分别表示模型的垂向深度和水平宽度,y和x分别表示模型中的垂向位置和水平位置.

长波长热扰动的引入可以在岩石圈地幔下方形成横向不均一的温度结构,这种温度异常可以形成大尺度地幔对流,在不施加边界驱动力的情况,增加必要的热扰动可以有效地驱动模型演化(King and Anderson,1998).热扰动引发的大尺度地幔对流与克拉通-活动带岩石圈厚度横向急剧变化引发的小尺度边界驱动对流相互作用,可以起到抑制小尺度边界驱动对流过快发展的作用(刘丹红和陈林,2019).经过反复测试后,我们发现：边界驱动对流在Tp小于Tlab的1%时完全占主导,大尺度地幔对流在Tp大于Tlab的8%时完全占主导.因此,选择Tp=40 ℃(约Tlab=1350 ℃的3%)来设计软流圈温度场,既能模拟大尺度地幔对流,又有利于形成边界驱动对流.

2 模型设计

基于克拉通的分层结构以及克拉通-活动带的岩石圈厚度差异,我们设计了一个初始模型(图2).模型具有4000 km的水平宽度和1000 km的垂向深度(图2a),通过由2001×381个节点组成的非均匀网格系统来划分模型空间.在水平方向上,模型是2 km的平均网格分辨率.在垂直方向上,模型0～200 km深度的分辨率为1 km,200～300 km深度的分辨率从1 km逐渐降至5 km,300～1000 km的分辨率为5 km.如此,既可以保证克拉通岩石圈地幔处的分辨率,又可以提高计算效率.

图2 初始模型(a)模型设计;(b)为(a)中放大的模型区域,白色实线为等温线,黑色和红色实线分别表示0.5 Myr时刻岩石圈垂向剖面的密度和黏度.组分色标：0,黏性空气;5和6,大陆上下地壳;9,克拉通上层地幔岩石圈;10,软流圈;11,MLD;12,活动带地幔岩石圈;14,克拉通下层地幔岩石圈.

模型左侧和右侧的块体分别代表克拉通和活动带,两者的地壳性质一样,由均为20 km的上、下地壳组成.地壳下覆由40 km的上层岩石圈地幔、20 km的MLD和60 km的下层岩石圈地幔组成,克拉通岩石圈地幔总共厚120 km,活动岩石圈地幔性质均一,厚度为40 km,克拉通-活动带厚度差异为80 km(图2b).我们采用岩石力学实验确定的不同矿物的流动律(flow law)来描述模型不同成分层的流变学性质(如, Chen, 2021),以“湿石英”和“斜长石”分别代表上、下地壳的流变学性质,以“湿橄榄石”代表MLD的流变学性质,以“干橄榄石”代表克拉通和活动带地幔及软流圈的流变学性质(具体参数见表1).克拉通的上地壳和下地壳的地温梯度分别为15 ℃/km和10 ℃/km,活动带的上地壳和下地壳的地温梯度分别为20 ℃/km和12.5 ℃/km,克拉通和活动带的岩石圈底界面温度Tlab均为1350 ℃,软流圈中的绝热地温梯度K=0.5 ℃/ km.

模型的上边界和左右两侧边界为自由滑边界条件,下边界为渗透性边界条件(Burg and Gerya,2005),计算区域下方500 km处将符合自由滑条件,模型底部(y=1000 km)采用的渗透性边界条件,可以解决模型底部的部分边界效应问题.模型上边界温度为常温(0 ℃),左右两侧边界为绝热边界,下边界为外部温度固定的边界条件(Burg and Gerya,2005),计算区域下方500 km处将符合固定温度条件,该区域的热流和温度能够随着演化而自行调整.模型演化全程无外力扰动,全由初始状态驱动.模型顶部设置有20 km厚的伪空气层,其特点为低密度(1 kg·m-3)和低黏度(1019Pa·s),可模拟减薄过程中的地形演化(Schmeling et al.,2008).

3 模拟结果

在数值模拟过程中,我们将活动带和软流圈的参考密度和参考黏度保持不变.在改变克拉通岩石圈地幔的密度结构时,通过岩石圈地幔各层厚度来加权调节克拉通上、下层岩石圈地幔的密度,保持整个克拉通岩石圈地幔的参考密度与软流圈的参考密度基本一致.其中,厚度加权计算公式为：

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

(21)

式中,ρuclm、ρMLD、ρlclm、ρmlm和ρam分别为克拉通上层岩石圈地幔、MLD层、克拉通下层岩石圈地幔、活动带岩石圈地幔和软流圈地幔;huclm、hMLD、hlclm、hmlm和ham分别为克拉通上层岩石圈地幔、MLD层、克拉通下层岩石圈地幔、活动带岩石圈地幔和软流圈地幔的厚度.

本文系统测试了在边界驱动对流的作用下,克拉通岩石圈地幔的垂向密度-黏度分层结构及MLD的分布样式对克拉通岩石圈稳定性的影响.通过对不同的Δρ0、Δη0及MLD样式A和B进行组合,我们以无MLD情况下的Δρ0=40 kg·m-3和Δη0=50的模型Model-11为实验的参考模型,在此基础上开展了42组数值模拟实验(表2),进而研究在边界驱动对流的条件下,克拉通岩石圈地幔的不同分层结构对其稳定性的影响.

3.1 参考模型(无MLD)

图3为参考模型Model-11的演化结果.长波长热扰动驱使软流圈呈现自西向东的流动,流经克拉通和活动带的交界处.由于克拉通和活动带岩石圈厚度的突变,模型在演化到10 Myr时克拉通-活动带岩石圈台阶处出现了明显的边界驱动对流(图3a1,b1).较重的克拉通下层岩石圈地幔受边界驱动对流的影响出现对流侵蚀,被拆沉坠入软流圈并滞留于地幔过渡带(图3a2,b2).当模型演化至160 Myr,由于边界驱动对流的侵蚀,较重的克拉通下层岩石圈地幔自东向西横向减薄了～200 km,而较轻的上层岩石圈地幔只是轻微变形,岩石圈-活动带间的过渡带被剥蚀而变得平缓,地幔对流也随之减慢(图3a3,b3,c).无MLD的克拉通岩石圈地幔分层模型结果显示：由于克拉通下层岩石圈地幔比软流圈更重更强,克拉通只出现了边缘减薄现象,其整体岩石圈拥有较长的寿命,具有长期的稳定性,总体上与刘丹红和陈林(2019)呈现的结果一致.

3.2 MLD连续分布对克拉通稳定性的影响

3.2.1 克拉通岩石圈密度分层结构的影响

这里不改变克拉通下层岩石圈的参考黏度,系统测试克拉通下层岩石圈的参考密度对克拉通稳定性的影响.图4显示了模型Model-23(Δρ0=40 kg·m-3,Δη0=50)和模型Model-22(Δρ0=20 kg·m-3,Δη0=50)的对比结果,两个模型的克拉通下层岩石圈参考黏度均为软流圈参考黏度的50倍,仅在克拉通上、下层岩石圈地幔的参考密度上存在差异(图4d).模型Model-23显示：边界驱动对流的下降流会拆掉克拉通下层岩石圈地幔并使其坠入地幔过渡带,而上升流所附带的热物质会使MLD发生部分熔融,克拉通下层岩石圈的拆沉会在克拉通内部发育出新的岩石圈台阶,从而使对流减薄自东向西持续,在80 Myr时产生～1450 km的横向减薄范围(图4a1,a2,c).模型Model-22显示：当克拉通下层岩石圈的参考密度减轻20 kg·m-3时,边界驱动对流对克拉通的减薄影响减弱,减薄规模也相对较小,80 Myr的克拉通横向减薄范围仅有～150 km(图4b1,b2,c).在克拉通内MLD连续分布的情况下,克拉通岩石圈地幔内部密度分层的差异对比表明：当克拉通下层岩石圈更重时,其受到的边界驱动对流影响更大.

表2 数值模拟实验Table 2 Numerical modelling experiments

3.2.2 克拉通岩石圈强度分层结构的影响

为了研究克拉通岩石圈地幔强度分层结构的影响,克拉通岩石圈地幔的参考密度将保持不变,测试克拉通下层岩石圈地幔的不同参考黏度对克拉通稳定性的影响.图5显示了模型Model-23(Δρ0=40 kg·m-3,Δη0=50)和模型Model-15(Δρ0=40 kg·m-3,Δη0=1)的对比结果,两个模型的克拉通岩石圈地幔都保持相同的密度分层结构,仅在克拉通下层岩石圈地幔的黏度上存在差异(图5d).当克拉通下层岩石圈与软流圈的参考黏度差Δη0=50(Model-23)时,强的克拉通下层岩石圈地幔容易受边界驱动对流侵蚀的影响,产生大规模的块状拆沉,75 Myr时的减薄范围为～1000 km(图5a1,a2).而当克拉通下层岩石圈与软流圈无参考黏度差(Model-15)时,在75 Myr时克拉通岩石圈自东向西的横向减薄范围仅为～300 km,但持续的对流会使致密下层岩石圈地幔出现重力失稳现象,从而在岩石圈底界面继续产生向西～900 km的10～15 km厚的滴落状剥离(图5b1,b2).在克拉通内MLD连续分布的情况下,克拉通岩石圈地幔内部黏度分层的差异对比表明：当克拉通下层岩石圈更强时,其受到的边界驱动对流影响更大.

3.3 MLD不连续分布对克拉通稳定性的影响

大量的地震观测表明(Ford et al.,2016;Sun et al.,2018;Chen et al.,2021),MLD可能呈横向不连续分布,因此,我们在克拉通岩石圈地幔中设置了不同间隔距离的不同宽度的多个MLD,以测试MLD的不连续分布对克拉通稳定性的影响.

图3 参考模型的成分场和温度场演化白色箭头和实线分别表示地幔流场的方向和大小.

图4 克拉通岩石圈密度分层结构对克拉通稳定性的影响(Δη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 克拉通岩石圈强度分层结构对克拉通稳定性的影响(Δρ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的分布样式对克拉通稳定性的影响(Δρ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 克拉通岩石圈密度分层结构对克拉通稳定性的影响(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 克拉通岩石圈强度分层结构对克拉通稳定性的影响(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 岩石圈的不同分层结构对克拉通稳定性的影响模式图不同颜色的圆形代表克拉通岩石圈减薄的程度：绿色表示轻微减薄,蓝色表示横向减薄范围不超过200 km,红色表示横向减薄范围超过200 km;不同线段组合的白色圆形代表克拉通岩石圈减薄的减薄方式.

3.3.1 MLD分布样式的影响

保持克拉通岩石圈地幔中的密度(Δρ0=40 kg·m-3)和黏度(Δη0=50)分层结构不变,以此测试不同间距的MLD分布对克拉通稳定性的影响,模型结果显示：MLD分布较为密集的模型Model-25(A=50 km,B=10 km)的横向减薄规模小,且减薄速度较为缓慢(图6a1,a2),在100 Myr时的克拉通岩石圈横向减薄范围仅为～200 km(图6a3),而MLD分布较为稀疏的模型Model-27(A=200 km,B=10 km)的横向减薄规模较大,且减薄速度也较快(图6b1,b2),在100 Myr时的克拉通横向减薄范围可达～850 km(图6b3).在不连续MLD分布情况下,更稀疏的MLD分布更有利于边界驱动对流对克拉通的减薄破坏(图6a4,b4).

3.3.2 克拉通岩石圈密度分层结构的影响

以MLD分布较为密集的模型(A=100 km,B=10 km)为基础,保持克拉通下层岩石圈地幔的参考黏度(Δη0=50)不变,测试克拉通岩石圈地幔的不同密度分层结构对克拉通稳定性的影响.模型Model-32(Δρ0=20 kg·m-3)同Model-26(Δρ0=40 kg·m-3)比较显示,更轻的克拉通下层岩石圈的横向减薄规模小,且减薄速度慢(图7a1—a2,b1—b2),在120 Myr时,Model-32的克拉通岩石圈横向减薄范围比Model-26少1000 km(图7a3,b3).在克拉通内MLD分布较为密集的情况下,克拉通岩石圈地幔内部密度分层的差异对克拉通的稳定性影响同连续MLD分布下的结果趋势一致：重的克拉通下层岩石圈受到的边界驱动对流影响更大(图7a4,b4).

3.3.3 克拉通岩石圈强度分层结构的影响

继续保持相同的MLD分布样式(A=100 km,B=10 km),在克拉通下层岩石圈地幔的参考密度(Δρ0=20 kg·m-3)不变的情况下,测试克拉通岩石圈地幔的不同黏度分层结构对克拉通稳定性的影响.模型结果显示：克拉通下层岩石圈更强的模型Model-32(Δη0=50)同相对较弱的模型Model-38(Δη0=20)对比,岩石圈横向减薄规模更大,且减薄速度更快(图8a1—a2,b1—b2),在160 Myr时Model-38的克拉通岩石圈横向减薄范围仅为～500 km(图8a3),而Model-32能达～950 km(图8b3).在克拉通内MLD分布较为密集的情况下,克拉通岩石圈地幔内部黏度分层的差异对克拉通的稳定性影响也同连续MLD分布下的结果趋势一致：更强的克拉通下层岩石圈更容易受到边界驱动对流的影响而发生拆沉(图8a4,b4).

4 讨论

4.1 岩石圈地幔分层性对稳定性的影响

克拉通岩石圈地幔分层结构在俯冲构造背景下的岩石圈拆沉以及拆沉后再分层的演化过程在模拟实验中已有研究(如,Peng et al.,2022;Cao and Liu,2023),而图9显示了边界驱动对流背景下的克拉通岩石圈地幔内的垂向密度-黏度分层和MLD分布样式对其稳定性的影响.在不考虑MLD的情况下,模拟结果和前人研究相似：相对重的下层岩石圈地幔更容易受到外力而发生岩石圈减薄,其强度增加则有利于保持稳定(刘丹红和陈林,2019;Currie et al.,2023).而考虑MLD时,弱MLD的存在不仅可以导致大范围的减薄发生,并且会使克拉通减薄从其边缘沿着MLD分布的方向进行(Wang et al.,2018;Liu et al.,2018).其中,只有当克拉通下层岩石圈地幔同软流圈保持同等密度时,更强的克拉通下层岩石圈地幔才更有利于保持克拉通的稳定;而若克拉通下层岩石圈地幔更重时(Δρ0≥20 kg·m-3),边界驱动对流会导致致密的克拉通下层岩石圈出现明显的横向减薄,并沿着MLD的分布方向持续进行,此时,克拉通下层岩石圈地幔的强度增加则不利于克拉通的稳定性.

不同构造背景(俯冲、地幔柱)下的克拉通岩石圈中的MLD位置分布与外扰动产生的变形会影响克拉通岩石圈的稳定性(Shi et al.,2021,2022),而且在边界驱动对流的构造背景下,克拉通岩石圈中MLD的存在会使其减薄破坏呈现一个前期慢速减薄和后期快速减薄的非均匀过程(图10).当克拉通下层岩石圈地幔更重(Δρ0≥20 kg·m-3)以及更强(Δη0≥20)时,在克拉通岩石圈减薄过程中,其减薄方式会发生转变,具体体现在长期(～60-100 Myr)的慢速的小规模丝状减薄向短期(～10-20 Myr)的快速的大规模块状拆沉的过渡,且克拉通岩石圈前期的横向减薄范围少(～200-400 km),后期的横向减薄范围多(～600-800 km)(如, Model-23),此时,若增强下层岩石圈地幔的强度,则会使减薄方式转变的时间提前,即克拉通岩石圈的大规模块状拆沉提前发生,甚至克拉通岩石圈的减薄方式一开始就是块状拆沉(如, Model-24).而当克拉通岩石圈的MLD呈不连续分布(B=10 km)时,与MLD连续分布的情况下的克拉通岩石圈地幔密度和黏度分层结构的演化趋势一致,会有一个减薄方式转变的过渡,但由于MLD的横向不连续,减薄方式的转变时间往后延长了～30 Myr,并且会降低克拉通的横向减薄速度,从而延长克拉通的寿命(图10),且更密集的MLD分布样式有利于克拉通岩石圈抵抗边界驱动对流的侵蚀(图6).

图10 不同岩石圈地幔分层的克拉通岩石圈减薄随时间的水平范围虚线和实线分别表示在无MLD(黑色)、连续MLD(蓝色)和不连续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);地震接收函数显示,现今华北克拉通西部的岩石圈尚存200 km,而东部的岩石圈仅有60～80 km(Chen et al.,2008).岩浆岩的性质和同位素示踪研究显示华北克拉通的减薄可能始于晚石炭纪-晚三叠世(300～220 Ma; Yang et al.,2008;徐义刚等,2009),在早白垩世(约130～120 Ma)达到减薄峰期(姜耀辉等,2005;Hou et al.,2007),但并没有减薄到现今岩石圈的厚度(徐义刚,2004,2006).这表明华北克拉通东部岩石圈的厚度自古生代以来至少减薄100 km,而西部岩石圈尚保持完整.目前,主流观点认为华北克拉通东部的破坏是由太平洋板块俯冲引起的岩石圈地幔熔/流体交代作用及其伴随的海沟后撤引发的强烈伸展作用共同作用的结果(吴福元等,2000;Zheng et al.,2006,2007;Xu,2007;朱日祥等,2011,2012).目前认为华北克拉通的破坏呈现时空不均匀性的特点(徐义刚,2004;Zheng et al.,2006,2007),且存在不止一期的岩石圈减薄破坏,华北克拉通的减薄很可能是一个循序渐进的漫长过程(徐义刚等,2009;刘丹红和陈林,2019).而最新的华北克拉通镁铁质岩浆的平衡压力研究显示华北克拉通破坏的时间集中在早白垩世(～120 Ma)期间,华北克拉通的部分东部岩石圈在极短的时间内(10 Ma)从200 km减薄到35 km(Chen et al.,2023),这种短时间内的大规模减薄难以用太平洋板块俯冲模型来解释.岩石地球化学证据表明：华北克拉通破坏前的岩石圈地幔根部的黏度是软流圈的2～8倍(Wang et al.,2022),在华北克拉通中部的岩石圈地幔中有上部的残留古老地幔和下部的改造地幔的双层特征(汤艳杰等,2021);地球物理证据表明：在现今的华北克拉通西部陆块和中部跨华北造山带的岩石圈地幔中分别观测到MLD的存在,且MLD存在的深度与华北克拉通东部陆块的LAB深度一致(Chen et al.,2014).综上所述,华北克拉通在破坏前存在明显的岩石圈地幔成分和结构的分层,并且破坏前的东部岩石圈地幔中可能存在MLD,但与中部岩石圈地幔中的MLD间隔较远.

MLD不连续分布的岩石圈分层模型结果显示(图9),在开始的很长一段时间内(90～130 Myr),克拉通下层岩石圈地幔只有其边缘向西200～400 km的局部区域发生减薄,之后仅约30 Myr就能完成继续向西约700 km的大规模减薄,这与上述的破坏时间集中的特征一致,这表明华北克拉通岩石圈地幔中很可能存在横向MLD不连续分布的垂向密度-黏度分层结构.因此,我们提出一个可能的华北克拉通破坏过程的新模型：(1)130 Ma之前减薄开始：华北克拉通的东部岩石圈由于受到周围的年轻大洋板块或造山带引起的边界驱动对流,其边缘的MLD刚开始受到的外力扰动影响不大,岩石圈产生规模较小、速度较慢且持续时间很长的丝状减薄(图11a);(2)130 Ma到120 Ma期间减薄进入峰期：当华北克拉通东部岩石圈中的横向宽度较窄的单个MLD被完全拆掉后,持续的对流开始剥离相邻且横向宽度较长的MLD,岩石圈减薄由丝状剥离转变为块状拆沉,导致华北克拉通岩石圈在短时间内产生大规模的减薄(图11b);(3)120 Ma至现今减薄减缓或结束：当岩石圈减薄到达相邻MLD间隔距离很大的华北克拉通中部时,对流减薄作用很难影响到华北克拉通西部内被厚岩石圈包围的MLD,岩石圈减薄方式又转变为丝状剥离,且减薄速度减缓,直至形成现今位于华北克拉通中部的平缓岩石圈台阶(图11c).由此,我们认为克拉通岩石圈地幔的垂向密度-强度分层结构和MLD的分布样式在华北克拉通破坏过程中发挥了重要作用.

5 结论

本文采用二维热-力学数值模拟方法探究在边界驱动对流条件下,克拉通岩石圈地幔的垂向密度-强度分层结构及其内部MLD的分布样式对克拉通稳定性的影响,获得的主要结论如下：

(1)当克拉通岩石圈内部不存在MLD时,克拉通岩石圈减薄的方式主要为丝状剥离,在下层岩石圈地幔较轻(Δρ0≤40 kg·m-3)且较强(Δη0≥1)的情况下,其减薄水平范围小于100 km,且主要发生在克拉通边缘处.

(2)当克拉通岩石圈内部存在一层连续的MLD时,在克拉通下层岩石圈地幔较重(Δρ0>0 kg·m-3)的情况下,克拉通岩石圈减薄的水平范围随着克拉通下层岩石圈地幔密度的增加而增大,若克拉通下层岩石圈较弱(Δη0=1),岩石圈地幔的减薄受对流侵蚀和重力失稳的共同影响;若克拉通下层岩石圈较强(Δη0≥1),岩石圈地幔的减薄仅受对流侵蚀影响,其减薄方式随着其强度的增加从缓慢丝状剥离转变为快速块状拆沉,甚至从开始阶段就越过丝状剥离而直接变成块状拆沉.

(3)当克拉通岩石圈内部存在一层不连续分布的MLD时,在克拉通下层岩石圈地幔的密度和黏度不变的条件下,单个MLD的宽度A=100 km且相邻两个MLD之间的间隙B=10 km的不连续分布样式能够将减薄方式改变的时间延后～30 Myr,并降低减薄的速度,有效延缓岩石圈地幔的减薄.间隙越稀疏岩石圈减薄越快,反之亦然,减薄方式早期表现为缓慢丝状剥离,中后期表现为快速块状拆沉.

(4)华北克拉通破坏具有持续时间长(>100 Ma)、峰期集中(130～120 Ma)的特点,并且其破坏的前缘位置与中部造山带下方MLD的东端对应.根据模拟结果,我们认为华北克拉通破坏的时空范围呈现的不均一现象,可能主要是由岩石圈地幔的垂向密度-强度分层结构和MLD不连续分布引发的岩石圈地幔逐段拆沉所导致.

图11 华北克拉通分期破坏过程的示意图(a) 130 Ma之前,华北克拉通岩石圈开始减薄,减薄规模小、速度慢; (b) 130 Ma到120 Ma期间,华北克拉通岩石圈进入减薄峰期,减薄规模大、速度快; (c) 120 Ma至现今,华北克拉通岩石圈减薄缓慢,甚至趋近于停止.

致谢谨以此文祝贺滕吉文先生90华诞暨从事地球物理工作70年.感谢瑞士苏黎世联邦理工学院Taras Gerya教授提供的I2VIS程序.感谢审稿人对本文提出的建设性意见.感谢中国科学院地质与地球物理研究所田小波研究员、颜智勇博士、崔晓娜博士、司润港博士对本文的建议和帮助.本文的模拟计算是在国家超级计算天津中心和北京超级云计算中心完成.