板块构造与地幔流的动力学：从局部尺度到全球尺度

G.Stadler M.Gurnis C.Burstedde L.C.Wilcox L.Alisic O.Ghattas

板块构造与地幔流的动力学：从局部尺度到全球尺度

G.Stadler M.Gurnis C.Burstedde L.C.Wilcox L.Alisic O.Ghattas

板块构造受控于集中在板块边界的推动力和阻力，然而由观测资料所约束的高精度全球地幔流模型仍面临着计算上的挑战。我们利用最新的自适应网格细化算法，将板块边界的分辨尺度下降到1km，通过分析相接板块的运动，在并行计算机上模拟全球地幔流。在上地幔中，当消减板块向下俯冲时，往往会产生弧后扩张和板片回退现象。下地幔中冷的热异常通过狭窄的高粘性的板块与洋壳相耦合，导致洋壳运动速度的下降。在海沟区域，弯曲岩石层内的粘性耗散占整个岩石层和地幔总耗散的5%至20%。

地幔对流和与其息息相关的板块构造主要影响和控制着地球的热量和地质过程的演变。这些过程恰恰是我们了解构造形变的起源和发展、地球上热量、地幔的成分状态的变化，以及最终的目的——了解地球演变过程的关键。板块的生成和运动很大程度上决定了固体地球的热量散失（Davies，1999），同时板块间的作用力决定着地质时期的能量耗散和热流失（Conrad and Hager，1999）。然而，尽管板块活动具有如此重要的意义，但推动或阻碍板块运动的力仍有很多基本的问题目前还没有确切的答案。

虽然学者们一致认为，俯冲板片在其内部热浮力的作用（Hager and O'Connell，1981），以及上地幔中固－固相变（Billen，2008）、板块从洋脊移动至海沟过程中逐渐冷却的扰动下，以每年1～10cm的速度移动，但对下地幔中板片无震扩张现象（Lithgow－Bertelloni and Richards，1998）在板块运动中的重要性却还没有明确的认识。板块俯冲产生的力可能控制着板块运动的速度。大部分受有效负浮力控制的庞大的板块均位于过渡带和下地幔，当发生俯冲时，如果俯冲板片较为坚固，推力将在海沟处直接与大洋板块的边缘相耦合（Zhong and Gurnis，1995）；但如果在开始俯冲并且海洋岩石层在海沟下方发生弯曲时海洋岩石层较为坚固，那么在狭窄的高粘度板块内会发生能量损失，这种能量的损失影响和限制了板块运动的速度（Conrad and Hager，1999）。为了解板块俯冲的物理过程，研究者利用二维和三维笛卡尔模型（Billen，2008；Zhong and Gurnis，1998； Cizkovaetal，2002；Schellart，2004；Faccennaetal，2007）以及包括观测结构的有限区域模型（Billenet al，2003；Jadamec and Billen，2010），论述了板块边界和板片强度在这个过程中的重要作用，但如何将真实流变特性引入含有狭窄板块和板块边界的模型，这又是全球地球动力学的一个难以解决的问题。俯冲板片到底是坚硬的还是脆弱的这个问题也没有确切的定论（Billen，2008）。

如果引入坚硬的板块，并用现实方法处理板块边界，那么观测约束模型的能力将会大大提高。观测数据约束了板片的形变，同时也证明了全球模型是可以利用的，比如包括通过深源地震获得的应变速率和板片内部的应力状态（Isacks and Molnar，1971），以及消减带发生的板片回退这一动力学现象和现今的弧后扩张现象（Dewey，1980）的模型。地球表面的一大部分 （全球约15%）并不符合刚性板块构造的模型，但却存在比板块边界处更为严重、相当于海沟的形变（Gordon and Stein，1992）。一些海洋板块的内部存在着弥散型形变，特别是印度—澳大利亚板块（Gordonetal，1998）。我们掌握了从地方尺度到区域尺度的丰富的测地、地形、重力和地震数据，通过这些数据可以约束形变。如果我们能够利用与这些数据具有相同精度的模型，并且考虑到真实的粘滞参数，则可以更深入地研究全球动力学。

