大陆岩石圈、地幔底部异常体与地幔对流相互作用的数值模拟

杨亭，傅容珊，黄川，班磊

蒙城地球物理国家野外科学观测研究站，中国科学技术大学地球和空间科学学院，合肥 230026

1 引言

在地表，存在着占其面积约30%的大陆岩石圈.不同于作为地幔对流顶部热边界层的大洋岩石圈 （McKenzie etal.，2005），大陆岩石圈，尤其是克拉通，具有较小的固有密度（intrinsic density），较低的水含量，较高的黏度等物理化学性质 （Artemieva，2009；Karato，2010；Khan etal.，2011；Peslier etal.，2010）.由于特殊的物理化学性质，大陆岩石圈一般不主动参与下方的地幔对流，但和地幔对流有着重要的相互影响.

为了探讨大陆岩石圈与周围地幔对流结构的相互影响，前人进行了一系列的数值实验.Gurnis（1988）第一次给出了大陆与地幔之间相互作用的数值模型.其研究发现，大陆在冷的下涌流区域聚合，形成超级大陆.由于超级大陆的形成，阻止了原下涌区的俯冲和地幔冷却.超级大陆下方的地幔逐渐增温，形成上涌流，由于拉张作用，将超级大陆分裂成小块大陆.这一过程周期性地进行.Phillips和Bunge（2005）在3D球域讨论了大陆的大小，地幔加热模式，下地幔黏度增加对于大陆运动的影响.其结果表明，大陆岩石圈的漂移速度为海洋岩石圈速度的1／3左右，这解释了大陆岩石圈和海洋岩石圈间速度的差异.Zhong等（2007）在3D球域探讨了大陆与地幔之间的相互作用并根据数值模型认为，当地幔具有中等强度的岩石圈，较弱的上地幔和高黏度的下地幔时，地幔在没有超级大陆时呈一阶对流模式；在有超级大陆时呈二阶对流模式.由于大陆对地幔对流的这种调节作用，使得地幔在两种对流模式间周期性地相互转换，进而导致超级大陆周期性地拼合与分裂.Coltice等（2007）研究了大陆岩石圈对地幔温度的影响.其结果表明，超级大陆的形成会导致地幔对流波长的增大，进而导致岩石圈下地幔温度高达100℃的增加，该现象被用来解释中大西洋玄武岩省.作者进而认为至少存在两种大陆玄武岩省：由热柱产生的和地幔全球变暖导致的.Yoshida（2010）考虑了大陆岩石圈组分和流变的横向不均匀性对其稳定性的影响.其结果表明，弱的大陆边缘能够保护克拉通岩石圈不被对流地幔所拉张，因此在岩石圈的稳定性方面具有重要作用.

上述研究为我们理解大陆岩石圈和周围地幔之间的相互作用提供了许多新的视角.随着观测技术的提高，人们逐渐认识到，在下地幔底部，存在着两块具有较高密度的热化学异常体（下称LSVPs）.这两块热化学异常体分别位于中太平洋和非洲底部（He and Wen，2009；Ishii and Tromp，1999；Kennett etal.，1998；Ni etal.，2002），共占据了CMB面积的约20%（Burke etal.，2008）.LSVPs受到周围地幔的影响而变形、运动（Zhang etal.，2010），同时其对地幔对流的结构也有重要影响（Davaille，1999）.

前人进行了一系列的数值和实验室实验来探讨LSVPs和地幔对流结构之间的相互影响.如Davaille（1999）在实验室中进行了一系列分层密度和黏度的热化学对流实验.其结果给我们两点启示：（1）即便是1%～2%的化学密度不均匀，也会对动力学产生显著的影响；（2）现在的地幔对流格局并不能够代表整个地球历史的地幔对流格局.特别地，作者认为在地球历史早期，地幔对流呈分层机制；但是，随着两个物质层的逐渐混合，其密度差异逐渐变小，导致地幔对流从分层机制转变为震荡的doming形式.作者认为dome和在dome顶端形成的小热柱是现今观测到的太平洋和非洲超级热柱和上方热点的成因.McNamara和Zhong（2005）以过去1.2亿年的板块运动历史作为上表面的速度边界条件，探讨了LSVPs物质形态的演化.其得到的LSVPs结构与地震观测的非洲和太平洋的热化学异常结构具有相似的形态.Lassak等（2010，2007）研究了LSVPs对CMB形态的影响.其结果显示，纯热起源的超级热柱（热柱丛）和热化学起源的超级热柱引起的CMB地形是不同的：对于热化学起源的超级热柱，其CMB地形在LSVPs下方相对平坦，且稍微升高，而在LSVPs的边缘，却剧烈上升，CMB地形呈平底锅形状；而热柱丛模型预测的CMB地形在每个热柱下面为高地形.Tan等（2011）在3D球域可压缩地幔对流条件下考察了热柱的位置和高密度高体积模量的LSVPs的水平运动.其研究表明，选取LSVPs的体积模量和密度在合理的参数范围内，能够保证dome结构存在几十亿年而不被地幔对流携带走.同时，热化学热柱更多地倾向于在dome的边缘而非顶部产生.这和Burke等（2008）的观测一致.

