俯冲板块的动力学性质对平板俯冲的影响：数值模拟

朱致远，吴本君

南京大学 地球科学与工程学院，南京 210023

根据弧后应力状态，科学家们将俯冲带分为马里亚纳型俯冲带和秘鲁—智利型俯冲带两类（Uyeda and Kanamori, 1979）。秘鲁、智利附近的俯冲带是秘鲁、智利型俯冲带的典型代表，纳兹卡板块向南美洲板块下方快速俯冲。在墨西哥中部、秘鲁和智利中部，俯冲板块在南美大陆岩石圈下形成水平的俯冲区域（Manea et al., 2017）。把俯冲板块在大陆岩石圈下出现水平片段的俯冲叫做平板俯冲（flat-slab subduction）。

平板俯冲与一系列特殊地质现象有着紧密联系，往往伴随着强烈的上覆大陆岩石圈挤压和地壳增厚，有规律的岩浆活动（Gutscher et al., 2000a;Kay and Mpodozis, 2001; Sigloch et al., 2008），以及俯冲带上覆岩石圈的热异常结构 （Dávila and Lithgow-Bertelloni, 2015）。例如，Gutscher等(2000a)认为南美洲西部出现的与年轻大洋俯冲相关的第四纪埃达克火山岩可以通过平板俯冲模型进行解释。Li等（2007）使用平板俯冲模型解释华南中生代的造山运动中浅海盆地的形成和相关的岩浆活动。Sigloch 等（2008）则通过长距离平板俯冲模型来解释北美洲西部75 Ma左右大范围岩浆活动和从海岸到大陆内部的长距离造山运动。此外，相较于高角度的俯冲，平板俯冲还会阻碍超高压岩石折返（Li et al., 2011）。

近年来，对平板俯冲形成机制与影响因素的研究认为浮力是影响平板俯冲形成的第一要素（Pilger,1981; Skinner and Clayton, 2011），年轻的板块与增厚的洋壳会使得俯冲板块的平均密度较小，使其难以下沉从而使得平板俯冲更易于发生。此外，特殊的上覆岩石圈热力学结构 （Rodríguez-González et al., 2012）、上覆岩石圈向海沟的运动（Lallemand et al., 2005）、俯冲板块的流变性质（Van Hunen et al.,2002; Billen and Hirth, 2007; Manea and Gurnis, 2007）、板块与地幔的耦合作用（Boutelier and Cruden, 2008）等都影响平板俯冲形成。

俯冲板块的厚度、粘度、密度等动力学性质是决定其应力应变状态的重要因素，很大程度上约束平板俯冲的形成（Huangfu et al., 2016）。然而，俯冲板块的动力学性质如何影响平板俯冲形成以及形态的变化仍需要进一步的深入研究。下沉的俯冲板块在上覆的大陆岩石圈下形成平板区域的长度和角度变化与俯冲板块的动力学性质的关系仍不明确。本文主要研究俯冲板块的动力学性质对于平板俯冲形态的制约。通过二维地球动力学数值模拟，着重探讨俯冲板块的厚度、粘度、密度等因素对俯冲过程中形成的上覆岩石圈下方的平板的长度和角度变化的影响，并结合智利中部的平板俯冲实例进行了讨论，从而加深了人们对于平板俯冲的动力学过程的认识。

1 数值方法和模型设置

1.1 数值模拟方法

本文建立的俯冲动力学的二维数值模型是基于不可压缩的斯托克斯流体，其控制方程为无量纲动量守恒方程和质量守恒方程：

其中，p是动压力，∇ρ是俯冲板块密度与地幔密度的差值，g是重力加速度，u是速度。应力张量τij与粘度η和应变率张量的关系如下：

其中，

简单起见，由于模拟的俯冲过程的时间比较短（小于35 Ma），笔者的模型实验设计中没有考虑温度和相变。使用联合拉格朗日粒子法的有限元模拟程序Underworld2（Moresi et al., 2019;Stegman et al., 2006）来求解上述动力学方程。

俯冲板块在俯冲过程中，受到的应力超过屈服应力时，会发生塑性形变。笔者的模型采用冯·米赛斯屈服准则，将粘性流变和塑性流变结合，形成实际的有效粘—塑性流变关系：

其中，地幔的参考粘度η0= 1×1020Pa·s，屈服应力τyield=48 Mpa，是应变率张量的第二不变量。

1.2 模型设置

图1 平板俯冲数值模型设置的示意图（虚线方框放大显示了初始俯冲带的特征）Fig. 1 Setup of the numerical model of flat-slab subduction (Dash-line box shows a three-layer subducting slab and a viscous overriding plate)

