大洋俯冲和大陆碰撞沿走向的转换动力学及流体-熔体活动的作用*

李忠海 许志琴

LI ZhongHai and XU ZhiQin

大陆构造与动力学国家重点实验室，中国地质科学院地质研究所，北京 100037

State Key Laboratory of Continental Tectonics and Dynamics，Institute of Geology，Chinese Academy of Geological Sciences，Beijing 100037，China

2015-06-01 收稿，2015-09-05 改回.

1 引言

自然界中，所有的会聚大陆板块都有侧向边界，并与周围的大洋板块或另一个大陆板块相邻。例如，特提斯构造带作为地球上最大的板块俯冲-碰撞造山带之一，它沿走向呈现一种陆-陆碰撞和洋-陆俯冲的相间结构，西起欧洲的阿尔卑斯陆-陆碰撞造山带，经过地中海复杂的板块俯冲带，向东延伸至西亚的扎格罗斯(Zagros)陆-陆碰撞造山带及其东部的马库元(Makrun)洋-陆俯冲带，再继续向东至喜马拉雅陆-陆碰撞造山带，而后与东印度洋板块向东南亚之下的洋-陆俯冲带相连。这种大洋俯冲和大陆碰撞沿走向的空间转换控制着特提斯构造带深部结构及动力学特征沿走向强烈的变化，也控制着其中每个造山带侧向边缘和造山带腹地不同的浅层地质特征和深部结构(Hatzfeld and Molnar，2010;Yin，2010;Xu et al.，2012;Li et al.，2013)。

在洋-陆俯冲和陆-陆碰撞的过程中，流体-熔体活动对其动力学过程具有非常重要的作用(Wang，2010;Faccenda，2014)。流体和熔体可以引起俯冲带物质的弱化及相变，从而很大程度上促进俯冲带中壳幔物质混合及相互作用，进而影响会聚板块之间的耦合、上覆板块岩浆作用以及隧道内高压-超高压变质岩石的折返(Zheng，2012)。我们最近通过二维数值模拟的方法对含水俯冲碰撞模型与无水模型进行了对比研究，发现无水模型中由于缺少流体-熔体活动而产生的弱化作用，上覆岩石圈的变形很小，在碰撞造山过程中的参与程度很低。而流体-熔体的活动可以加强上覆岩石圈的地形隆升和造山作用(Li et al.，2015)。

本文通过建立一系列三维大尺度、高分辨率的洋-陆俯冲和陆-陆碰撞沿走向转换的数值模型，对比研究俯冲碰撞带中流体-熔体活动的作用和影响，并与中东特提斯构造带的深部结构观测进行对比，进而探讨其动力学机制。

2 数值模拟方法

对于板块俯冲碰撞相关的动力学数值模型，一般对三组控制方程进行求解，包括斯托克斯流体动力学方程、物质守恒方程以及热量守恒方程。本文中，对这些方程在不规则的欧拉网格节点上采用有限差分算法进行离散化和近似求解。数值模拟的基本方法可见Gerya(2010)，具体的控制方程、流变模型及参数等可见Li et al.(2013)，本文中不作详细介绍。由于本文侧重点在于对比研究俯冲-碰撞过程中流体-熔体的活动和影响，下面详细介绍数值模型中采用的流体活动模型和部分熔融模型。

2.1 脱水-水化作用和水的运移

俯冲带中的水主要包括两部分。一种是自由水，一种是矿物水。对于自由水，基于前人的系统研究(Gerya and Meilick，2011;Li et al.，2013)，假设沉积物和大洋上地壳玄武岩在地表含2%(质量百分含量)的水，并且随着深度的增加而减少，直至75km 的自由水含量降为零。对此，假设线性的定量化关系式:

其中，XH2O(P0)代表地表的自由水含量2%，Δz 是深度(0～75km)，对于其他的岩石属性(大陆地壳，大洋下地壳，地幔)，假设其自由水含量为零。