实际的地幔，同时包含大陆岩石圈，周围地幔和LSVPs.因此，在考虑大陆岩石圈和地幔对流的相互影响时，应该将LSVPs同时包含进来.但是，迄今为止，同时包含上述三种成分的地幔对流还未深入展开.本文修改了计算地幔对流的有限元代码Citcom（McNamara and Zhong，2004；Moresi and Solomatov，1995），以使得能够计算包含三种组分（大陆岩石圈，周围地幔，LSVPs）的地幔对流.并探讨大陆岩石圈和周围地幔及LSVPs的相互影响.

2 模型和方法

2.1 控制方程

本文假定地幔为满足Boussinesq近似和无限大普朗特数假设的不可压缩流体.在2D直角坐标下讨论地幔的演化.满足上述假设的无量纲的质量、动量和能量守恒方程为：

另外，对于热化学对流，化学组分的运动由平流方程所控制：

上述方程中，ui为速度，P为动力学压强，μ为黏度，T为温度，C为化学组分，gi＝ （gx，gz）＝ （0，1）为重力加速度，X，i为X的空间导数，X·为X的时间导数，Ra，Rc分别为底部加热和化学瑞利数.由下述公式给出：

公式右侧的变量是有量纲的，其含义及数值在表1中给出.

Rc与Ra之比为浮力数B：

浮力数表征化学浮力与热浮力的相对大小.

表1 文中物理参数的数值Table 1 Values of physical parameters in this article

2.2 模型参数

为了初步探讨大陆，LSVPs与周围地幔间的相互作用，本文在2D直角坐标下设置了两类模型：包含大陆和周围地幔两种物质的模型；同时包含大陆，周围地幔和LSVPs三种物质的模型.其中，大陆岩石圈的平均密度比周围地幔低2.9%，LSVPs的平均密度比周围地幔高2.9%（图1）.

地幔的黏度随着温度、组分、压力、应变率、矿物相等复杂地变化.上述因素对地幔黏度的影响究竟如何，至今仍存在很大不确定性（Tackley，2011）.作为初步研究，本文简单地将无量纲的黏度设置为随着深度和组分而变化.

其中，η0是依赖于深度的前置因子，对岩石圈（0～100km），软流圈（100～200km），过渡区（200～660km），下地幔（660～2890km）分别设置为1，0.01，0.1，1.即拥有高黏度的岩石圈和下地幔及低黏度的软流圈.黏滞比ηc表示在同一深度下，化学异常物质（包括大陆岩石圈和LSVPs）和周围地幔的黏度之比.对于大陆岩石圈，将这一值设置为100.LSVPs由于密度较高，导致其停留在CMB上方而不参与上层地幔的对流，这使其温度较周围地幔高.根据Arrhenius公式，其黏度应比周围地幔低.但另一方面，Solomatov等（1996）认为在下地幔，活化能可能取负值，因此，高温的LSVPs可能黏度更高（Korenaga，2005；Solomatov，1996）.钙钛矿到过钙钛矿的相变具有正的克拉贝隆斜率（Murakami etal.，2004；Oganov and Ono，2004），在下涌区，由于温度低，地幔物质更可能经相变转化为过钙钛矿相.因此，低温的下涌区具有较厚的过钙钛矿相（Nakagawa and Tackley，2005）.而过钙钛矿相的黏度比钙钛矿小3个量级（Ammann etal.，2010），这导致高温的LSVPs可能具有较高的黏度.鉴于此，我们分别探讨了LSVPs的ηc取为25或0.04对于模型结果的影响（图1）.

本文使用粒子比方法来模拟化学成分随着时间的演化（Tackley and King，2003）.初始时刻，每一单元拥有16个追踪粒子.追踪粒子属于以下三种类型中的一个：-1，0，1，分别用以表示大陆岩石圈，周围地幔，LSVPs三种物质.随着地幔的演化，追踪粒子会随着地幔流扩散.某一时刻，每一单元的组分由下面公式所确定：

式中，Ce，NC，NM，NL分别为某一时刻，某一单元的组分，以及单元中大陆岩石圈类型粒子，周围地幔类型粒子和LSVPs类型粒子的数目.从单元组分公式可知，单元组分Ce在-1～1之间变化.当某一单元全部为-1类型粒子时，单元组分为-1，此时该单元属于大陆岩石圈；当某一单元中全部为1类型粒子时，单元组分为1，此时该单元属于LSVPs；当某一单元中全部为0类型粒子时，单元组分为0，此时该单元属于周围地幔.在结果分析时，本文将组分场小于-0.8的区域判定为大陆岩石圈，组分场大于0.8的区域判定为LSVPs区域，组分场在-0.8～0.8之间的区域判定为周围地幔.

