2016, Vol. 59
2. 蒙城地球物理国家野外科学观测研究站, 合肥 230026
2. Mengcheng National Geophysical Observatory, Hefei 230026, China
克拉通是大陆上相对稳定的构造单元，包括地壳和岩石圈地幔(Rudnick and Fountain，1995;Schubert et al.，2001).地质、地球物理和地球化学研究显示：克拉通有200多公里的岩石圈根，有的厚达300到400 km；克拉通区域表面热流值小；岩石圈根相对较冷、含水量低，并且显示出是部分熔融分异后的残留物；岩石圈根的密度相对较小、黏性较高；克拉通形成后可长期保持稳定，有的可长达40亿年(Boyd et al.，1985;Carlson et al.，2005;King，2005;Pearson，1999;Pearson et al.，1995).但是有些克拉通(如美国的Wyoming克拉通、我国的华北克拉通)在其形成后遭遇过改造(Carlson et al.，2004;Gao et al.，2002).其中华北克拉通是目前研究最为广泛和深入的地区.在华北克拉通发现了38亿年的陆壳残留物(Liu et al.，1992).在经历了早元古代(约25—18亿年间)的增生聚合碰撞事件后随着岩石圈厚度增大，华北克拉通作为一个完整的大陆块体长时期处于相对稳定的状态(Zhao et al.，2001).金伯利岩研究显示约在早古生代(~4.6亿年前)华北克拉通地幔岩石圈根的厚度约为200 km(路凤香等，1991).大约在中侏罗世(~1.6亿年)中国东部开始大规模火山喷发.目前的初步研究认为华北克拉通岩石圈减薄发生在晚中生代，在早白垩世(1.2—1.3亿年)达到高潮，岩石圈的厚度从200 km减薄到80～120 km左右.伴随着这一过程，大量的火山活动使地幔岩石浸入地壳，地壳运动非常活跃，导致古老的华北克拉通的解体或破坏(Griffin et al.，1998;Menzies et al.，1993;Wu et al.，2005;Xu，2001;邓晋福等，2003;范蔚茗和Menzies，1992;李江海等，2004;邵济安等，2002;吴福元等，2003;许文良等，2004;赵越等，2004).
尽管目前对造成克拉通岩石圈减薄的机制和过程仍有许多争论，但普遍认为对流减薄在克拉通岩石圈减薄过程中具有重要影响.克拉通岩石圈减薄的概念性模型包括拆沉、对流侵蚀、热侵蚀、橄榄体和熔体的相互作用、岩浆提取、岩石圈地幔水化、以及机械拉张等(吴福元等，2008).不同机制导致的岩石圈减薄后的后果和地表反应存在区别(吴福元等，2008;朱日祥等，2012).各种概念模型都有相应的地质、地球物理和地球化学观测等数据支撑(Carlson et al.，2005;Duggen et al.，2003;King，2005;Schurr et al.，2006;Zhao and Zheng，2007;吴福元等，2008;许文良等，2004).但对华北克拉通，目前的研究结果显示，由于岩石圈俯冲的影响，华北克拉通区域及其下覆地幔受到影响，软流圈地幔和岩石圈地幔的相互作用是导致华北克拉通被破坏的主因(朱日祥等，2011，2012).因而研究软流圈地幔对流对岩石圈地幔的影响具有重要意义(朱日祥等，2012).
研究地幔对流对岩石圈的影响一直是地球动力学研究所关心的话题(Buck and Toksöz，1983;Davies，1994;Houseman et al.，1981;Huang et al.，2003;Marotta et al.，1999;Morency et al.，2002;Olson et al.，1988;Schott and Schmeling，1998;Shapiro et al.，1999;Yuen and Fleitout，1985;乔彦超等，2012，2013).数值模拟结果显示，加厚岩石圈减薄的过程和时间与模型参数密切相关.不同模型由于关注的对象不同，模型参数各异，导致得到的岩石圈减薄时间尺度变化很大，从~10 Ma(Marotta et al.，1999)到几百个Ma(Buck and Toksöz，1983;Houseman et al.，1981)、甚至能够长期保持稳定(Shapiro et al.，1999).因而，定量研究加厚岩石圈的对流减薄和模型参数的关系，显得十分重要.
Morency等(2002)通过数值计算和理论分析的方法，定量分析了加厚岩石圈根在对流地幔中被移除的时间.其结果显示，加厚岩石圈根的移除时间与其大小(即其二维模型中的加厚岩石圈根的厚度和宽度)和黏性参数密切相关，并给出了估算公式.这对讨论增厚岩石圈的减薄过程及其对地表的影响等具有重要意义.然而在Morency等(2002)的模型中并未考虑克拉通地幔岩石圈和普通地幔岩石圈之间可能存在的由化学组分不同引起的密度和黏性差异，而这正是克拉通岩石圈和一般地幔岩石圈间的差异所在.现有研究都显示，克拉通岩石圈的密度和黏性差异对其稳定性具有重大影响(Lenardic and Moresi，1999;Lenardic et al.，2003;O'Neill et al.，2008;Wang et al.，2014;Yoshida，2012).因而本文将在Morency等(2002)的定量化方法的研究基础上，通过在数值模型中增加考虑由化学组分引起的岩石圈的密度和黏性变化的影响，进一步研究克拉通岩石圈在对流地幔中的减薄过程和减薄时间，为讨论克拉通岩石圈的对流减薄提供数值参考模型.
2 数值模型克拉通岩石圈的对流减薄可采用热化学地幔对流方程组来描述，该方程组包括质量守恒方程、动量守恒方程、能量守恒方程和化学组分方程(Moresi and Solomatov，1995;Zhong et al.，2000).在Boussinesq近似下，该方程组的无量纲形式为：
模型假定地幔黏性与温度相关，其无量纲形式为：
数值模型是3×1的二维模型，方程(1)采用二维有限元程序Citcom解算.模型划分为289×97个有限单元网格，在靠近上下边界处，适度进行了网格加密.上表面给定温度，下表面热流为0，上表面固定，下表面和侧面边界自由滑移.本文模型深度取值1000 km，一方面能反映出软流圈与下地幔的黏性跳变，另一方面则是增大分辨率，模型选用的其他参数值见表 1.根据(2)式，模型瑞利数Ra=1.3×107.由于我们研究的是拟稳状态下的对流模型，考虑到地球的冷却效应，为保证地幔温度保持大约1250 K，模型选取的内生热略大于地球内生热(Turcotte and Schubert，2002).
 | 表 1 模型参数值 Table 1 Physical parameters |
为了定量研究克拉通岩石圈根的对流减薄时间，这里参考Morency等(2002)的做法，即：(1)首先计算一个参考模型，使其达到拟稳状态，此时模型的表面热流、内部温度以及均方根速度等基本保持不变.由于E很大，模型表现为近似刚性盖层下的小尺度对流形式(图 1(a—d)).该边界层厚度，Zq，用以刻画岩石圈厚度，其值按横向平均温度T≤0.9Ti计算，Ti定义为无量纲模型深度在0.3到0.8之间的温度平均值.(2)克拉通岩石圈根通过增加参考模型局部区域岩石圈厚度来模拟(图 1(e—h)).具体方法是在模型中部选取一定宽度(x)的区域将该区域厚度从Zq增加到Z0，Z0=(γ+1)Zq.由于边界层基本处于热传导状态，内部温度近似成线性关系，增厚边界层内的温度随深度的变化关系和增厚前一样，只是系数减小1/(γ+1).此时，增厚岩石圈相比于初始时，厚度增加了γZq，这里我们将γ称之为增厚倍数.(3)本文的一个重要任务是确定岩石圈对流减薄的时间.图 1(e，i—l)分别给出了模型caseB3在岩石圈加厚后的温度场和流场的时间演化过程.从图中可以看到，经过一定时间的演化，增厚岩石圈厚度从(γ+1)Zq基本重新减薄到Zq.增厚岩石圈根移除时间的计算方法仍采用Morency等(2002)的方法，即计算每一时间步加厚岩石圈区域在深度Z=0.9Zq处的横向平均温度进而给出时间变化曲线(图 2)，然后利用温度演化公式(Morency et al.，2002)：
 | 图 1 参考模型在拟稳状态下的等温线和流线(a)以及横向平均温度(b)、黏性(c)和均方根速度Vrms(d).模型caseB3在初始状态下的等温线和流线(e)以及横向平均温度(f)、黏性(g)和Vrms(h) (f—h)中虚线是整个模型的横向平均值，实线是岩石圈加厚部分的横向平均值.(i—l)显示加厚岩石圈演化过程的温度场和流场.(a，e，i—l)中，实线代表等温线，相邻等温线间距为无量纲温度0.1，黑色点线表示流函数值小于零的流线，黑色虚线表示流函数值大于等于零的流线. Fig. 1 Isotherm, streamline (a) and radial dependences of horizontal average temperature (b), viscosity (c) and Vrms (d) at the quasi stable state for the reference model. Isotherm, streamline (e) and radial dependences of horizontal average temperature (f), viscosity (g) and Vrms (h) at the initial state of the Model caseB3 The dashed line are displayed for the horizontal average value of the whole model and the solid line are displayed for the horizontal average of the thickened lithosphere in (f—h). (i—l) The temperature and flow field evolution of the thickened lithosphere. Solid line are displayed for isotherm, the spacing of adjacent isotherm is dimensionless temperature 0.1, black dotted line for positive streamline and black dashed line for negative streamline in (a, e, i—l). |
 | 图 2 深度Z=0.9Zq处的横向平均温度的时间演化曲线 不同线型(点线除外)表示在不同的增厚倍数γ下的数值计算值，而通过对应线型的点线代表对应数值计算值的拟合曲线. Fig. 2 Horizontal average temperature evolution with time at Z=0.9Zq The line (except dotted line) in different linear are displayed for the numerical values of the different thickening factors γ and the dotted line corresponds to the line in different linear are displayed the fitting curves corresponds to the numerical values. |
| | 表 2 模型参数和结果 Table 2 Model parameters and results |
| | 表 3 模型参数和结果 Table 3 Models parameters and results |
增厚岩石圈的对流减薄与岩石圈的物性参数密切相关(Huang et al.，2003;Morency et al.，2002).瑞利数和活化能等对岩石圈减薄过程具有重要影响(Huang et al.，2003;Morency et al.，2002)，但由于Morency等(2002)在其工作中详细分析和讨论了Ra(瑞利数)和E(活化能)对增厚岩石圈移除时间的影响，为能充分讨论克拉通岩石圈组分差异导致的密度和黏性的影响，本文将不再考虑Ra和E变化的影响，而是将其作为模型给定参数(参见表 1).根据表 1的模型参数，本文参考模型在处于拟稳状态时的热边界层厚度Zq =0.12.下面首先将给出B=0和ηc=1的一组模型的结果，该组模型忽略克拉通岩石圈固有的密度和黏性差异.然后通过改变B和ηc值研究克拉通岩石圈固有的密度和黏性差异的影响.即在讨论并总结一般岩石圈对流减薄规律的基础上，进一步讨论分析克拉通岩石圈对流减薄的规律.
3.1 加厚岩石圈根的对流减薄对一般加厚岩石圈，在本文模型中即假定B=0和ηc=1.该组模型与Morency等(2002)的模型基本一致.但应该注意，与Morency等(2002)模型不同的是，为和实际地球更加接近，本文考虑了660 km处的黏性跃变(Hager and Richards，1989).表 2给出了γ变化从0.5到2和x变化从0.4到3的共20个算例的结果.
图 1(e，i—l)显示了模型caseB3的演化过程，该模型增厚岩石圈的宽度为1.5，增厚倍数为1.0.图 1e是模型caseB3在初始时刻的等温线和流线.从图中可以看到，增厚岩石圈的减薄过程受到侧面和底部对流的影响(图 1i和1j).在两侧以及底部对流的同时作用下，增厚的岩石圈逐步被移除、抹平(图 1(i—k)).在0.0073时刻，增厚部分的岩石圈基本被移除(图 1k).这与表 2给出的、依据公式(4)和(5)计算得到的模型caseB3的移除时间0.0072基本一致.下面将详细分析增厚岩石圈移除时间与γ和x的关系.
图 3a给出了对流移除时间与增厚倍数γ的对数关系图.图中相同形状表示x相同的模型，连接相同形状的直线是通过线性拟合得到的相同x值下增厚岩石圈移除时间和增厚倍数的关系.从图中可以看到：对相同的x，对流移除时间与增厚倍数基本成幂函数关系，对流移除时间随增厚倍数增大而增大；对不同的x，拟合直线的斜率略有差异，但基本一致.计算结果显示：
 | 图 3 增厚岩石圈对流移除时间与增厚倍数γ(a)和宽度x(b)的关系；(c) 模型计算结果和拟合结果的比较 Fig. 3 Convective removal duration of the thickened lithosphere scale with the thickening factor γ (a) and width x (b)；(c) Comparison of the model results and the fitting results |
图 3b给出了对流移除时间与增厚宽度x的对数关系图.图中相同形状表示γ相同的模型.连接相同形状的直线是线性拟合得到的相同增厚倍数时，移除时间和增厚岩石圈宽度的关系.从图中可以看出，对流移除时间和增厚宽度也基本成幂函数关系，具体计算结果显示：
图 3c给出了根据(8)式计算的数值和表 2中给出的模型结果的比较.从图中可以看到，式(8)可以很好地表示模型的计算结果，只是在移除时间很大和很小时略微存在差异.
根据(8)式可以方便地估算一个增厚岩石圈被对流减薄的时间.对于初始平衡厚度为120 km的岩石圈，将其厚度增加到300 km(即γ=1.5)，当增厚岩石圈宽度为300 km(x=0.3)时，增厚岩石圈根被移除的时间为225 Ma；而当加厚岩石圈宽度为1500 km(x=1.5)时，则其增厚岩石圈根的移除时间是342 Ma；而对于更大的加厚岩石圈宽度，如10000 km(x=10.0)，其根的移除时间是559 Ma.
尽管本文模型的黏性结构与Morency等(2002)的模型略有差别，但两者的结果仍具有一致性.Morency等(2002)通过数值计算和理论分析认为，加厚岩石圈的减薄包括侧面和底部的侵蚀影响，当加厚岩石圈宽度很小时，以侧面侵蚀为主；而当加厚岩石圈宽度很大时，以底部侵蚀为主.忽略次要影响的前提下，他们给出的无量纲关系式为：
如前面提到的，本文采用在深度大于660 km的下地幔黏性增加30倍来表示黏性随深度的变化，而Morency等(2002)采用活化体积V*来表示黏性随深度的变化.要达到岩石圈下覆软流圈的黏性与下地幔的平均黏性相差约30倍，V*的取值大约应为0.75 cm3·mol-1，相应的无量纲的活化体积V=2.1621.根据表 1及表 2的数据，通过式(12)和(13)计算得到ΔTe=0.183.将此V及ΔTe值和本文模型的其他参数值代入(9)—(11)式，并取γ=1.5可计算得到无量纲的时间：
从上面的结果中可以看到，Morency等(2002)的结果与本文模型结果存在一定差异.这主要是由于Morency等(2002)在分析数据时的理论模型过于简单所导致的.如图 1所示，增厚岩石圈根的移除是侧面和底部对流侵蚀同时作用的结果，但在Morency等(2002)的模型中，以特征宽度为界限，宽度小时忽略了底部的作用，而宽度大时又忽略了侧面的作用，所以他们的结果显示在宽度小时，移除时间与宽度成正比，而宽度较大时，移除时间与宽度无关.
3.2 克拉通岩石圈根的对流减薄本文关注的焦点是克拉通岩石圈根的对流减薄，这就涉及到克拉通岩石圈的密度和黏性参数的变化.因此，研究克拉通岩石圈根的对流减薄，即使不考虑Ra和E的影响，模型至少也包含有与克拉通岩石圈根的大小(γ，x)和化学组分导致的密度与黏性变化(B，ηc)有关的四个参量.前面已经研究了不考虑化学组分导致的密度与黏性变化情况下的增厚岩石圈根在热对流下的对流减薄问题，并得到了对流减薄时间和增厚的宽度以及厚度的量化关系.这里将在此基础上，通过改变B和ηc来模拟化学组分变化导致的克拉通岩石圈的密度和黏性变化，以讨论它们对克拉通岩石圈根的对流减薄的影响.表 3给出了这组模型的计算结果.
比较表 2和表 3中移除时间的大小，我们发现在模型中增加克拉通岩石圈的密度和黏性变化后，增厚岩石圈根更难以被移除.但是，即使考虑克拉通岩石圈的密度和黏性的影响，当B和ηc较小时，岩石圈根的移除过程和3.1节模型所示的普通岩石圈根移除过程基本类似(如图 4).岩石圈根的移除仍然是侧面和底部对流减薄共同作用的后果.图 4显示了模型caseB3_7的温度、流场和克拉通岩石圈组分随时间演化的过程.
 | 图 4 caseB3_7在不同时刻的(a—f)等温线及流线和(g—l)组分场演化 (a—f)中实线代表等温线，相邻等温线间距为无量纲温度0.1.黑色点线表示流线为负值，黑色虚线表示流线为正值. Fig. 4 The isotherm, streamline (a—f) and composition (g—l) evolution of the model caseB3_7, displayed at different time Solid line for isotherm, the spacing of adjacent isotherm is dimensionless temperature 0.1. Black dotted line for positive streamline and black dashed line for negative streamline in (a—f). |
从图 4可以看到模型温度场的演化和普通增厚岩石圈模型的(图 1(e，i—l))演化相似，底部和侧面的对流侵蚀着岩石圈根，使其慢慢变薄，直至重新达到稳定.从组分场(图 4(g—l))也可以看到，克拉通岩石圈组分被侧面和底部的对流带到了下部地幔并混入其中，增厚的克拉通岩石圈根逐渐被移除.表 3显示caseB3_7的移除时间为0.0266，从图 4中也可以看到，图 4f和图 4l对应的时间为0.0259，此时加厚的克拉通岩石圈根几乎被完全移除.
但当B或ηc较大时，移除时间将显著增长(表 3).图 5显示的是caseB3_14、caseB3_15、caseB3_16、caseB3_20和caseB3_23的组分时间演化图.该组模型的克拉通岩石圈根大小(γ=1.0，x=1.5)相同，克拉通组分变化导致的黏性变化ηc以及浮力数B的数值如图所示.从图 5中可以明显看到，当B=0.2时，随着ηc的增大，增厚的岩石圈根的移除时间不断加长(图 5中的第1，2，3行)，这也清楚地显示在岩石圈内温度场的演化上(图 6).图 6给出模型caseB3_14-17、caseB3_20和caseB3_23在深度Z=0.9Zq处的横向平均温度演化曲线.一般来说，该处的温度演化基本满足(4)式的指数规律(图 2).图 6中显示，在B=0.2且黏性比较小时(caseB3_14和caseB3_15)，平均温度曲线基本满足(4)式.拟合得到的增厚的岩石圈根移除时间分别为0.0282(894 Ma)和0.0503(1.60 Ga)(表 3).但当黏性比较大时，在0.1的时间内，Z=0.9Zq处的横向平均温度一直处于上升阶段(图 6中caseB3_16，ηc=10).这表明，在0.1的时间内克拉通岩石圈并没有被充分移除，这也反映在图 5的组分场上(图 5中的第3行).
 | 图 5 不同组分黏性比以及浮力数下，增厚的克拉通岩石圈根的组分场演化.增厚克拉通岩石圈根的增厚倍数γ=1.0，增厚宽度x=1.5 Fig. 5 The composition of the thickening cratonic lithosphere, after stretching factor γ=1.0, width x=1.5, displayed at different viscosity ratio and buoyance number |
 | 图 6 模型caseB3_14-17、caseB3_20和caseB3_23在深度Z=0.9Zq处的横向平均温度演化曲线 实线、虚线、点划线、点线以及较粗的实线、虚线分别表示caseB3_14—17、caseB3_20和caseB3_23. Fig. 6 Horizontal average temperature evolution with time at Z=0.9Zq Solid，dashed，dot-dashed，dotted line and the thicker solid，dashed line are displayed for caseB3_14—17，caseB3_20和caseB3_23. |
图 5还显示，当B或ηc足够大时(如ηc=10；B=0.3和0.4)，克拉通岩石圈的对流减薄过程与前述过程表现的非常不一样(图 5中的第4，5行).图 5显示，此时很少克拉通岩石圈组分被带入其下的地幔中，克拉通岩石圈组分在漫长的演化过程中，向两边扩散，这也会使得克拉通岩石圈得到部分减薄(图 5中的第4，5行)，但这一过程非常漫长.从模型caseB3_20和caseB3_23在深度Z=0.9Zq处的横向平均温度演化曲线(图 6中的较粗的实线和虚线)上看，平均温度在0.1时间内仍处于上升阶段.很显然，如果平均温度一直处于上升，就无法利用(4)式进行拟合，从而无法确定对流减薄时间.无量纲0.1的时间大约相当于3.17 Ga，这个时间非常长，如果模型在这么长时间内无法计算对流移除时间，我们就不将这类模型纳入后面的统计分析中.
为分析克拉通岩石圈密度变化的影响，这里先定义一个参数等效密度变化Δρtc*，即
图 7a给出的是增厚区域宽度一定条件(x=1.0)下对流移除时间随着岩石圈增厚倍数γ的对数关系图，其中相同形状表示相同的组分黏性比ηc(圆圈表示ηc=1，方形表示ηc=3，菱形表示ηc=10)，而相同的线型则代表着相同的浮力数B(实线表示B=0，点划线表示B=0.05，虚线表示B=0.1)，通过相同形状的线条是其线性拟合.从图中能够看到：总体的趋势是对流移除时间随着增厚倍数的增大而增大，这和热对流的结果一致.同时，还可以看到，图 7a中实线相连的菱形(caseB4_2、caseC4_2和caseD4_2)，点划线相连的菱形(caseB4_4、caseC4_5和caseD4_5)以及虚线相连的菱形(caseB4_7、caseC4_8和caseD4_8)对应的三条拟合直线几乎平行，其斜率相差不大，平均值约为0.33.实线相连的方形(caseB4_1，caseC4_1和caseD4_1)，点划线相连的方形(caseB4_3，caseC4_4和caseD4_4)以及虚线相连的方形(caseB4_6，caseC4_7和caseD4_7)对应的三条拟合直线也几乎平行，其斜率相差不大，平均值约为0.56.实线相连的圆圈(caseA4、caseB4、caseC4和caseD4)以及虚线相连的圆圈(caseB4_5、caseC4_6和caseD4_6)对应的两条拟合直线几乎平行，其斜率相差不大，平均值约为0.76.这说明对相同的ηc，对流移除时间随着增厚倍数的增大而增大，对流移除时间与增厚倍数基本成幂函数关系，B的影响不大(图 7a中的相同形状的实线、虚线和点划线的斜率大致相同)；但对流移除时间与增厚倍数的关系的幂指数随着ηc的增大而减小(图 7a中圆圈、方形和菱形对应的斜率，即θ，具体值见表 4)，这说明对流移除时间与增厚倍数的关系与ηc相关.根据表 4的结果，可以得到对流移除时间与增厚倍数的关系的幂指数θ与黏性变化ηc的关系(图 7b).
 | 图 7 (a)对流移除时间和岩石圈根增厚倍数γ的量化关系；(b)量化关系参数θ随着组分黏性比ηc的变化关系；(c)去除了增厚倍数γ的影响后，对流移除时间随着组分黏性比ηc的量化关系；(d)去除了增厚倍数γ和组分黏性比ηc的影响后，对流移除时间随着Δρtc的量化关系；(e)去除增厚倍数γ，组分黏性比ηc以及Δρtc的影响后，对流移除时间随着增厚宽度x的量化关系；(f) 模型计算结果和拟合结果的比较 Fig. 7 (a) Convective removal duration scale with the lithospheric root thickening factor γ; (b) Scaling law parameters θ varies with the composition viscosity ratio ηc; (c) Convective removal duration scale with the composition viscosity ratio ηc without the effect of the lithospheric root thickening factor γ; (d) Convective removal duration scale with the Δρtc without the effect of the lithospheric root thickening factor γ and the composition viscosity ratio ηc; (e) Convective removal duration scale with the width x without the effect of the lithospheric root thickening factor γ, the composition viscosity ratio ηc and the Δρtc; (f) Comparison of the model results and the fitting results |
| | 表 4 γ的量化关系参数 Table 4 Scaling law parameters for γ |
由图 7a我们知道增厚倍数γ对于对流移除时间的影响很大，且当岩石圈的组分黏性比ηc不同时，影响的幅度也不同，所以在下面的分析数据中，我们都将增厚倍数对移除时间的影响利用(20)式去除后再来研究.去除增厚倍数的影响后，我们发现岩石圈的组分黏性比对移除时间的影响占主导地位(图 7c—7e).图 7c给出的是去除了增厚倍数γ的影响后，对流移除时间随着组分黏性比ηc的对数关系图，图中的圆圈包含表 3中所有的B≤0.1的模型和表 2中的caseB1、caseB2、caseB3、caseB4、caseC4和caseD4.图 7c中包含18组数据.每组增厚倍数、增厚区域的宽度以及等效密度变化都相同.从图 7c可以看到增厚倍数、增厚区域的宽度以及等效密度变化的影响几乎难以区分.因此在获取对流移除时间和ηc的幂指数时，我们采取将所有相同的ηc的模型的对流移除时间求平均，然后进行拟合，据此得到α=0.52，即
和3.1节一样，在讨论对流移除时间与等效密度变化Δρtc的关系时，也先将对流移除时间中包含的γ和ηc影响先行去除.图 7d给出了去除了增厚倍数γ以及组分黏性比ηc影响后，对流移除时间随着等效密度变化Δρtc的变化规律，图中的圆圈同样包含了表 3中所有的B≤0.1的模型和表 2中的caseB1、caseB2、caseB3、caseB4、caseC4和caseD4.如图 7d所示，随着热边界层总的等效密度变化的增大(浮力数的减小)，整体上的变化趋势是相同的，都表现出了对流移除时间在减小，这里采用对所有数据点进行线性拟合的方法确定(18)式中的β.拟合结果给出β=-0.21，故
最后，在得到(20)—(22)式的基础上，通过将对流移除时间进行归算，即扣除增厚倍数γ、组分黏性比ηc以及等效密度变化Δρtc的影响后来研究它和增厚宽度x的关系，如图 7e所示.图 7e给出的是去除了增厚倍数γ、组分黏性比ηc以及等效密度变化Δρtc的影响后，对流移除时间随着增厚区域宽度x的变化规律，图中的圆圈同样包含了表 3中所有的γ=1.0且B≤0.1的模型和表 2中的caseB1、caseB2、caseB3和caseB4.从图 7e中可以看到：利用(20)—(22)式归算后，对固定的x，对流移除时间基本一致.因此，在获取(18)式的μ值时，我们再一次采用将相同x时的数据取平均的方法.根据平均值数据拟合得到μ=0.04，即
图 7f给出了根据(24)式计算的移除时间和数值模型计算得到的移除时间之间的比较，从图 7f中可以看到，(24)式可以很好地估算模型的计算结果.
式(24)可用来估算克拉通岩石圈被对流减薄的时间.采用3.1节相同的假定，即假定300 km(即γ=1.5)厚的克拉通岩石圈，在对流的作用下减薄到初始平衡厚度120 km.若克拉通岩石圈的密度比周围地幔的密度低0.4%(即B=0.1或Δρtc=0.082)且组分黏性比ηc=10，在克拉通岩石圈宽度为300 km(x=0.3)时，其对流减薄时间大约为1.11 Ga；在克拉通岩石圈宽度为1500 km(x=1.5)时，其对流减薄时间增加到大约1.18 Ga.对比3.1节结果可以看到，0.4%的组分密度变化和10倍的组分黏性变化，使得岩石圈的移除时间从225~342 Ma增加到1.11~1.18 Ga.这显示，克拉通岩石圈组分变化对移除时间具有重要影响.
保持其他参数不变，仅改变克拉通岩石圈组分浮力或黏性比，可以估算组分浮力或黏性比的影响.如仍假定300 km厚(即γ=1.5)、1500 km宽(x=1.5)、组分黏性比ηc=10，但若克拉通岩石圈的密度只比周围地幔的密度低0.2%(即B=0.05或Δρtc=0.132)，对流减薄时间略微减小为1.07 Ga.这说明组分密度的变化影响不大.但如仍假定300 km厚(即γ=1.5)、1500 km宽(x=1.5)、克拉通岩石圈的密度比周围地幔的密度低0.4%(即B=0.1或Δρtc=0.082)，而组分黏性比ηc=1，则依据(24)式，对流减薄时间约为426 Ma.这说明组分黏性的变化具有很大影响.
同样也可估算岩石圈厚度变化的影响.如假定240 km(即γ=1.0)厚的克拉通岩石圈，在对流的作用下减薄到初始平衡厚度120 km，根据(24)式，在克拉通岩石圈宽度为1500 km(x=1.5)和克拉通岩石圈的密度比周围地幔的密度低0.4%(即B=0.1或Δρtc=0.082)时，如组分黏性比ηc=10，克拉通岩石圈对流减薄时间大约是1.03 Ga.这显示，克拉通岩石圈的厚度也具有一定影响.
早期的研究显示，克拉通岩石圈根仅仅只有化学浮力并不能保持长期稳定，其根的强度(即克拉通岩石圈根的黏性，本文指组分黏性比ηc)是其保持长期稳定的关键因素(Lenardic and Moresi，1999).最近，Wang等(2014)在研究克拉通的稳定性时，也指出在组分黏性比ηc=3时，具有浮力的克拉通岩石圈根就能在对流的地幔中保持长达亿年不受侵蚀；而当组分黏性比ηc增大到10，即使克拉通岩石圈具有的浮力更小抑或没有本身固有的浮力，都能在地幔中长期的稳定下来.这些都与式(24)的结果相一致.
根据目前的研究结果，华北克拉通岩石圈减薄发生在晚中生代，在早白垩世(120—130 Ma)达到高潮(吴福元等，2008)，并且岩石圈的厚度从200多公里减薄到80～120 km左右，说明华北克拉通的减薄持续时间约为140 Ma.为能应用(24)式，这里假定华北克拉通从约200 km(即γ=0.67)厚经过对流减薄后减薄到120 km，空间尺度为1500 km(即x=1.5).如果根据一般克拉通都具有高黏性、低密度的特点，假定华北克拉通岩石圈的组分黏性比ηc=10(华北克拉通岩石圈的黏性是周围地幔黏性的10倍)，密度比周围地幔的密度低0.4%(即B=0.1或Δρtc=0.082)，那么根据(24)式，对流减薄时间约897 Ma.这个时间明显比140 Ma长.华北克拉通要在140 Ma得到减薄，必须发生交代置换或水化等作用使得密度或黏性产生变化.若仅发生密度变化，如交代作用使得克拉通岩石圈密度和周围地幔的密度相同(即B=0或Δρtc=0.182)，则对流减薄时间约为760 Ma；即使外推到克拉通岩石圈的密度比周围地幔高0.4%(即B=-0.1或Δρtc=0.282)，对流减薄时间也需约692 Ma.这说明如果化学交代等作用仅使克拉通岩石圈密度增大，就不足于在约140 Ma的时间内使得华北克拉通岩石圈的厚度减薄到120 km.但若考虑交代、熔融或水化等作用可以减小克拉通岩石圈的组分黏性比时，华北克拉通就可以在200 Ma内得到减薄.根据(24)式，即使假定华北克拉通岩石圈的密度仍比周围地幔的密度低0.4%，当组分黏性比减小到ηc=3时，对流减薄时间减小到约446 Ma；当组分黏性比减小到ηc=1，对流减薄时间减小到约227 Ma；如果组分黏性比降到ηc=0.5，那么克拉通岩石圈对流减薄时间就只需145 Ma.这与目前研究得到的华北克拉通可能由于太平洋的俯冲，经过脱水作用，华北克拉通岩石圈经历了强烈的水化，而大大降低了岩石圈的强度的结果一致(Xia et al.，2013;朱日祥等，2012).
4 结论本文采用二维有限元数值模拟的方法，通过增厚(γ)一定宽度(x)内的拟稳状态下的模型边界层厚度来模拟加厚的岩石圈，然后考虑由化学组分引起的岩石圈的密度(B或Δρtc)和黏性变化(ηc)，定量研究了增厚的克拉通岩石圈在对流地幔中的移除时间，结果显示：
(1)当B=0和ηc=1，即忽略克拉通岩石圈固有的密度和黏性差异时，加厚岩石圈的移除时间基本上是随着增厚岩石圈的宽度(x)增大而增大，同时也随着增厚岩石圈的厚度(γ)增大而增大，并且通过计算分析，得到移除时间和增厚岩石圈的宽度以及厚度的量化关系t2D=0.0073γ0.70x0.26.这显示，将一个初始厚度为120 km的岩石圈增厚到300 km，根据增厚的宽度(300 km或1500 km)的不同，其岩石圈根的对流移除时间在225 Ma或342 Ma不等.
(2)对克拉通岩石圈，即顾及组分密度(B或Δρtc)和组分黏性(ηc)的变化，其移除时间总体上随着增厚岩石圈的宽度(x)、厚度(γ)、密度变化(B或Δρtc)以及黏性变化(ηc)的增大而增大.当克拉通岩石圈的B以及ηc增大到一定量时，增厚岩石圈根的组分向两边蔓延并长期保持稳定；而当加厚克拉通岩石圈根的B以及ηc较小时，通过数据计算分析，得到量化关系t2Dc=0.0057ηc0.52Δρtc-0.21γ0.78ηc-0.36x0.04.这显示克拉通岩石圈的组分密度变化和宽度范围对对流减薄时间影响较小，而组分黏性比和其厚度，特别是组分黏性比具有重大的影响.对流移除时间随着组分密度的增大而增大，但增大的幅度很小，只在100 Ma量级；但随着组分黏性比的增大，对流移除时间急速地增大，能从几百Ma增大到几个Ga.
(3)对华北克拉通，如果其200 km厚的岩石圈是在140 Ma减薄到120 km厚，那么需要交代、熔融或水化作用减小其黏性.
致谢 作者感谢李曙光老师以及Shijie Zhong教授对本研究提出的宝贵建议.| [1] | Boyd F R, Gurney J J, Richardson S H. 1985. Evidence for a 150-200 km thick Archaean lithosphere from diamond inclusion thermobarometry. Nature, 315(6018): 387-389. |
| [2] | Buck W R, Toksöz M N. 1983. Thermal effects of continental collisions: thickening a variable viscosity lithosphere. Tectonophysics, 100(1-3): 53-69. |
| [3] | Carlson R W, Irving A J, Schulze D J, et al. 2004. Timing of Precambrian melt depletion and Phanerozoic refertilization events in the lithospheric mantle of the Wyoming craton and adjacent Central Plains Orogen. Lithos, 77(1-4): 453-472. |
| [4] | Carlson R W, Pearson D G, James D E. 2005. Physical, chemical, and chronological characteristics of continental mantle. Rev. Geophys., 43(1), doi: 10.1029/2004RG000156. |
| [5] | Davaille A, Jaupart C. 1994. Onset of thermal convection in fluids with temperature-dependent viscosity: application to the oceanic mantle. Journal of Geophysical Research, 99(B10): 19853-19866. |
| [6] | Davies G F. 1994. Thermomechanical erosion of the lithosphere by mantle plumes. Journal of Geophysical Research, 99(B8): 15709-15722. |
| [7] | Deng J F, Su S G, Zhao H L, et al. 2003. Deep processes of Mesozoic Yanshanian lithosphere thinning in North China. Earth Science Frontiers (in Chinese), 10(3): 41-50. |
| [8] | Duggen S, Hoernle K, Van Den Bogaard P, et al. 2003. Deep roots of the Messinian salinity crisis. Nature, 422(6932): 602-606. |
| [9] | Fan W M, Menzies M A. 1992. Destruction of aged lower lithosphere and accretion of asthenosphere mantle beneath eastern China. Geotectonica et Metallogenia (in Chinese), 16: 171-180. |
| [10] | Gao S, Rudnick R L, Carlson R W, et al. 2002. Re-Os evidence for replacement of ancient mantle lithosphere beneath the North China craton. Earth and Planetary Science Letters, 198(3-4): 307-322. |
| [11] | Grasset O, Parmentier E M. 1998. Thermal convection in a volumetrically heated, infinite Prandtl number fluid with strongly temperature-dependent viscosity: Implications for planetary thermal evolution. Journal of Geophysical Research, 103(B8): 18171-18181. |
| [12] | Griffin W L, Andi Z, O'reilly S Y, et al. 1998. Phanerozoic evolution of the lithosphere beneath the Sino-Korean Craton.// Mantle Dynamics and Plate Interactions in East Asia. American Geophysical Union, 27: 107-126. |
| [13] | Hager B H, Richards M A. 1989. Long-wavelength variations in earth's geoid-physical models and dynamical implications. Philos. Trans. Roy. Soc. A, 328(1599): 309-327. |
| [14] | Houseman G A, Mckenzie D P, Molnar P. 1981. Convective instability of a thickened boundary layer and its relevance for the thermal evolution of continental convergent belts. Journal of Geophysical Research, 86(B7): 6115-6132. |
| [15] | Huang J S, Zhong S J, Van Hunen J. 2003. Controls on sublithospheric small-scale convection. Journal of Geophysical Research, 108(B8), doi: 10.1029/2003JB002456. |
| [16] | King S D. 2005. Archean cratons and mantle dynamics. Earth and Planetary Science Letters, 234(1-2): 1-14. |
| [17] | Lenardic A, Moresi L N. 1999. Some thoughts on the stability of cratonic lithosphere: Effects of buoyancy and viscosity. Journal of Geophysical Research, 104(B6): 12747-12758. |
| [18] | Lenardic A, Moresi L N, Muhlhaus H. 2003. Longevity and stability of cratonic lithosphere: Insights from numerical simulations of coupled mantle convection and continental tectonics. Journal of Geophysical Research, 108(B6), doi: 10.1029/2002JB001859. |
| [19] | Li J H, Niu X L, Kusky T, et al. 2004. Neoarchean plate tectonic evolution of North China and its correlation with global cratonic blocks. Earth Science Frontiers (in Chinese), 11(3): 273-283. |
| [20] | Liu D Y, Nutman A P, Compston W, et al. 1992. Remnants of ≥3800 Ma crust in the Chinese Part of the Sino-Korean craton. Geology, 20(4): 339-342. |
| [21] | Lu F X, Han Z G, Zheng J P, et al. 1991. Characteristics of Paleozoic mantle-lithosphere in Fuxian, Liaoning Province. Geological Science and Technology Information (in Chinese), 10(S1): 2-20. |
| [22] | Marotta A M, Fernàndez M, Sabadini R. 1999. The onset of extension during lithospheric shortening: a two-dimensional thermomechanical model for lithospheric unrooting. Geophys. J. Int., 139(1): 98-114. |
| [23] | Menzies M A, Fan W M, Zhang M. 1993. Palaeozoic and Cenozoic lithoprobes and the loss of >120 km of Archaean lithosphere, Sino-Korean craton, China. Geological Society, London, Special Publications, 76(1): 71-81. |
| [24] | Morency C, Doin M P, Dumoulin C. 2002. Convective destabilization of a thickened continental lithosphere. Earth and Planetary Science Letters, 202(2): 303-320. |
| [25] | Moresi L N, Solomatov V S. 1995. Numerical investigation of 2D convection with extremely large viscosity variations. Phys. Fluids, 7(9): 2154-2162. |
| [26] | Olson P, Schubert G, Anderson C, et al. 1988. Plume formation and lithosphere erosion: A comparison of laboratory and numerical experiments. Journal of Geophysical Research, 93(B12): 15065-15084. |
| [27] | O'Neill C J, Lenardic A, Griffin W L, et al. 2008. Dynamics of cratons in an evolving mantle. Lithos, 102(1-2): 12-24. |
| [28] | Pearson D G, Snyder G A, Shirey S B, et al. 1995. Archean Re-Os age for Siberian eclogites and constraints on Archaean tectonics. Nature, 374(6524): 711-713. |
| [29] | Pearson D G. 1999. The age of continental roots. Lithos, 48(1-4): 171-194. |
| [30] | Qiao Y C, Guo Z Q, Shi Y L. 2012. Numerical simulation of a possible thinning mechanism of the North China Craton-Gravitational instability delamination induced by lower crust eclogite. Chinese J. Geophys. (in Chinese), 55(12): 4249-4256, doi: 10.6038/j.issn.0001-5733.2012.12.036. |
| [31] | Qiao Y C, Guo Z Q, Shi Y L. 2013. Thermal convection thinning of the North China Craton: Numerical simulation. Science China: Earth Sciences, 56(5): 773-782. |
| [32] | Rudnick R L, Fountain D M. 1995. Nature and composition of the continental crust: a Lower crustal perspective. Rev. Geophys., 33(3): 267-309. |
| [33] | Schott B, Schmeling H. 1998. Delamination and detachment of a lithospheric root. Tectonophysics, 296(3-4): 225-247. |
| [34] | Schubert G, Turcotte D L, Olson P. 2001. Mantle Convection in the Earth and Planets. Cambridge, UK: Cambridge University Press, 940. |
| [35] | Schurr B, Rietbrock A, Asch G, et al. 2006. Evidence for lithospheric detachment in the central Andes from local earthquake tomography. Tectonophysics, 415(1-4): 203-223. |
| [36] | Shao J A, Zhang L Q, Li D M. 2002. Three Proterozoic extensional events in North China Craton. Acta Petrologica Sinica (in Chinese), 18(2): 152-160. |
| [37] | Shapiro S S, Hager B H, Jordan T H. 1999. Stability and dynamics of the continental tectosphere. Lithos, 48(1-4): 115-133. |
| [38] | Solomatov V S, Moresi L N. 2000. Scaling of time-dependent stagnant lid convection: Application to small-scale convection on Earth and other terrestrial planets. Journal of Geophysical Research, 105(B9): 21795-21817. |
| [39] | Turcotte D L, Schubert G. 2002. Geodynamics. Cambridge: Cambridge University Press. |
| [40] | Wang H L, Van Hunen J, Pearson D G, et al. 2014. Craton stability and longevity: The roles of composition-dependent rheology and buoyancy. Earth and Planetary Science Letters, 391: 224-233. |
| [41] | Wu F Y, Ge W C, Sun D Y, et al. 2003. Discussions on the lithospheric thinning in eastern China. Earth Science Frontiers (in Chinese), 10(3): 51-60. |
| [42] | Wu F Y, Lin J Q, Wilde S A, et al. 2005. Nature and significance of the Early Cretaceous giant igneous event in eastern China. Earth and Planetary Science Letters, 233(1-2): 103-119. |
| [43] | Wu F Y, Xu Y G, Gao S, et al. 2008. Lithospheric thinning and destruction of the North China Craton. Acta Petrologica Sinica (in Chinese), 24(6): 1145-1174. |
| [44] | Xia Q K, Liu J, Liu S C, et al. 2013. High water content in Mesozoic primitive basalts of the North China Craton and implications on the destruction of cratonic mantle lithosphere. Earth and Planetary Science Letters, 361: 85-97. |
| [45] | Xu W L, Wang Q H, Wang D Y, et al. 2004. Processes and mechanism of Mesozoic lithospheric thinning in eastern North China Craton: evidence from Mesozoic igneous rocks and deep-seated xenoliths. Earth Science Frontiers (in Chinese), 11(3): 309-313. |
| [46] | Xu Y G. 2001. Thermo-tectonic destruction of the Archaean lithospheric keel beneath the Sino-Korean Craton in China: Evidence, timing and mechanism. Phys. Chem. Earth A, 26(9-10): 747-757. |
| [47] | Yoshida M. 2012. Dynamic role of the rheological contrast between cratonic and oceanic lithospheres in the longevity of cratonic lithosphere: A three-dimensional numerical study. Tectonophysics, 532-535: 156-166. |
| [48] | Yuen D A, Fleitout L. 1985. Thinning of the lithosphere by small-scale convective destabilization. Nature, 313(5998): 125-128. |
| [49] | Zhao G C, Wilde S A, Cawood P A, et al. 2001. Archean blocks and their boundaries in the North China Craton: lithological, geochemical, structural and P-T path constraints and tectonic evolution. Precambrian Res., 107(1-2): 45-73. |
| [50] | Zhao L, Zheng T Y. 2007. Complex upper-mantle deformation beneath the North China Craton: implications for lithospheric thinning. Geophys. J. Int., 170(3): 1095-1099. |
| [51] | Zhao Y, Xu G, Zhang S Z, et al. 2004. Yanshanian movement and conversion of tectonic regimes in East Asia. Earth Science Frontiers (in Chinese), 11(3): 319-328. |
| [52] | Zhong S J, Zuber M T, Moresi L, et al. 2000. Role of temperature-dependent viscosity and surface plates in spherical shell models of mantle convection. Journal of Geophysical Research, 105(B5): 11063-11082. |
| [53] | Zhu R X, Chen L, Wu F Y, et al. 2011. Timing, scale and mechanism of the destruction of the North China Craton. Science China: Earth Sciences, 54(6): 789-797. |
| [54] | Zhu R X, Xu Y G, Zhu G, et al. 2012. Destruction of the North China Craton. Science China: Earth Sciences, 55(10): 1565-1587. |
| [55] | 邓晋福, 苏尚国, 赵海玲等. 2003. 华北地区燕山期岩石圈减薄的深部过程. 地学前缘, 10(3): 41-50. |
| [56] | 范蔚茗, Menzies M A. 1992. 中国东部古老岩石圈下部的破坏和软流圈地幔的增生. 大地构造与成矿学, 16: 171-180. |
| [57] | 李江海, 牛向龙, Kusky T等. 2004. 从全球对比探讨华北克拉通早期地质演化与板块构造过程. 地学前缘, 11(3): 273-283. |
| [58] | 路凤香, 韩柱国, 郑建平等. 1991. 辽宁复县地区古生代岩石圈地幔特征. 地质科技情报, 10(S1): 2-20. |
| [59] | 乔彦超, 郭子祺, 石耀霖. 2012. 数值模拟华北克拉通岩石圈减薄的一种可能机制——下地壳榴辉岩重力失稳引起的拆沉. 地球物理学报, 55(12): 4249-4256, doi: 10.6038/j.issn.0001-5733. 2012.12.036. |
| [60] | 乔彦超, 郭子祺, 石耀霖. 2013. 数值模拟华北克拉通岩石圈热对流侵蚀减薄机制. 中国科学: 地球科学, 43(4): 642-652. |
| [61] | 邵济安, 张履桥, 李大明. 2002. 华北克拉通元古代的三次伸展事件. 岩石学报, 18(2): 152-160. |
| [62] | 吴福元, 葛文春, 孙德有等. 2003. 中国东部岩石圈减薄研究中的几个问题. 地学前缘, 10(3): 51-60. |
| [63] | 吴福元, 徐义刚, 高山等. 2008. 华北岩石圈减薄与克拉通破坏研究的主要学术争论. 岩石学报, 24(6): 1145-1174. |
| [64] | 许文良, 王清海, 王冬艳等. 2004. 华北克拉通东部中生代岩石圈减薄的过程与机制: 中生代火成岩和深源捕虏体证据. 地学前缘, 11(3): 309-313. |
| [65] | 赵越, 徐刚, 张拴宏等. 2004. 燕山运动与东亚构造体制的转变. 地学前缘, 11(3): 319-328. |
| [66] | 朱日祥, 陈凌, 吴福元等. 2011. 华北克拉通破坏的时间、范围与机制. 中国科学: 地球科学, 41(5): 583-592. |
| [67] | 朱日祥, 徐义刚, 朱光等. 2012. 华北克拉通破坏. 中国科学: 地球科学, 42(8): 1135-1159. |