正如目前所看到的一样，在研究过程中，我们遇到的最大困难并不是获得观测结果，而是数据的处理和计算。如果将真实的粘滞系数和当地的地质情况引入模型中进行处理，从以往的经验来讲是行不通的，因为没有较好的数学方法可用，同时也没有足够的计算机资源。总而言之，目前这一代模型无法充分利用现今地壳形变的观测数据来约束结果。例如，板块运动模型不能使用板块边缘和板块内部的形变观测结果，因为如果将全球地幔对流模型的尺度缩小至与板块边界断层这么大，就会造成时间量程和长度量程过大，从而无法处理。使用自适应网格细化（AMR），并利用并行式计算机处理数据，则可以将真实的粘滞系数和精密标度的观测结果整合进全球地幔的对流模型，再将观测数据带入模型，从而逐步接近板块运动动力来源以及固体地球能量散失等未解问题的谜底。

1 并行自适应网格细化

为了找出板块边界中粘滞性变化巨大的规律，需要用到网格边长约为1km的局部处理方法。如果用全球统一大小的网格将边长设置为1km，则需要个数为1012数量级的网格单元。这不仅超出了现今最强大的超级计算机的处理能力，同时也超出了数学方法处理数据的能力范围。倘若使用自适应网格细化，我们可以在接近板块边缘的区域使用1km精度的处理方法，而在热边界层内则使用相对粗糙的5km精度处理（包括大洋岩石层），在剩余部分的地幔使用15～50km精度处理。自适应网格细化相对于统一网格规格的处理方法，节省了超过103倍的空间。正是由于节省了这些空间，使得这项工作可以在千万亿次超级计算机上实现。

但如果将自适应网格细化方法分配到几千个处理器中进行计算，则是具有挑战性的（Diachinetal，2006）。因为自适应网格细化需要动态地、不规则地更新网格之间的拓扑关系，但在并行计算机上处理时，每一个处理器中只存储一小部分网格信息。所以为了保证每个处理器处理负载保持平均，在每一步细化或粗化过程后都需要重新刷新每个处理器的网格信息。为此我们设计了用来控制网格分配的可扩展算法，以及在我们的自适应网格细化有限单元框架中用于处理地幔流等式的数学方法（Bursteddeetal，2008，2009，2009）。我们的并行自适应网格细化算法基于一个由自适应八叉树组成的森林。这些八叉树使用了曲面多立方体来组成整体的几何构型，球壳用于表示地幔，它又由24级用来表示粘滞特性的不同颜色的立方体组成（图1）。每个立方体可以通过八叉树结构自适应地再次进行细分，这样就可以利用运行速度较快的算法控制网格细化，并构成数字模拟中需要的网格连通信息。这些算法分配在20多万个处理器中运行，适用于多精度、范围广阔并且要求局部区域高精度处理的（或过程动态变化的）科学和工程问题，例如涉及近锋面、不连续面、材质分界面、内隅角（凹角）、边界和内部层等方面的问题。

图1 （a）地球的地幔分解为若干个24个粘性等级的曲面立方体。图例中的颜色表示不同的有效粘度。狭窄的低粘度区域，即板块边界，在图中地球的表面用红色曲线表示。（b）澳大利亚板块转枢区〔即图c中用白色矩形圈出的区域〕的放大图。图中显示的自适应网格具有1km的高分辨度。（c）剖面图，图示了在板块边界处以及对非线性粘滞性动态反应的改良。从图c中还能看出，在西南太平洋板块从新赫布里底到汤加的区域内出现了塑性破裂。其中澳大利亚板块标为AUS，新赫布里底板块标为NH，汤加板块标为TO，太平洋板块标为PAC

2 流变学和同化约束

我们求解瞬时地幔流时，将其等同于求解一个质量和动量相平衡的过程（斯托克斯流），在这个过程中，假定固体地球内部温度（密度）分布已知，则可以确定流速、应变速率和应力状态。岩石层和地幔的扩散蠕变、位错蠕变以及塑性屈服点决定了它们的粘度。实验的约束条件为我们提供了本构关系的形式（资料和方法见www.sciencemag.org/cgi/content/full/329/5995/1033/DC1）和未知参数的取值范围，我们可以在取值范围内进行数值的变动。根据这些关系，可以判断位错蠕变主要出现在上地幔区域。在整个地幔范围内都会出现扩散蠕变，但它仅在下地幔区域引起形变。在分析中同时也结合了塑性屈服点的数值，并根据塑性屈服点给定介质所能承受的最大强度。