初始时刻大陆岩石圈设为厚度0.1（有量纲厚度为289km），位于地表中央的长方形块体（图1）.为了探讨大陆岩石圈的尺度对地幔对流的影响，我们将初始时刻大陆岩石圈的大小设置为不同的值（表2）.考虑到高黏度的LSVPs可能保持一定的块状结构而不因重力作用塌陷铺满CMB，因此，对于高黏度的LSVPs，我们将其初始时刻设置为一个方形块体（图1a，1b）.块体的大小等于现今非洲和太平洋块体热化学异常的总和.对于低黏度的LSVPs，我们将其在初始时刻设置为厚度为0.0865（250km）的平层（图1c，1d）.模型的初始温度场由具有相同参数的纯热对流计算得到的稳态温度场的径向平均温度加上随机温度扰动得到.模型各个边界均为自由滑移速度边界条件.顶部和底部采用等温边界条件，侧边界采用绝热边界条件.

由于真实的地幔对流没有侧向边界，在使用2D模型进行仿真时，模型的侧边界可能对对流结构具有影响.为了减小侧边界的影响，本文将模型的长宽比设为6.模型将计算区域剖分为128×384（垂向单元数×水平向单元数）个长方形单元，每一单元的边长为22.6km（垂向）×45.2km（水平）.

3 结果及分析

为了便于讨论，我们根据地表的组分场将计算区域划分为三个部分（见图2）：大陆区域（离开洋陆边界500km以上的地表组分场为大陆物质的区域），大陆边缘（距离洋陆边界500km以内的区域），海洋区域（离开洋陆边界500km以上的地表组分场为海洋物质的区域）.其中陆洋边界指地表组分场为0.5的等组分线.之后我们分别讨论了三个区域岩石圈及下伏地幔各项物理参数随时间的演化及在稳态时的平均值（表2）.这些参数包括：

vzCB：大陆边缘地幔平均垂向速度，单位cm／a，通过计算大陆边缘区域从岩石圈底部到CMB地幔垂向速度的平均值得到.

Δvz：陆洋区域地幔平均垂向速度之差，单位cm／a.

ΔTU：陆洋上地幔平均温差，单位K，其中，上地幔温度的取值深度范围为200～660km.

ΔTL：陆洋下地幔底部平均温度之差，单位K，下地幔底部温度的取值深度范围为1890～2890km，即地幔最底部的1000km.

ΔTLAB：表示陆洋LAB平均温度之差，单位K；

图1 模型CaseB1（a，b）和模型CaseC1（c，d）的无量纲的初始密度场（a，c）和黏度场（b，d）（a），（c）中颜色表示初始时刻无量纲密度.（b），（d）中颜色表示初始时刻无量纲黏度的常用对数.横坐标表示无量纲长度，纵坐标表示无量纲高度.每1单位代表2890km.Fig.1 Non-dimensional density（a，c）and viscosity（b，d）distribution at initial time for model CaseB1（a，b）and CaseC1（c，d）X-axis stands for non-dimensional length and Y-axis stands for non-dimensional height.The scalar for length and depth is 2890km.

图2 大陆区域、大陆边缘和海洋区域的分布示意图颜色表示组分.-1表示大陆岩石圈物质，0表示周围地幔物质.O代表海洋区域，CB代表大陆边缘区域，C代表大陆区域.横轴表示无量纲长度，纵轴表示无量纲深度.每1单位代表2890km.Fig.2 Schematic view of the distribution of the Continental region，Oceanic region and Continental margin region Color represent the composition.-1stands for continental lithosphere material and 0stands for surrounding mantle material.O stands for Oceanic region，CB stands for Continental margin region and C stands for Continental region.X-axis stands for non-dimensional length and Y-axis stands for nondimensional depth.The scalar for length and depth is 2890km.

ΔdLAB：陆洋LAB平均深度之差，单位km.其中，某一时刻大陆区域的平均LAB深度及温度由下述方法计算：将该时刻所有位于大陆区域的温度剖面进行平均，得到平均的温度剖面（如图3右），从上往下，取该温度剖面的第一个局部极大值作为大陆区域平均的LAB温度，其对应的深度为大陆区域平均的LAB深度.海洋区域平均LAB深度及温度的计算与大陆区域相同.

Rvx：洋陆岩石圈平均水平速率之比.

Rca：陆洋区域LSVPs物质平均丰度比.其中，LSVPs物质的选取范围为1890～2890km，即最底部的1000km.

3.1 包含大陆岩石圈和周围地幔的地幔对流

CaseA1，CaseA2，CaseA3三组模型只包含大陆岩石圈和周围地幔.三组模型的差别在于初始时刻大陆的无量纲长度分别为1，3和0.5（见表2）.