对于矿物水，采用了PerpleX(Connolly，2005)所计算的四种典型岩性的含水量随温度和压力的变化，这些数据是依据吉布斯自由能最小化原理，对特定岩石成分和热力学数据计算而得(Gerya and Meilick，2011)。随着俯冲的进行，矿物水析出后保存在新生成的示踪点中，并独立的运移，当其遇到由于水化或部分熔融作用而产生的水不饱和岩性的时候，该活动的水就会被吸收。水的运动速度由该位置处的压强梯度进行计算，具体算法如下:

其中，vx(water)、vy(water)、vz(water)分别是水在x、y 和z 方向的速度。vx、vy、vz是该位置处的岩石速度。vy(percolation)代表流体-熔体向上的渗透速度，根据前人的研究，参考模型中采用vy(percolation)=10cm/a(Gorczyk et al.，2007;Peacock，1990)。在这里，模型假设水生成之后的水平速度与岩石速度相当，由于水和岩石的密度差，导致其速度差异主要体现在垂向上。同时，我们还设计了不同的模型对流体-熔体渗透速率的变化进行了探讨。

2.2 岩石部分熔融作用

数值模型中包含了多种岩石类型的部分熔融计算(Li et al.，2013)。对于水化蚀变的地幔橄榄岩，采用了Katz et al.(2003)的算法和参数。对于其他岩石类型，基于实验岩石学的约束条件，作为一种近似算法，假设部分熔融的体积比例与温度之间存在一种线性关系:

其中，Tsolidus和Tliquidus分别代表特定岩性的固相线温度和液相线温度。

3 初始模型和边界条件

基于要研究的洋-陆俯冲和陆-陆碰撞沿走向转换的动力学背景，我们设计了一组三维大尺度、高分辨率动力学数值模型(图1)。初始模型中，上覆板块是一个均一的大陆板块，而俯冲板块沿走向平分为两个部分，一侧为均一的大洋板块，而另一侧长度为450km 的大洋板块之后与大陆板块相连。大陆板块的上地壳厚度为20km，下地壳厚度为15km;而大洋板块的上地壳玄武岩厚度为3km，下地壳辉长岩厚度为5km。值得注意的是，在地壳表面之上，与自由滑动的模型顶界面之间，设计有一层相对高黏滞度的伪空气层，其与上地壳的接触面用以模拟地形起伏面，该地形起伏面包含了近似的剥蚀和沉积作用(如Li et al.，2013)。对于大陆岩石圈的温度结构，地表是0℃，岩石圈底部为约1325℃，其间采用线性插值的方法。大洋岩石圈的年龄约为30Ma。岩石圈之下的软流圈地幔的温度梯度为0.5℃/km。

图1 初始模型设计(a)三维初始模型，空间尺度为1000 ×680 ×656km，分辨率为2×2 ×2km. 为了更加清楚的显示内部结构，最顶层的20km 被剪切掉. 颜色代表不同的岩石类型，色标如(b):0-伪空气层;1-水;2、3-沉积物;4-部分熔融沉积物;5、6-大陆上、下地壳;7、8-水化的大陆上、下地壳;9、10-部分熔融的大陆上、下地壳;11、12-大洋上、下地壳;13-部分熔融的大洋地壳;14-岩石圈地幔;15-软流圈地幔;16、17-水化、蛇纹石化的地幔;18-部分熔融的地幔. 水化的和部分熔融的岩石类型在图1 的初始模型中不存在，将随着模型的演化而产生Fig.1 Initial model configuration(a)the spatial scale of 3-D initial model is 1000 ×680 ×656km with the resolution of 2 × 2 × 2km. Colors indicate different rock types as in (b). The top layer (y ＞-20km)is cut off for clarity;(b)the colorgrid for different rock types，with:0-air;1-water;2，3-sediment;4-partial molten sediment;5/6-upper/lower continental crust;7/8-hydrated upper/lower continental crust;9/10-partial molten upper/lower continental crust;11/12-upper/lower oceanic crust;13-partial molten oceanic crust;14-lithospheric mantle;15-asthenospheric mantle;16/17-hydrated/serpentinized mantle;18-partially molten mantle. The hydrated and partially molten rocks are not shown in Fig.1，but will appear during the evolution of the model