我们整合板块边界的方法与前人无异，都是将其看作一个狭窄的区域，并用折合粘度来描述它。在本构关系中，通过一个空间依赖的前因子来表示（Billen，2008）。在处理洋脊和转换边界时，我们将其看作浅层垂直位面，这些洋脊和转换边界的位置参考其他文章（Bird，2003）。对于消减带，我们使用全局优化方法处理逆冲地震的震源机制，将其限制为一个与板块边界走向方向一致的发生连续形变的平面（www.sciencemag.org/cgi/content/full/329/5995/1033/DC1 Materials and Methods）。消减带的深度有时深于地震发生的最大深度，一直延伸到地幔楔，因为从俯冲板片中释放出的水可能会在地幔楔中板片的上部产生倾斜的低粘度区域（Hebertetal，2009）。

全球模型完全由温度变化所控制，而温度变化又基于3个分量。一个分量基于地表热年龄，一个基于确定有地震活动的板片，另一个基于层析地震成像（www.sciencemag.org/cgi/content/full/329/5995/1033/DC1Materials and Methods），这与最近低分辨度模型的输入类似（Ghoshetal，2010；Yoshida and Nakakuki，2009）。我们在建立上地幔的结构时，之所以使用地震活动的资料而非层析成像的数据，是为了最大程度将上地幔描述得更为精细。总体来讲，在目前的层析成像模型中，板片基本无法进行分辨，因为板片和地幔楔常常连成模糊的一片。要正确解决板块运动问题，将板片看做是高粘度的应力导层是非常重要的（Zhongetal，1998）。

3 板块运动和板性

即使使用了由实验确定的流变系数，并通过地震活动确定了地幔结构，但也不能保证刚性板块具有非线性和塑性流变的特性。因此，除了使用往常惯用的单欧拉极来表示单个板块，我们也要评估板块的表面速度与刚性板块的表面速度有多接近，即板块的板性。然而我们发现，为了使模型符合板块运动，需要将屈服应力的实验值大幅度调低（调至大约100MPa），从先前的研究结果也可以看出这种调整并不出人意料（Zhonget al，1998）。在正常情况下，只考虑上地幔板片的浮力和板块冷却产生的板块运动，与在无净旋转参考系中的板块运动大体上相一致（Argus and Gordon，1991）。但是我们发现，无论在小尺度还是大尺度的情况下，表面粘度和应变速率上的变化都非常剧烈，在消减板块边缘，应变速率甚至在局部从10－14/s变化到10－17/s量级（图2a）。消减板块的边界与低应变率的上覆板块相邻，高应变速率的区域围绕着消减板块的边界分布。尽管从整体看来较小（图2a），但是这种变化在板块间断层与在地表显示相比，产生的波长更长。一些模型设定的倾斜区域的规模甚至能和板块的厚度相当。在这种模型中，消减板块在穿过软弱区时渐渐熔化，而实际上这是由于模型的精度过低而产生的人为假象。真实的情况是，岩石层由于断层而滑动，因此板片的粘度和应变速率才会动态地发生变化。所以，板片和板块的自然连接与之前的模型有本质的区别。在应变速率和粘度中出现的大尺度的变化是由板块弯曲引起的，同时发生形变区域的宽度是由消减板片和地幔深层的非线性流变性质和两者之间的耦合程度所决定的（图3d和图4d）。

从图2（b）可以看到，正常情况下，与其他一些线性模型相比，板块的运动速度更符合观测结果（Lithgow－Bertelloni and Richards， 1998； Becker and O'Connell，2001）。主要大洋板块的运动速度比上覆板块快，因此与线性模型的对称汇聚不同，真实的板块汇聚是非对称过程（Lithgow－Bertelloni and Richards，1998；Becker and O'Connell，2001）。由于俯冲板片具有复杂的流变特性，故上地幔中的板块可以将大洋板块渐渐拉向海沟。尽管最后的结果与线性模型中的参数化法相差不多（Conrad and Lithgow－Bertelloni，2002），但是浮力与弯曲板块中非线性阻力相作用不断地达到自洽的平衡状态，从而对大洋板块产生拉力这个过程，使我们更清楚地看到了板块运动中推力和阻力的产生和相互作用。