当大陆面积较小时，大陆边缘往往伴随着剧烈的俯冲（vzCB＜0，见表2，CaseA3）.由于大陆面积较小，大陆下方的地幔也受边缘俯冲的影响，而产生下涌.因此，与海洋区域相比，大陆下方的地幔总体向下流动（Δvz＜0）.由于质量守恒，在海洋区域，地幔总体上处于上涌状态.相对应的，海洋区域的地幔温度大于大陆区域的地幔温度（ΔTU＜0，ΔTL＜0）.

表2 本文计算的模型及模型结果.L0代表无量纲的初始大陆长度，ηc代表LSVPs与周围地幔黏度比，NaN表示该值没有意义.其他符号意义及单位见第3节Table 2 Simulated models in this article and model results.L0represents the initial length of the continental lithosphere.ηcrepresents the viscosity ratio between the LSVPs and the surrounding mantle.NaNrepresents Not a Number.Check other parameters in Section 3

图3 CaseA1模型某一时刻大陆（虚线）与海洋（实线）区域的平均垂向速度（左）和温度剖面（右）左图中，速度向上为正，向下为负.右图中实短横线和虚短横线分别代表海洋与大陆的平均LAB.Fig.3 Average vertical velocity（left）and temperature（right）profile for model CaseA1 at certain time after reaching steady state Positive is for upward velocity while negative is for downward velocity in the left panel.Horizontal ticks in the right panel represent LAB.Dashed（solid）line is for continent（Oceanic）region.

随着大陆面积的增大，大陆边缘俯冲逐渐减弱，大陆区域地幔的下涌速度逐渐减小，温度逐渐升高.对于超级大陆的情况（CaseA2），大陆边缘的下涌非常弱（vzCB幅值较小），大陆区域的垂向速度甚至稍大于海洋区域的垂向速度（Δvz＞0）.这表明，在此时，大陆内部以上涌为主（Gurnis，1988）.对应的，超级大陆下方地幔的温度也较高（ΔTU＞0）.在超级大陆下方的上涌流和高温，可能为超级大陆的裂解、伴随超级大陆裂解产生的大量的高温火山喷发（Condie，2004）提供了物质和能量上的支持.

随着大陆面积的增加，大陆与海洋岩石圈软流圈界面（LAB）深度和温度的差异逐渐减小.但是，在我们的模型中，大陆LAB的深度总是大于海洋区域的LAB深度.这表明大陆岩石圈的热学厚度总是大于海洋岩石圈的热学厚度.对于CaseA2的超级大陆模型，大陆区域上地幔的温度比海洋区域上地幔的温度高19K（ΔTU＝19），但大陆区域LAB的温度仅比海洋下LAB的温度高1K（ΔTLAB＝1）.这表明，虽然超级大陆下方上地幔的温度比海洋区域高，但是大陆岩石圈本身的温度却和海洋岩石圈温度相当.由于大陆岩石圈中的craton因亏损水而具有较高的黏度（Peslier etal.，2010）.因此，相对于海洋岩石圈，大陆的craton部分在超级大陆时期可能也具有较高的强度.

随着大陆面积的增加，大陆内部的水平应力从压应力为主转变为张应力为主（图4）.本文模型中，大陆岩石圈黏度没有横向变化，因此，大陆显示出较为均匀的拉张或压缩（图4）.在真实的地球中，超级大陆的拉张（分裂）可能更倾向于位于大陆的薄弱地区.超级大陆内部的拉张应力，可能促进了超级大陆的裂解（Yoshida，2010）.

大陆岩石圈的水平速度，始终小于海洋岩石圈的水平速度（Rvx＞1）.这和观测一致.

我们以CaseA1模型为例探讨大陆及海洋区域各项物理参数随时间的演化及其相互关系.首先，大陆区域地幔的平均垂向运动速度，海洋区域地幔的平均垂向运动速度，海洋区域地表热流和大陆区域上地幔温度间存在很好的关联性.如大陆区域地幔的平均垂向运动速度和海洋区域地幔的平均垂向运动速度，大陆区域上地幔温度反相关（图5a和图5e）.大陆区域地幔的下涌速度大的时候，海洋区域地幔的上涌速度大.在大陆区域地幔下涌流达到最大值的一段时间后，海洋区域地表热流也达到最大值（图5a和图5c）.其次，大陆区域地幔垂向运动速度和岩石圈在地表的水平运动速度，大陆区域地表热流，LAB温度，下地幔底部温度间似乎没有明显的相关性.如下地幔底部的温度变化和大陆区域地幔的下涌速度间没有明显相关性（图5a和图5f）.其原因可能在于下地幔底部的垂向运动较小（图3左），导致上方的短周期变化没有足够的时间来影响下地幔底部.值得注意的是，大陆地表热流随时间变化的频率和幅度远小于海洋地区地表热流随时间变化的频率和振幅.由于大陆地区地表热流由大陆地区的近地表地温梯度来确定（热传导），这意味着，地幔中短周期的温度、运动速度等的变化，不会对浅表地温梯度造成影响.