如图1所示，此次研究的模型长度为9000 km、深度为1000 km，被分为720×80个网格。每个网格内随机分布20个粒子，保存了物质的密度、粘度等性质，并可用于追踪模型中物质的运动历史。上地幔的粘度为ηum=1×1020 Pa·s，与参考粘度η0相同。模型中，下地幔从660 km深度延伸至1000 km，其粘度ηlm为上地幔粘度的100倍。模型中的左右两侧边界是周期性边界，下边界是无滑移边界。上边界的俯冲板块和上覆板块上施加了相对运动的速度。俯冲板块以3 cm/a速率，上覆板块以4.5 cm/a速率同时向海沟方向运动。主要的模型参数见表1。

参考模型中俯冲板块的长度为4500 km，厚度为70 km，对应大洋岩石圈年龄约为35 Ma。俯冲板块密度比地幔密度多60 kg/cm3，大洋岩石圈在负浮力的作用下，向大陆岩石圈下方俯冲。开始时俯冲板块和上覆板块之间放置了50 km粘度与上地幔一致的薄弱区域，以便俯冲启动。上覆板块的长度为3500 km，厚度为120 km。假设上覆板块的密度与地幔密度相同，粘度为上地幔粘度的4000倍。由于俯冲板块在俯冲时受到强烈的拉张和挤压作用，将俯冲板块设置为三层：上下两层为粘塑性层，中间为较硬的粘性层（图1b），从而更好地模拟了板块的应力集中区域发生的塑性形变（Lyu et al., 2019; OzBench et al., 2008）。

表1 模型参数Table 1 Model parameters

2 结果

通过数值模拟实验，此次研究了俯冲板块的动力学性质对平板俯冲形态的影响。首先，分析了参考模型中平板俯冲形态随时间的变化。在参考模型的基础上，考察了俯冲板块厚度、粘度、密度等因素对于海沟运动速度，俯冲板片的平板长度和平板角度等参数的影响。数值模型的主要参数设置如表2所示。

2.1 参考模型

在参考模型REF_MODEL中，俯冲板块的厚度为70 km，粘度为上地幔粘度2000倍。俯冲板块密度比地幔密度大，Δρsp=60 kg/m3。为了更好的描述在上覆板块下方形成的平板的形态，笔者定义了一些变量（表1）：vtr、tfs、Lfs、α、β和θ。vtr是海沟运动速度，向左运动为正，表示海沟后撤。tfs为平板俯冲发生时间，即在俯冲板片上首次出现较为明显的平板区域的时间。平板区域拥有比它两侧更加平缓的倾角以及一定长度。Lfs是平板区域的水平投影距离，简称平板长度。α是平板与水平面的角度，顺时针旋转为正，简称平板角度。β和θ分别为俯冲角度和板片前段弯折角度，顺时针旋转为正。图2是定义Lfs、α、β和θ的示意图。

表2 数值模型Table 2 Numerical models

图3显示了参考模型REF_MODEL中平板俯冲的形成和演化的过程。在推动力与负浮力共同作用下，俯冲板块首先向着上覆板块下方俯冲（图3b）。在约22 Ma时，如图3c所示，俯冲板片前段在上地幔底部滑动，俯冲板块中部出现一段与水平线的角度较小的区域，认为平板俯冲在这时发生。从图4a中可以看到，此时海沟后撤的速率vtr较大，约为4.0 cm/a。平板俯冲刚发生时，平板长度Lfs约为230 km（图4b）。图4还显示了这时平板角度α约为20°，板块俯冲角度β为35°，θ为44°。平板区域的角度与俯冲板片的前后段的弯折角度相差约15°～20°。

随着时间的变化，平板俯冲进一步发展。如图4中实线所示，平板角度α变得更加平缓，平板区域的长度Lfs缓慢增加。海沟运动的速率vtr逐渐减小。板片俯冲角度β和板片前段弯折角度θ变化不大，约为5°～6°。图3d显示了32 Ma时平板俯冲的形态。这时，海沟运动速率减小至3.45 cm/a。平板区域长度Lfs约为245 km，变化不大。平板角度接近0°，与板片前后段有较明显的角度差异（相差约为40°）。

图2 模型参数示意图Fig. 2 Definition sketch of the flat-slab

图3 参考模型REF_MODEL的平板俯冲演化过程Fig. 3 The evolution of the reference model REF_MODEL

2.2 板块厚度的影响