图2 以180°W为中心3种模型下的应变速率、板块运动速度和板性。（a，b和d）模型1：只包括板块冷却和上地幔板片。（c）模型2：与模型1大致相同，但没有下地幔横向结构。（e，f）模型4：与模型2类似，只是将n设定为3.5。（a）为第二应变速率不变量。（b，c和f）表示在NNR参考系中板块运动的速度。其中绿色箭头表示从参考文献（Argus and Gordon，1991）得到的速度矢量；黑色箭头表示预测速度矢量；实际板块的边界用红色、灰色和蓝色实线表示。（d，e）分别用两种方式表现了太平洋板块的板性，其中图（d）为用右侧色标表示的计算速度和由最佳适配欧拉极 P2（www.sciencemag.org/cgi/content/full/329/5995/1033/DC1Materials and Methods）所得到的速度的矢量差；图（e）是由右侧色标表示的旋转量ω的半径为20°的球冠内的欧拉极独立推测出的板性。Nuvel1－NNR模型的极点位于图中红色三角形处，太平洋板块内计算速度的最佳适配极点在黑色矩形内

当坚硬的板块俯冲至下地幔中时，板块运动会发生本质的变化，通过层析成像得出的过渡带和下地幔温度变化幅度增大的现象可以证明这一点。这种变化会引起整个地幔流运动速度和大洋板块运动速度的下降（图2b 和 2c；见 www.sciencemag.org/cgi/content/full/329/5995/1033/DC1Materials and Methods），然而这个现象是与线性模型截然相反的（Lithgow－Bertelloni and Richards，1998）。另外，随着板块间缝隙变窄、板块边缘介质粘度变高，地表应变速率的范围会有所下降。这可能是板片俯冲至过渡带，以及下地幔中介质的粘度与温度相关（见www.sciencemag.org/cgi/content/full/329/5995/1033/DC1Materials and Methods）等原因综合造成的结果。智利的消减边界的横截面（图3d）显示了在俯冲板片和下地幔之间的粘度状态并不连续，存在着一个间断。当过渡带和下地幔温度发生变化时，高粘度的板片嵌入到高粘度的下地幔中，深层粘度更高的介质便阻碍地表板块发生运动。这时的下地幔结构将会阻碍或终止海沟区域的非对称汇聚过程（图2c）。尽管俯冲板片具有非线性的流变性质，但由于俯冲板片的粘度很高，它和下地幔高粘度介质产生强烈的耦合过程会增加或降低海洋板块的运动速度。这就意味着在深层地幔运动和板块运动之间存在着某种平衡关系。

通过研究本构关系中未知参数的取值范围 （包括上地幔和下地幔的活化能、火成岩粒度、屈服应力和非线性指数幂n），我们可以获得刚性板块的运动过程（图2d）。除了有较大板内应变的新西兰南部区域外，在太平洋板块（PAC）我们可以找到几乎没有板内应变的近乎理想的板块运动。太平洋板块的刚性是对非线性流变性比较敏感的函数。如果流变性的非线性过强，太平洋板块的板性将会降低到不可接受的程度。例如，我们假设n＝3.5，同时增大活化体积和下地幔的浮力，在模拟的结果中太平洋板块出现了大量的应变，导致了汇聚边界的形变延伸超过103km，并且太平洋板块南部和北部的欧拉极之间间隙很大（图2e）。n的取值越接近3.5，越接近上述实验的结果（Hirth and Kohlstedt，2003）。我们发现，当n的取值大于3时，板内应变增大得很快，而且板性与观测数据明显不符。我们还发现，当降低上地幔的粘度，与通用的模型一致时，板性会有所改善。印度－澳大利亚板块被一块广阔的形变区域分成两部分，这与观测数据相一致。但与太平洋板块相比，印度－澳大利亚板块的板性对非线性流变性较不敏感。

4 板片形变和板间耦合