为了考察不同时期的地幔结构和上述物理量之间的关系，我们给出了模型CaseA1在图5中t1，t2，t3三个时刻的温度场（图6）.在t1时刻大陆区域地幔的下涌速度达到瞬时极大值（图5a）.此时，大陆区域同时具有较低的上地幔温度（图5e）.而温度场显示此时大陆的两侧都伴随着剧烈的俯冲（图6a）.t2时刻大陆区域地幔中的下涌流微弱（图5a）.温度场也显示此时大陆两侧均没有强烈的俯冲（图6b）.t3时刻大陆边缘存在着剧烈的俯冲（图5a，图6c）.但是，由于俯冲板块在下地幔并非处于大陆下方（图6c），此时大陆区域下地幔底部温度反而处于瞬时极大值（图5f）.可见，对于较小的大陆而言，某一时刻，大陆区域地幔底部的温度和大陆区域是否存在俯冲并没有必然的联系.

3.2 包含大陆岩石圈，周围地幔和LSVPs的地幔对流

CaseB1，CaseB2，CaseB3以及CaseC1，CaseC2，CaseC3六组模型同时包含大陆岩石圈和LSVPs.其中CaseB1，CaseB2，CaseB3中，LSVPs用体积相当于现今太平洋和非洲LSVPs体积之和的高黏度块体来表示（图1a，1b）.CaseC1，CaseC2，CaseC3中，LSVPs用厚度为250km的低黏度的热化学层来表示（图1c，1d）.其他模型参数参见表2.

不论LSVPs是高黏度的块体还是低黏度的层，包含LSVPs都不会改变各物理量（如大陆边缘区域地幔垂向速度）随着大陆体积的增大而变化的趋势（表2）.其特别之处在于：（1）LSVPs的出现，会减小大陆边缘的俯冲（表2）.尤其是对于非超级大陆的情况，这种减小更加强烈.（2）LSVPs物质的丰度在大陆区域与在海洋区域之比Rca随着大陆体积的增加而增加.在地球历史早期，地幔底部LSVPs物质的含量可能高于今天（Davaille，1999），但是，由于大陆岩石圈的总体体积较小（Piper，2010），当时LSVPs可能主要位于海洋区域.

图4 CaseA1，CaseA2，CaseA3模型在同一时间的水平方向正应力场拉应力为正，压应力为负.线条所围区域为大陆岩石圈.Fig.4 Normal stress in the horizontal direction shown in the same time for CaseA1，CaseA2，CaseA3，respectively Positive for tensional and negative for compressional stress，respectively.Regions encircled by the solid lines represent the continental lithosphere regions.

图5 CaseA1模型相关物理参数随时间的演化虚线表示大陆区域，固线表示海洋区域.（a）为地幔平均垂向速度，向上为正.（b）为平均地表热流.（c）为上地幔平均过绝热温度.（d）为地幔底部1000km平均过绝热温度.Fig.5 Evolution of physical quantities with time for model CaseA1 Dashed（solid）line is for continent（Oceanic）region.（a）Averaged vertical velocity in the mantle.Upward direction is positive.（b）Averaged surface heat flux.（c）Averaged superadiabatic temperature in the upper mantle.（d）Averaged superadiabatic temperature in the lowermost mantle.

图6 CaseA1模型不同时刻的温度场（a），（b），（c）分别对应着图5中t1，t2，t3三个时刻的温度场.图中向下箭头所指为大陆岩石圈的左右两个端点.Fig.6 Snapshots of the temperature field for CaseA1at different time Panels（a），（b）and（c）are the temperature snapshots corresponding to time t1，t2and t3in Fig.5，respectively.Black arrows point the boundary of the continental lithosphere.

我们以CaseC1模型为例，探讨存在LSVPs的情况下，大陆及海洋区域各项物理参数随时间的演化及其相互关系（图7）.在总体上，各项物理参数随时间的演化及其相互关系与CaseA1类似.如大陆区域的地幔平均垂向速度随时间波动剧烈而海洋区域随时间波动较小；海洋区域热流随时间波动剧烈而大陆区域随时间波动小等等.下面我们考察大陆和海洋区域LSVPs丰度随时间的演化（图7f）.t1时刻到t2时刻之间，大陆区域LSVPs的丰度一直保持在低值.地幔结构（图8a，8b）显示，此时LSVPs主要有两个，它们均位于远离大陆的海洋区域.从t2到t3时刻，大陆下方的LSVPs物质丰度迅速达到瞬时极大值，其原因在于大陆岩石圈向左运移，使得自身位于原先的LSVPs物质上方（图8b，8c）.t3到t4时刻，大陆区域LSVPs的丰度开始从先前的高值迅速向之后的低值转化.其主要原因可能在于大陆周围的下涌流导致了LSVPs沿着CMB向右滑移（图8d）.这可能意味着当大陆水平漂移时，地幔底部LSVPs的分布与地表大陆的位置无关.但当大陆位置长时间保持不变时，其通过控制下涌流的位置控制了地幔底部的LSVPs的空间分布.值得注意的是，CaseC1大陆区域LSVPs物质的丰度随时间的变化与大陆区域下地幔底部温度随时间的变化一致（图7d，7f）.当大陆下方LSVPs物质含量多时，大陆下方下地幔的平均温度高，而当大陆下方LSVPs物质含量少时，大陆下方下地幔的平均温度低.CaseC1的海洋平均热流约90mW·m-2，大陆区域的海洋平均热流约为30mW·m-2（图7c），考虑到大陆区域岩石圈内放射性内生热产生的地表热流约为37mW·m-2（Turcotte and Schubert，2002）.本文计算的海洋地表热流（90mW·m-2）与大陆地表热流（30+37＝67mW·m-2）均与观测地表热流（分别为101mW·m-2和65mW·m-2）较为一致.这表明本文选取的物理参数与地球真实的物理参数较为一致.