对于数值模型的速度边界条件，顶部(y=0)和前后两个侧面(z=0 和z=656km)是自由滑动边界。底部边界是渗透性边界，采用近无限深度的外部自由滑动边界条件(Li et al.，2013)，在模型底边界(680km)之下的200km 处，满足自由滑动条件(Gerya，2010;Li et al.，2010)。与普通的自由滑动边界条件相同，外部自由滑动条件也将满足计算区域内的物质守恒。模型的左右两个边界是会聚边界，自初始条件开始，两个边界之间的会聚速率为5cm/a;在10Ma 到15Ma之间，伴随着大陆俯冲碰撞的发生，会聚速率逐渐减低至2.5cm/a;而后保持2.5cm/a 的会聚速率不变。

对于数值模型的温度边界条件，模型顶部为固定温度(0℃);四周的侧面边界条件为水平方向温度梯度为零(即零热流)。底部边界采用的是外部边界固定温度条件，在模型底边界(680km)之下的200km 处，假设一个固定的地幔温度，这样就可以使得在底部渗透边界处，温度和热流都可以据模型演化而动态调整(Gerya，2010;Li et al.，2010)。

4 模拟结果

4.1 参考模型

基于前人研究，流体-熔体向上渗透速率的估计值约为vy(percolation)=10cm/a(Gorczyk et al.，2007;Peacock，1990)，参考模型中我们采用了这个渗透速率，具体算法如前文所述。

数值模拟结果显示在洋-陆俯冲带和陆-陆碰撞带产生了差异巨大的演化特征。在大陆碰撞带，先期俯冲的大洋岩石圈板块在碰撞后发生断离，断离深度约为150km，断离时间为18Ma 左右(图2a)。本文中，板块断离的判定准则为，俯冲板块岩石圈的物质场错断一定距离(＞10km)，其间被软流圈物质充填。板块断离后，残留的俯冲板块向上弯曲，造成大陆板块俯冲角度变小，并推动了俯冲大陆地壳的折返以及部分熔融物质在上覆岩石圈底部的积聚(图2b)。在大洋俯冲一侧，海沟持续的后撤，伴随着俯冲大洋板块的脱水，以及地幔楔中的水化和部分熔融作用。同时，大洋海沟的后撤产生了弧后的伸展以及边缘海盆地的拉张(图2c)。

三维数值模型显示了沿走向的空间差异性俯冲碰撞模式，该模式造成了上覆地壳物质从大陆碰撞一侧向大洋俯冲一侧逃逸(图2b，c)。同时，造成了深部结构的巨大差异性，导致周围的地幔物质穿过板块的断离窗口，从板下地幔进入地幔楔中，而后从大陆俯冲带向大洋俯冲带侧向流动(图2d)。并且，深部地幔流动的方向与地表物质的挤出方向一致。

4.2 不同的流体-熔体渗透速率

为了探讨流体-熔体的渗透速率对于俯冲-碰撞过程的影响，我们建立了两个新的数值模型，其中一个采用渗透速率vy(percolation)=5cm/a，另一个采用渗透速率vy(percolation)=0。

把流体-熔体的渗透速率从参考模型中的vy(percolation)=10cm/a 降低至5cm/a，三维模型的大尺度演化模式变化不大，大洋板块一侧持续俯冲进入地幔之中，而大陆俯冲碰撞一侧发生俯冲板块的断离，伴随着俯冲引发的部分熔融作用以及上覆地壳物质的侧向挤出(图3)。但是，相对于参考模型中18Ma 时、150km 深度处的板块断离(图2a)，降低渗透速率，导致板块断离的时间较晚(～22Ma)，并且深度较大(～200km)(图3a)。

进而，为了探讨一个流体-熔体活动的端元模型，我们将流体-熔体的渗透速率降低为零，即模型中包含了脱水作用，但是脱出的水只能随固体一起运移，而不发生向上的渗透。数值模拟结果显示，大陆碰撞一侧的板块断离的时间和深度进一步增大，达到26Ma 和300km 左右(图4a)。同时，由于俯冲和上覆大陆板块之间较强的耦合作用，上覆岩石圈地幔也有向下俯冲的趋势，从而倾向于形成一种双向对冲碰撞模式(图4a)。