俯冲板块的厚度直接影响了由于负浮力导致的板片拉力。笔者测试了不同的板片厚度对平板俯冲形态的影响（图4，5）。改变厚度时，俯冲板块内部分层厚度保持比例。模型H_50中，俯冲板块的厚度较薄，Hsp=50 km，由于板片负浮力产生的拉力较小，板片俯冲的速度较慢（图5ac）。而且在上覆板块的挤压过程中，薄的俯冲板块发生了较大变形，不利于平板俯冲的发育。板片在约20 Ma时开始发生平板俯冲，平板区域长度Lfs只有约180 km，且随时间变化有较大起伏（图4b）。由于薄板片易变形，平板的倾斜角度α在32 Ma时接近-10°（图4c，负号表示此段板片向上弯折）。图4e显示，板片前段的弯折角度θ也较大，在55°～ 64°之间变化。

图4 俯冲板块厚度对板片平俯冲的影响Fig. 4 Influence of subducting plate thickness to flat slab subduction

图5 模型H_50（a-c）和模型H_90（d-f）随时间的演化过程Fig. 5 Time evolution for the numerical models H_50 (a-c) and H_90 (d-f)

当俯冲板块厚度为90 km时（模型H_90），较厚的板片产生了较大的负浮力，导致板块的俯冲和后撤的速度也更快。因此，俯冲过程中没有形成平板区域（图5d-f）。同时，海沟运动的速率变化也不大，约为4 cm/a（图4a）。如图4c所示，17～ 33 Ma的时间里，板片倾斜角度α维持在30°左右。板片俯冲角度β随时间基本没有变化，与参考模型REF_MODEL相比，差别也不大。随着板片前段在上地幔底部滑移，板片弯折角度θ从41°变化到30°。

通过以上分析，笔者发现太薄和太厚的板片都不利于平板区域的产生。参考模型REF_MODEL中的70 km厚度的俯冲板块较有助于平板俯冲的发生。

2.3 板块粘度的影响

俯冲板块的粘度直接决定了板片的强度，控制了板片变形的难易程度。笔者考察了俯冲板块粘度分别为600ηum，1000ηum和4000ηum时，平板俯冲形态的变化。虽然随着俯冲过程的演化，所有模型中海沟运动速率都在减小，如图6a所示，但是俯冲板块的粘度直接影响了海沟后撤的速率。板片粘度越小，海沟后撤速率越大。当ηsp/ηum=600时（模型V_600），在32 Ma时海沟运动速率为3.7 cm/a，比模型V_4000（ηsp/ηum= 4000）中的速率大0.4 cm/a。

此外，从20 Ma到30 Ma的平板俯冲演化过程中，四个模型的平板长度都随着时间而增加（图6b）。模型V_4000中，平板长度Lfs从235 km增加到了275 km。而模型V_600中的平板长度随时间仅有少量变化，从190 km 增加到了200 km。不同粘度的模型中，平板角度α都随平板俯冲的演化而趋于水平；板片俯冲角度β和前段弯折角度θ随时间变化不大，约在5°以内。

2.4 板块密度差的影响

根据半空间冷却模型，板块厚度和密度都会随着年龄增加而增大。考虑到板片的密度同时受到应力、组分等多方面因素的影响，在未考虑温度场和相变的情况下，将密度和厚度分别进行讨论。前人的研究表明，俯冲板块的密度比地幔密度重，差值约为50～100 kg/m3（Skinner and Clayton,2011; OzBench et al., 2008）。 图7显示了俯冲板块与地幔的密度差对平板俯冲演化过程的影响。从图7a可以看出，如果发生了平板俯冲，如模型D_50和D_70，俯冲板片的密度差的变化对海沟后撤的速率影响不大。板块密度差越小，越容易发生平板俯冲。如图8a-c所示，当Δρsp=50 kg/m3（模型D_50）时，平板俯冲在约15 Ma就发生了，随着时间的演化，在32 Ma时，平板长度Lfs达到313 km。但是，当板片密度差很大时，以Δρsp=80 kg/m3（模型D_80）为例，由于负浮力增大，板片难以被地幔流托起。如图8d-f所示，俯冲过程中没有形成平板区域，板片中段的倾斜角度α在25°～32°之间变化。与模型H_90类似，海沟运动速率的变化也不大，在3.9 cm/a上下浮动。图7d和7e显示不同的俯冲板片密度差对板片的俯冲角度β和前段弯折角度θ基本没有影响。

图6 俯冲板块粘度的影响Fig. 6 Influence of subducting plate viscosity

图7 俯冲板块密度的影响Fig. 7 Influence of subducting plate density

图8 模型D_50（a-c）和模型D_80（d-f）随时间的演化过程Fig. 8 Time evolution for the numerical models D_50 (a-c) and D_80 (d-f)

这些结果与前人的研究一致（Espurt et al.,2008; Huangfu et al., 2016）。