图7 CaseC1模型相关物理参数随时间的演化虚线表示大陆区域，固线表示海洋区域.（a）为地幔平均垂向速度，向上为正.（b）为平均地表热流.（c）表示下地幔平均过绝热温度.（d）表示下地幔底部平均组分（用来表示LSVPs的平均丰度）.Fig.7 Variation of physical quantities with time for model CaseC1 Dashed（solid）line is for continent（Oceanic）region.（a）Averaged vertical velocity in the mantle.Upward direction is positive.（b）Averaged surface heat flux.（c）Averaged superadiabatic temperature in the lowermost mantle.（d）Averaged composition in the lowermost mantle.

图8 CaseC1模型不同时刻的温度场（a），（b），（c），（d）四幅图分别对应着图7中t1，t2，t3，t4四个时刻的温度场.（e）为t1时刻的流函数图，红色代表逆时针方向运动，蓝线代表顺时针方向运动.图中，顶部粗蓝线包围区域为大陆岩石圈，底部红色粗线包围区域为LSVPs.Fig.8 Snapshots of the temperature field for CaseC1at different times Panels（a），（b），（c）and（d）correspond to time t1，t2，t3and t4in Fig.7，respectively.Panel（e）shows the streamline corresponding to time t1.Thick blue and red lines corresponds to the boundaries of continental lithosphere and LSVPs respectively.

LSVPs及地表大陆岩石圈的出现，对地幔对流的结构具有显著的影响.在下地幔底部，沿CMB的水平流动受到LSVPs的阻碍而转为向上的流动（图8e）.这可能是在LSVPs的边缘形成尖锐及陡峭界面（Ni etal.，2002）产生许多热柱（Burke etal.，2008）的原因.由于周围地幔流的影响，在LSVPs顶部形成凸凹的界面，在界面的高峰处，往往存在许多热柱（图8）.在上地幔，水平流动在大陆边缘转化为向下的流动，形成俯冲带.为了进一步比较地表的大陆岩石圈和地幔底部的LSVPs，我们考察了CaseC1模型CMB热流随时间的演化（图9）.CMB热流在LSVPs区域较低，而在其他区域较高.这反映了地幔底部LSVPs对地核的保温作用，与大陆岩石圈类似.

图9 CaseC1模型CMB热流随时间的演化固线表示LSVPs区域，虚线表示非LSVPs区域.Fig.9 CMB heat flux varying with time for CaseC1 Solid and dashed lines represent for LSVPs and non-LSVPs regions，respectively.

4 结论与讨论

在与地球相似的物理参数下，本文在2D直角坐标域计算了同时包含大陆岩石圈，LSVPs的地幔对流，以探讨大陆岩石圈与周围地幔、LSVPs间的相互影响.本文主要有以下结论：

（1）大陆边缘俯冲的强弱和大陆的体积呈反相关.当大陆体积较小时，大陆边缘常伴随着较强的俯冲.随着大陆体积的增大，大陆边缘的俯冲逐渐减弱.与之对应的，大陆区域地幔在大陆体积小时呈现下涌状态，在体积大时呈现上涌状态.大陆岩石圈在体积较小时水平方向受到压应力，而当体积大时，水平方向受到张应力.

需要指出的是，本文在地表只设置一个大陆岩石圈.由于当今地球存在多个大陆岩石圈，它们之间及它们与地幔流之间存在复杂的相互作用，上述结论可能不适用于当今地球.但是，上述结论可用来研究超级大陆和地幔间的相互影响.本文结果显示，对于Pangea（距今时间最近的超级大陆，形成于～320Ma（Li and Zhong，2009））这样的超级大陆而言，其周围俯冲较弱，内部存在较强的上涌流，大陆下方的上地幔具有较高的温度，大陆岩石圈以拉张应力为主.最早的超级大陆可能在太古代已出现（Piper，2010），当时大陆面积可能比现今大陆面积小（Condie，2000；Piper，2010）.根据本文研究，我们认为太古代时期，超级大陆区域地幔的上涌流应弱于Pangea超级大陆区域的上涌流，其岩石圈内水平张力也应小于Pangea.这可能使得Archean时期的超级大陆在更长时间内不被裂解（Piper，2010）.