图2 参考数值模型结果该模型中流体-熔体的向上渗透速率vy(percolation)=10cm/a. (a、b)大陆碰撞一侧的模型演化;(c)大洋俯冲一侧的模型演化;(d)等黏滞系数面和速度场分别展示三维模型的内部结构和地幔流动特征Fig.2 Results of the reference numerical modelIn this model，the upward percolation velocity of fluid and melt is prescribed as vy(percolation)=10cm/a. (a，b)evolution of the continental collision side;(c)evolution of the oceanic subduction side;(d)iso-viscosity surface (1022Pas)and velocity field demonstrate the inner structure and mantle flow characteristics of the 3-D model

4.3 无流体-熔体活动

为了进一步对比研究流体-熔体活动的作用，我们又设计了一个无流体-熔体活动的模型，即在参考模型中(图2)取消所有关于水的活动以及部分熔融作用的计算，从而成为一个单纯的热动力学模型。

数值模拟的结果显示，该模型中大陆碰撞一侧虽然也能最终发生板块的断离(图5a)，但时间更晚(35Ma)并且深度更大(400 ～500km)。由于缺少流体和熔体的润滑作用，会聚板块之间具有强烈的耦合作用，从而无论洋-陆俯冲带还是陆-陆碰撞带，上覆岩石圈地幔都会被拖曳而向下俯冲，从而形成一种双向对冲碰撞模式(图5)。

5 讨论

扎格罗斯-喜马拉雅大陆碰撞造山带及其邻区的大洋俯冲带显示了陆-陆碰撞和洋-陆俯冲的相间结构(图6a)，这种大地构造特征成为控制该区域板块会聚动力学的关键。针对该构造背景，我们建立了一系列三维大尺度、高分辨率的数值模型，探讨洋-陆俯冲和陆-陆碰撞沿走向的转换动力学。模型结果显示了强烈的沿走向的差异性，大洋俯冲带中，岩石圈板块持续的俯冲进入地幔之中;而在大陆碰撞一侧，俯冲板块在一定深度处发生断离。在板块会聚速率相同的情况下，板块断离的深度与俯冲带流体-熔体的活动具有负相关关系。在含流体-熔体模型中，断离深度随着流体-熔体渗透速度vy(percolation)的减小，而逐渐增大:vy(percolation)= 10、5、0cm/a 时，对应的断离深度依次是大约150km、200km 和300km。在无流体-活动模型中，板块断离深度最大，达到400～500km。其动力学机制在于，流体-熔体活动可以很大程度的降低周围岩石的流变强度，造成俯冲板块上部及地幔楔黏滞系数的减小，同时能够降低俯冲板块和上覆板块之间的耦合作用。在俯冲大洋板块拉力一定的情况下，流变强度的降低将对板块断离产生促进作用，从而使得俯冲板块在较早的时间以及较浅的深度处发生断离。

图4 流体-熔体渗透速率为零的数值模拟结果该模型中流体-熔体的向上渗透速率vy(percolation)=0. (a)大陆碰撞一侧的模型演化结果;(b)大洋俯冲一侧的模型演化结果Fig.4 Numerical model result with zero upward percolation velocityIn this model，the upward percolation velocity of fluid and melt is zero，i. e. vy(percolation)=0cm/a. (a)evolution result of the continental collision side;(b)evolution result of the oceanic subduction side

流体-熔体活动对俯冲-碰撞带的表层和深部构造样式产生了很大的影响，但沿走向的巨大差异性是所有模型中共存的。这种走向差异性促使周围地幔物质发生复杂的流动，从板下地幔穿过板块断离窗口进入地幔楔，而后从大陆碰撞带向大洋俯冲带侧向运移。这种绕过大洋俯冲板块侧向边缘的地幔流动，可以促使地幔变形以及地幔矿物发生晶格优选定向，进而产生地震波各向异性(Li et al.，2014)，并能够通过地震台站在地表进行观测(例如，Wang et al.，2008)。同时，深部的地幔流动方向与浅部物质从造山带向俯冲带的侧向逃逸方向是大体一致的，这也显示了二者之间可能存在一种相互促进的关系。