标准模型实际上过高地估计了纳斯卡板块（NAZ）向东运动的速度以及南美板块（SAM）向西运动的速度（图2b）。通过线性模型，可以预测向秘鲁—智利海沟的对称汇聚运动，但这与通过板块动力学无净旋转参考系得出的南美板块几乎静止不动的结论相违背。加入上面所述的下地幔结构，会使两个模型发生退化，产生两个板块的减速现象，使其向海沟处的对称汇聚（图2c）。位于南美板块下方的纳斯卡板块的横截面显示，应力状态也与深源地震的情况不符——在深源地震中，最深处发生的地震往往是由于倾斜挤压造成的断裂引起的，而200～300km深处的地震则是由拉张断裂引起的。在岩石层，尽管海沟附近的消减板块处于强烈的扩张状态中，但南美板块却仅仅受到应力主轴为南北方向的轻微挤压作用。沿智利海沟发生的特大地震表明这个区域是一个耦合作用很强的板块边界。通过增加纳斯卡板块和南美板块之间软弱区的粘度，我们同时改变了大尺度板块运动以及下地幔板片和上覆板块的应力状态。与只有板块运动和板块冷却的正常情况相比，增加板间耦合将降低纳斯卡板块向东移动的速度，同时也会停止南美板块向西运动，并使安第斯山脉内的挤压更加剧烈。除此以外，全球板块运动的过程并没有发生改变。

图3 （a）模型1的速度矢量和表面应变速率。（b）模型6中哥伦比亚、秘鲁和智利海沟处软弱区粘度增加产生的影响。白色箭头是Nuvel1－NNR模型的结果（Argus and Gordon，1991），黑色箭头是预测的结果。（c，d）用红线表示的剖面位置。（c）为根据矩心矩张量解得出的下倾压缩（蓝色）和下倾扩张（红色）分布（www.sciencemag.org/cgi/content/full/329/5995/1033/DC1Materials and Methods）。（d）模型6中预测的压缩轴。其中NAZ指纳斯卡板块，SAM指南美板块

除了与海陆板块边界的形变保持一致，我们发现我们的模型还与形变、在海沟处发生的快速 “回退”和一些洋—洋消减带特有的现象（例如弧后扩张）相吻合。全球或局部模型不能较好地预测这些形变。一般的海沟回退模型往往将上覆板块移除，将地幔延伸到地表，使俯冲板片被很热的地幔包裹起来（Schellart，2004；Faccenna et al，2007）。但是我们所使用的方法与一般模型不同，其中包含上覆板块。我们发现新赫布里底俯冲板片和新赫布里底（NH）板块几乎总是发生回退现象（俯冲板片并没有穿入转换带，而是垂直地插入地幔，产生拉应力作用于相邻的澳大利亚板块），俯冲板片总是处于斜倾的拉张状态中，与观测到的地震性质相一致。

汤加—克马德克地区因发生快速回退现象而闻名。从一些模型里可以看出，在它下方的贝尼奥夫带存在明显的下倾压缩。这个现象对上覆板块的流变特性结构可能有明显影响。早期的模拟分析使用两个数据集作为输入条件——板块边界位置（Bird，2003）和板块年龄（Mulleretal，2008）。我们在分析时发现，一些板块边界的位置与互相不重叠的若干个山脊不符，并且在劳乌海盆和豪夫利海槽中最年轻的岩石层的面积并没有想象中的那么大（www.sciencemag.org/cgi/content/full/329/5995/1033/DC1Materials and Methods）。我们同时发现克马德克（KE）和汤加（TO）微板块与澳大利亚板块（AUS）互相闭锁，不管如何变化n、屈服应力以及上地幔、下地幔的粘度系数，汤加海沟都没有出现回退。但是如果矫正山脊位置和上覆板块的年龄等数据，使其与当地地质资料相符 （www.sciencemag.org/cgi/content/full/329/5995/1033/DC1Materials and Methods），克马德克和汤加与澳大利亚板块之间的耦合将会减弱，再加上克马德克和汤加两个板块快速向东运动，汤加—克马德克海沟就很有可能出现快速回退（图4a）。具体来说，克马德克板块围绕着位于北岛和新西兰附近的欧拉极向东转动，就像实际观测到的那样。