（2）岩石圈与软流圈边界面（LAB）在大陆下方较深，温度较低；在海洋区域较浅，温度较高.随着大陆体积的增大，陆洋之间LAB深度、温度的差异逐渐减小.

（3）大陆岩石圈的水平漂移速度比海洋岩石圈的水平运动速度小.这和观测到的新生代以来的板块运动速度一致（Gordon and Jurdy，1986）.

（4）海洋区域地表热流高，并且随着时间的波动较大，而大陆区域的地表热流低，且随着时间的波动较小.LSVPs区域CMB热流低.

由于大陆区域地表热流随时间波动较小，故海洋区域地表热流随时间的变化，可用来近似估计全球地表热流随时间的演化.

（5）大陆下方LSVPs物质的丰度与大陆的体积呈正相关.

当大陆体积小时，大陆区域地幔经常处于下涌状态，而下地幔的LSVPs堆总是伴随着上涌流（图8e），因此，大陆与LSVPs趋向于彼此远离.这导致大陆下方LSVPs物质丰度比海洋区域小.随着大陆体积的增大，上涌区及其对应的LSVPs倾向于分布在大陆内部或海洋区域.这导致超级大陆区域LSVPs物质的丰度，和海洋区域接近.考虑到超级大陆Pangea大致以现今非洲为中心.本文认为非洲LSVPs和太平洋LSVPs可能分别源于Pangea超级大陆时期超级大陆和海洋区域下方的LSVPs.

（6）LSVPs和大陆岩石圈对于地幔对流结构具有重要的影响.在LSVPs的周围易于形成上涌流和热柱，在大陆岩石圈周围易于形成俯冲带.

本文对低黏度和高黏度LSVPs的模型结果显示，地幔底部LSVPs黏度的高低不影响上述结论.

本文采用两种方法验证包含三种成分热化学对流程序的正确性（结果未列出）.首先，本文所谓化学异常，是指大陆岩石圈、LSVPs和周围地幔间密度和黏度差异.如果我们将两者密度和黏度差异改为零，则本文热化学对流的结果应和纯热对流的计算结果一致.我们使用和Blankenbach等（1989）文中Case1C相同的参数，计算结果与其一致.其次，将大陆的体积设置为0，采用和Van Keken等（1997）包含两种物质对流模型相同的参数，本文程序计算结果与其一致.

需要指出的是，本文的结论是在特定的物理参数下得到的.如只考虑了深度和组分相关的黏度，而真实地幔的黏度比这要复杂得多.另一方面，我们的模型采用二维直角坐标系，且只包含一个大陆岩石圈，而地球本身是由多个大陆岩石圈组成的三维球体.因此，本文的数值模型只是最初步的结果，要得到更接近地球实际的结果，需要在3D计算平台下，进行包含多个大陆的更广泛的数值实验.

致谢 在程序使用方面，作者得到了钟时杰教授、黄金水教授及冷伟教授的帮助；两位匿名专家提供了宝贵的意见，在此表示感谢.

Ammann M W，Brodholt J P，Wookey J，etal.2010.Firstprinciples constraints on diffusion in lower-mantle minerals and a weak D layer.Nature，465（7297）：462-465.

Artemieva I M.2009.The continental lithosphere：reconciling thermal，seismic，and petrologic data.Lithos，109（1-2）：23-46.

Blankenbach B，Busse F，Christensen U，etal.1989.A benchmark comparison for mantle convection codes.Geophysical Journal International，98（1）：23-38.

Burke K，Steinberger B，Torsvik T H，etal.2008.Plume generation zones at the margins of large low shear velocity provinces on the core-mantle boundary.Earth and Planetary Science Letters，265（1-2）：49-60.

Coltice N，Phillips B R，Bertrand H，etal.2007.Global warming of the mantle at the origin of flood basalts over supercontinents.Geology，35（5）：391-394.

Condie K C.2000.Episodic continental growth models：afterthoughts and extensions.Tectonophysics，322（1-2）：153-162.

Condie K C.2004.Supercontinents and superplume events：distinguishing signals in the geologic record.Physics of the Earth and Planetary Interiors，146（1-2）：319-332.

Davaille A.1999.Simultaneous generation of hotspots and superswells by convection in a heterogeneous planetary mantle.Nature，402（6763）：756-760.

Gordon R G，Jurdy D M.1986.Cenozoic global plate motions.Journal of Geophysical Research，91（B12）：12389-12406.

Gurnis M.1988.Large-scale mantle convection and the aggregation and dispersal of supercontinents.Nature，332（6166）：695-699.

He Y M，Wen L X.2009.Structural features and shear-velocity structure of the“Pacific Anomaly”.Journal of Geophysical Research，114（B2）：B02309.