图5 无流体-熔体活动的数值模拟结果该模型中取消了所有流体和熔体活动的计算. (a)大陆碰撞一侧的模型演化结果;(b)大洋俯冲一侧的模型演化结果Fig.5 Numerical model result with no fluid and melt activityIn this model，the whole calculations of fluid and melt activity are canceled. (a)evolution result of the continental collision side;(b)evolution result of the oceanic subduction side

图6 中-东特提斯构造带的大地构造特征及深部结构(a)中-东特提斯构造带及邻区地形，其中红色涂抹线表示大陆碰撞边界，蓝色涂抹线代表大洋俯冲边界;(b-d)喜马拉雅-东印度洋板块会聚边界的深部结构层析成像结果(Replumaz et al. ，2004);(e-f)分别指示扎格罗斯大陆碰撞带及其东部的马库元大洋俯冲带的深部结构层析成像结果(Amaru，2007)，剖面位置如图中浅蓝色和深蓝色直线所示;(g)安纳托利亚地区的深部结构(Gessner et al. ，2013)Fig.6 Tectonics and deep structures of the middle-eastern Tethyan belt(a)topography of the middle-eastern Tethyan belt and its surrounding regions. On the map，red charcoal lines indicate the continental collision boundaries，whereas blue charcoal lines indicate the oceanic subduction boundaries;(b-d)deep structure of the Himalaya and Eastern Indian Oceanic plate convergent regions from seismic tomography (Replumaz et al. ，2014);(e-f)demonstrate the deep structure of Zagros continental collision belts and Makran oceanic subduction zones，respectively，based on the seismic tomography (Amaru，2007);(g)deep structure of the Anatolia region according to the systematic geological and geophysical reviews (Gessner et al. ，2013)

模型所展示的大尺度深部结构特征与中-东特提斯构造带具有很好的相似性，在此我们进行三组对比。首先，对于喜马拉雅大陆碰撞造山带及东部的东印度洋大洋板块俯冲带，层析成像的结果表明，在浅部(例如，200km 深度处)，大陆碰撞带和大洋俯冲带都显示了高地震波速的俯冲板块的存在(如图6b);但是，到一定深度处(如440km)，喜马拉雅带的高速体缺失，意即该处俯冲板块发生断离(如图6c);而到很深的位置处(如1175km)，存在沿走向连续的高速体，解释为先期俯冲的特提斯大洋板块(如图6d)。这种三维空间结构与我们的数值模型相吻合，即大陆碰撞一侧发生俯冲板块的断离，形成一个断离窗口，而大洋俯冲板块具有垂向的连续性。而后，对于扎格罗斯大陆碰撞造山带及其东部的马库元大洋俯冲带，前人的深部探测的结果也揭示大陆碰撞带板块断离的可能存在(如图6e)，但相邻的大洋俯冲带却显示连续的大洋俯冲板块高速体(如图6f)。从而也显示该区域的深部结构和动力学特征与数值模型结果的相似性。最后，在扎格罗斯造山带西部的安纳托利亚地区，也存在深部结构沿走向的巨大差异性(如图6g)。基于我们的数值模型可见，大洋俯冲和大陆碰撞沿走向的空间转换是该地区大地构造特征的主控因素。

我们的数值模型提供了俯冲板块沿走向差异性的一个端元模型，即大陆板块和大洋板块相邻，这也是特提斯构造带的典型状况。如果对此进行拓展，可以联想到，在自然界中，会聚的大陆板块(也包含某些大洋板块)可能沿走向不是均一属性的，因此在碰撞的过程中就会造成深部结构及表层造山过程的差异性，这可以为地质、地球物理观测的解译提供一定的思路。

6 结论

(1)俯冲-碰撞带中，流体-熔体活动可以降低周围岩石的流变强度以及会聚板块之间的耦合作用，并且能够很大程度上促进大陆碰撞带中俯冲板块的断离。

(2)大洋俯冲和大陆碰撞的空间转换数值模型揭示深部结构巨大的沿走向的差异性，大陆碰撞带发生俯冲板块断离，而大洋俯冲板块持续下插。同时，地壳物质发生从陆-陆碰撞带向洋-陆俯冲带的侧向逃逸。