图4 （a）新赫布里底—汤加区域基于无净旋转参考系的表面应变率和板块的运动速度。预测结果用黑色箭头表示，Nuvel1－NNR模型的结果用白色箭头表示（Argus and Gordon，1991），同时还包括Bird（2003）中的微板块运动。其中澳大利亚板块用AUS标示，新赫布里底板块为NH，汤加板块为TO，克马德克板块为KE，太平洋板块为PAC。（b，c和d）图（a）中用红线表示的剖面。（b）预测速度和粘度的分布。（c）由矩心矩张量解得出的下倾压缩（蓝色）和下倾拉张（红色）分布。（d）模型6中预测的压缩轴

另外，海沟回退模型中显示的应力状况与从震源机制中推测出的应力状况相符合。模型中的板片倾斜地挤压穿入过渡带，在过渡带和上地幔之间，板片的顶端附近承受着强大的挤压作用，而在转枢区则承受拉张力。以上这些都与震源机制一致（图4c～d）。在转枢区附近，介质的粘度在100km内降低了约10倍，而沿着海沟的轴线下降得更剧烈，这与克马德克海沟和外槽壁出现的板块局部刚度大幅度下降现象不谋而合（Billen and Gurnis，2005）。由于地幔楔的粘度很低，再加上板片附近的压力梯度很大，因此地幔楔的流速约是板片的5倍，这也与最近的局部地区研究成果相符（Jadamec and Billen，2010）。

5 全球模型的能量耗散

固体地球内部能量耗散的位置目前还不能确定。一些分析将能量耗散的位置定在转枢区 （Conrad and Hager，1999；Buffett and Rowley，2006），而其他一些文章的作者却不这样认为（Wu et al，2008；Davies，2009）。但是先前的模型预先假设了板块运动和转枢区的粘度。通过分析震源机制和海沟回退，应力状态和应变速率可以通过观测来约束，而不依赖于板块运动。

我们发现在俯冲板片内部、地幔楔和下地幔中的板片中，能量的耗散很高（图5）（www.sciencemag.org/cgi/content/full/329/5995/1033/DC1Materials and Methods）。在俯冲板片中，耗散最高的地点是转枢区，但是整个板片的耗散更高，转枢区并不是其耗散的主要来源；另一个耗散很高的地方是过渡带内部，与地震活动发生时能量释放的峰值相同（Abe and Kanamori，1979）。板块间的低粘度区也是能量耗散巨大的区域。出人意料的是，尽管地幔楔的粘度非常低（Billen et al，2003），但由于其应变速率很高（图4b）（Jadamec and Billen，2010），因此能量的耗散也非常高（图5c）。

图5 在小圆的17°S截取剖面的逐点能量耗散。（a）为模型6（无下地幔结构）的情况，（b）为模型8（有下地幔结构）的情况。汤加地区用大写字母T表示，智利北部用Ch表示。（c）模型6汤加地区的放大图。（d）模型8中汤加地区的放大图

我们发现与整个地幔和岩石层的能量耗散相比，转枢区的能量耗散只占一小部分。即使假设下地幔中不存在浮力，转枢区的能量耗散也仅占地幔内部粘滞耗散不到20%的份额（可见www.sciencemag.org/cgi/content/full/329/5995/1033/DC1Materials and Methods）。下地幔是能量耗散的主要源区 （图 5）（www.sciencemag.org/cgi/content/full/329/5995/1033/DC1Materials and Methods），占总能量耗散的70%。下地幔能量耗散的量由下地幔的粘度及其地层走向的横向变化决定，其中横向变化又依赖于随温度而变化的粘度。如果是这样的话，转枢区的能量耗散则低于5%。转枢区只是地幔和地壳中由若干成分紧密连接的系统的一个部分，如果考虑整个地幔和地壳，这些统计数值都要低于早期文章的估计结果（Conrad and Hager，1999；Buffett and Rowley，2006），但是这与最新的研究论证却是一致的（Leng and Zhong，2010）。

译自：Science.2010.329：1033～1038

原 题：The dynamics of plate tectonics and mantle flow：from local to global scales

（中国地震局兰州地震研究所 赵雯佳译；云南省地震局 余丰晏校；吕春来复校）