Ishii M，Tromp J.1999.Normal-mode and free-air gravity constraints on lateral variations in velocity and density of Earth′s mantle.Science，285（5431）：1231-1236.

Karato S I.2010.Rheology of the deep upper mantle and its implications for the preservation of the continental roots：a review.Tectonophysics，481（1-4）：82-98.

Kennett B L N，Widiyantoro S，Van Der Hilst R D.1998.Joint seismic tomography for bulk-sound and shear wave speed in the Earth′s mantle.Journal of Geophysical Research，103（B6）：12469-12493.

Khan A，Boschi L，Connolly J A D.2011.Mapping the Earth′s thermochemical and anisotropic structure using global surface wave data.Journal of Geophysical Research，116（B1）：B01301.

Korenaga J.2005.Firm mantle plumes and the nature of the coremantle boundary region.Earth and Planetary Science Letters，232（1-2）：29-37.

Lassak T M，McNamara A K，Garnero E J，etal.2010.Coremantle boundary topography as a possible constraint on lower mantle chemistry and dynamics.Earth and Planetary Science Letters，289（1-2）：232-241.

Lassak T M，McNamara A K，Zhong S J.2007.Influence of thermochemical piles on topography at Earth＇s core-mantle boundary.Earth and Planetary Science Letters，261（3-4）：443-455.

Li Z X，Zhong S J.2009.Supercontinent-superplume coupling，true polar wander and plume mobility：Plate dominance in wholemantle tectonics.Physics of the Earth and Planetary Interiors，176（3-4）：143-156.

McKenzie D，Jackson J，Priestley K.2005.Thermal structure of oceanic and continental lithosphere.Earth and Planetary Science Letters，233（3-4）：337-349.

McNamara A K，Zhong S J.2004.The influence of thermochemical convection on the fixity of mantle plumes.Earth and Planetary Science Letters，222（2）：485-500.

McNamara A K，Zhong S J.2005.Thermochemical structures beneath Africa and the Pacific Ocean.Nature，437（7062）：1136-1139.

Moresi L N，Solomatov V S.1995.Numerical investigation of 2D convection with extremely large viscosity variations.Physics of Fluids，7（9）：2154-2162.

Murakami M，Hirose K，Kawamura K，etal.2004.Post-perovskite phase transition in MgSiO3.Science，304（5672）：855-858.

Nakagawa T，Tackley P J.2005.The interaction between the postperovskite phase change and a thermo-chemical boundary layer near the core-mantle boundary.Earth and Planetary Science Letters，238（1-2）：204-216.

Ni S D，Tan E，Gurnis M，etal.2002.Sharp sides to the African superplume.Science，296（5574）：1850-1852.

Oganov A R，Ono S.2004.Theoretical and experimental evidence for a post-perovskite phase of MgSiO3in Earth′s D″layer.Nature，430（6998）：445-448.

Peslier A H，Woodland A B，Bell D R，etal.2010.Olivine water contents in the continental lithosphere and the longevity of cratons.Nature，467（7311）：78-81.

Phillips B R，Bunge H P.2005.Heterogeneity and time dependence in 3Dspherical mantle convection models with continental drift.Earth and Planetary Science Letters，233（1-2）：121-135.

Piper J D A.2010.Protopangaea：Palaeomagnetic definition of Earth′s oldest（mid-Archaean-Palaeoproterozoic）supercontinent.Journal of Geodynamics，50（3-4）：154-165.

Solomatov V S.1996.Can hotter mantle have a larger viscosity？Geophysical Research Letters，23（9）：937-940.

Tackley P J.2011.Dynamics and evolution of the deep mantle resulting from thermal，chemical，phase and melting effects.Earth-Science Reviews，110（1-4）：1-25.

Tackley P J.King S D.2003.Testing the tracer ratio method for modeling active compositional fields in mantle convection simulations.GeochemistryGeophysicsGeosystems，4（4）：8302.

Tan E，Leng W，Zhong S J，etal.2011.On the location of plumes and lateral movement of thermochemical structures with high bulk modulus in the 3-D compressible mantle.Geochemistry Geophysics Geosystems，12（7）：Q07005.

Turcotte D L，Schubert G.2002.Geodynamics.Cambridge：Cambridge University Press.

Van Keken P E，King S D，Schmeling H，etal.1997.A comparison of methods for the modeling of thermochemical convection.Journal of Geophysical Research，102（B10）：22477-22495.

Yoshida M.2010.Preliminary three-dimensional model of mantle convection with deformable，mobile continental lithosphere.Earth and Planetary Science Letters，295（1-2）：205-218.

Zhang N，Zhong S J，Leng W，etal.2010.A model for the evolution of the Earth′s mantle structure since the Early Paleozoic.Journal of Geophysical Research，115（B6）：B06401.

Zhong S J，Zhang N，Li Z X，etal.2007.Supercontinent cycles，true polar wander，and very long-wavelength mantle convection.Earth and Planetary Science Letters，261（3-4）：551-564.