(3)这种三维空间中沿走向的差异性俯冲-碰撞模式与中-东特提斯构造带三个典型地区相吻合，揭示了其深部结构的控制机制;如喜马拉雅碰撞带及东印度洋俯冲带、扎格罗斯碰撞带及马库元俯冲带，以及安纳托利亚陆-陆碰撞和洋-陆俯冲转换带。

致谢 感谢罗纲教授、王世民教授对稿件提出的一系列建设性的意见。

Amaru ML. 2007. Global travel time tomography with 3-D reference models. Geologica Ultraiectina，274:49 -79

Connolly JAD. 2005. Computation of phase equilibria by linear programming:A tool for geodynamic modeling and its application to subduction zone decarbonation. Earth and Planetary Science Letters，236(1 -2):524 -541

Faccenda M. 2014. Water in the slab:A trilogy. Tectonophysics，614:1-30

Gerya TV. 2010. Introduction to Numerical Geodynamic Modelling. New York:Cambridge University Press，1 -345

Gerya TV and Meilick FI. 2011. Geodynamic regimes of subduction under an active margin:Effects of rheological weakening by fluids and melts. Journal of Metamorphic Geology，29(1):7 -31

Gessner K，Gallardo LA，Markwitz V，Ring U and Thomson SN. 2013.What caused the denudation of the Menderes Massif:Review of crustal evolution，lithosphere structure，and dynamic topography in Southwest Turkey. Gondwana Research，24(1):243 -274

Gorczyk W，Willner AP，Gerya TV，Connolly JAD and Burg JP. 2007.Physical controls of magmatic productivity at Pacific-type convergent margins:Numerical modelling. Physics of the Earth and Planetary Interiors，163(1 -4):209 -232

Hatzfeld D and Molnar P. 2010. Comparisons of the kinematics and deep structures of the Zagros and Himalaya and of the Iranian and Tibetan plateaus and geodynamic implications. Review of Geophysics，48(2):RG2005

Katz RF， Spiegelman M and Langmuir CH. 2003. A new parameterisation of hydrous mantle melting. Geochemistry，Geophysics，Geosystems，4(9):1073

Li ZH，Gerya TV and Burg JP. 2010. Influence of tectonic overpressure on P-T paths of HP-UHP rocks in continental collision zones:Thermomechanical modeling. Journal of Metamorphic Geology，28(3):227 -247

Li ZH，Xu ZQ，Gerya TV and Burg JP. 2013. Collision of continental corner from 3-D numerical modeling. Earth and Planetary Science Letters，380:98 -111

Li ZH，Di Leo JF and Ribe NM. 2014. Subduction-induced mantle flow，finite strain and seismic anisotropy:Numerical modeling. Journal of Geophysical Research:Solid Earth，119(6):5052 -5076

Li ZH，Liu MQ and Gerya T. 2015. Material transportation and fluidmelt activity in the subduction channel:Numerical modeling.Science China (Earth Sciences)，58(8):1251 -1268

Peacock SA. 1990. Fluid processes in subduction zones. Science，248(4953):329 -337

Replumaz A，Kárason H，van der Hilst RD，Besse J and Tapponnier P.2004. 4-D evolution of SE Asia’s mantle from geological reconstructions and seismic tomography. Earth and Planetary Science Letters，221(1 -4):103 -115

Wang CY，Flesch L，Silver P，Chang LJ and Chan W. 2008. Evidence for mechanically coupled lithosphere in Central Asia and resulting implications. Geology，36(5):363 -366

Wang Q. 2010. A review of water contents and ductile deformation mechanisms of olivine:Implications for the lithosphere-asthenosphere boundary of continents. Lithos，120(1 -2):30 -41

Xu ZQ，Ji SC，Cai ZH，Zeng LS，Geng QR and Cao H. 2012.Kinematics and dynamics of the Namche Barwa Syntaxis，eastern Himalaya:Constraints from deformation，fabrics and geochronology.Gondwana Research，21(1):19 -36

Yin A. 2010. Cenozoic tectonic evolution of Asia:A preliminary synthesis. Tectonophysics，488(1 -4):293 -325

Zheng YF. 2012. Metamorphic chemical geodynamics in continental subduction zones. Chemical Geology，328:5 -48