3 讨论

3.1 智利中部的平板俯冲

南美俯冲带是纳兹卡板块向南美洲板块西缘俯冲而形成的。俯冲过程的一个构造结果是板间和板内地震的发生。这些地震形成了贝尼奥夫带，描述了俯冲板片的形状。通过对南美洲板块下方的贝尼奥夫带的研究，科学家们发现俯冲板片在秘鲁和智利下方形成了数百米的近水平区域（Gutscher et al., 2000b）。Manea等（2017）总结和比较了秘鲁和智利平板俯冲带，认为海沟后撤、俯冲板块的构造演化和上覆板块结构的不连续性明显有助于平板俯冲的形成。Chen等（2019）的纳兹卡板块重建模型发现，纳兹卡俯冲在80 Ma前先从南美洲北部（5°S）开始。然后向南传播，大约55 Ma前，纳兹卡俯冲到达南美洲南部（40°S）。该研究还发现，前渊沉积和安第斯山脉的形成与纳兹卡板片和下地幔之间的相互作用有关。

在秘鲁附近，由于俯冲开始时间较早，俯冲板片已经进入下地幔并堆叠（Chen et al., 2019）。在智利中部，地震层析成像显示纳兹卡板片还停留在下地幔顶部（图9），与笔者的参考模型REF_MODEL得到的俯冲板片的形态相似。参考模型REF_MODEL预测的平板片长度约为240 km，与前人研究结果相近（Manea et al., 2017）。智利中部的纳兹卡板块俯冲时的年龄约为38 Ma （Manea et al., 2017），对应的板块厚度约为75 km。由前面的分析可知，这个板片厚度有利于形成平板俯冲。纳兹卡板片前段的弯折角度约为48°（Li et al.,2019），参考模型REF_MODEL得到的弯折角度θ与之相一致。本文的数值模拟实验可以为研究南美洲智利地区的平板俯冲过程和机制提供参考。

3.2 板块运动速度和海沟后撤

皇甫鹏鹏等（2016）的研究发现俯冲大洋板块和上覆大陆板块相对运动的速率对平板俯冲的形成有着决定性的影响。上覆板块朝向海沟的运动速率越大越有利于形成平板俯冲。Coltice等（2017）根据HS3-NUVEL1A模型（Gripp and Gordon, 2002）的结果，统计了全球主要俯冲带两侧板块的运动速度和海沟的运动速度。其中，在智利附近的俯冲带，俯冲的纳兹卡板块的运动速率约为3.3 cm/a，而上覆板块的运动速率约为4.6 cm/a。本文数值模型所用的板块运动速度边界条件也是根据以上结论设置的。

根据Li等（2019）的总结，海沟后撤的模型中俯冲板块与地幔转换带呈较小角度，最终导致板片滞留在地幔转换带处。本文的模型设置中，上覆板块的向海沟逆冲速率大于俯冲板块的运动速率，海沟后撤的速率约为3～4 cm/a，与前人的研究结果一致（Coltice et al., 2017）。值得注意的是，模型H_90和模型D_80没有出现水平板片，它们的海沟运动速率在3.9～4.0 cm/a 之间浮动。其他发生平板俯冲的模型中，海沟运动的速率都随着平板片的形成而减小，有的模型（V_4000）的海沟后撤速度减小到3.2 cm/a。

本文主要关注了俯冲板块的厚度、粘度和密度差等性质分别对平板片形态的影响。简单起见，温度相变等因素并没有被考虑。上覆板块的厚度粘度等性质也有可能对平板俯冲的动力学过程有显著影响（Manea et al., 2012; Huangfu et al., 2016），笔者将在后续的工作中进行更加深入的研究。

4 结论

本文研究了平板俯冲的动力学特征，以及俯冲板块的动力学性质对俯冲平板片形成的影响。通过数值模拟实验，得出了以下初步结论：

（1） 俯冲板块的厚度和密度差对平板片的形成有着决定性的影响。厚度为70 km左右的板片最容易发生平板俯冲。太厚的俯冲板块（Hsp=90 km）比较难发生变形，阻碍平板俯冲的发生。俯冲板块与上地幔的密度差越小，越容易发生平板俯冲。但是，当板片密度差很大时（Δρsp=80 kg/m3），俯冲过程中不能形成平板区域。

图9 智利中部的板片平俯冲Fig. 9 Central Chile flat slab

（2） 俯冲板块的粘度也是平板片形态的主要影响因素。较大的板片强度有助于形成平板俯冲。板片粘度越大，形成的平板片长度越长。

（3）平板俯冲的过程伴随着海沟后撤速率的减小。