智利三联点南部含洋脊海洋板块俯冲过程的数值模拟

2022-01-25原一哲郭长升胡才博魏东平

地球物理学报 2022年1期

原一哲，郭长升，胡才博*，魏东平

1 中国科学院大学地球与行星科学学院，北京 100049 2 中国科学院计算地球动力学重点实验室，北京 100049

0 引言

三联点是三个板块边界的交点(Mckenzie and Morgan,1969).按板块边界的类型(洋脊R、转换断层F、海沟T)与几何形状划分为十六种不同的类型(王振山和魏东平，2018).智利南部扩张洋脊的俯冲是影响南美板块边界的主要构造活动之一，也是地球动力学领域的重要研究内容(Atwater,1970;Kay et al.,1993;McCrory et al.,2009).

智利三联点(CTJ)位于南极洲板块、纳兹卡板块与南美板块的交界处，由南极洲-纳兹卡板块之间的智利洋脊俯冲至智利海沟形成.随着纳兹卡板块的东向俯冲，纳兹卡板块范围逐渐变小，智利三联点(CTJ)同时向北移动.大约14 Ma前至今，智利三联点从54°S不断向北移动至46°30′S(Cande and Leslie,1986)(如图1).以智利三联点对应纬度(46°30′S)为南北分界线，北部区域，纳兹卡板块以7.8 cm·a-1的速率快速俯冲于南美板块下方，南部区域，南极洲板块以2 cm·a-1速率缓慢俯冲于南美板块之下(DeMets et al.,1994).由于智利三联点北部区域智利洋脊段并未俯冲至南美板块之下，而南部区域智利洋脊段SCR-1则在距今6 Ma俯冲至南美板块之下，SCR0洋脊段则在距今3 Ma俯冲至南美板块之下(如图1).根据智利三联点北部贝尼奥夫地震带的分布，智利三联点北部洋脊俯冲至南美板块的倾角为30°(Cahill and Isacks,1992)，而智利三联点南部尚未观测到明显的贝尼奥夫地震带(Breitsprecher and Thorkelson,2009).

图1 智利三联点地形图(修改自Agurto-Detzel et al.,2014;Lagabrielle et al.,2015)NZ：纳兹卡板块；AN：南极洲板块；SA：南美洲板块；CTJ：智利三联点；SCR:洋脊；红色五角星:5.9 Ma智利三联点位置；黑色五角星:2.2 Ma智利三联点位置(Maksymowicz et al.,2012)；白色线段表示海沟位置；黄色线段代表转换断层带(实线)和预测转换断层带(虚线)；红色实线代表智利洋中脊位置的扩张中心(智利三联点南部包括SCR-1和SCR0；北部包括SCR1、SCR2和SCR3)；红色三角形代表第四纪火山群位置.Fig.1 Topographic map of the Chile triple junction (Modified from Agurto-Detzel et al.,2014;Lagabrielle et al.,2015)NZ:Nazca plate,AN:Antarctic plate;SA:South American plate;CTJ:Chile triple junction;SCR:Ocean Ridge;Red star:Chile triple point location at 5.9 Ma;Black star:Chile triple point location at 2.2 Ma (Maksymowicz et al.,2012);White lines denote the trench;The yellow lines denote the fracture zone (solid line)and the predicted fracture zone (dashed line);Red segments show the spreading centers of the Chile ridge (The southern part of the CTJ includes SCR-1 and SCR0;The northern part includes SCR1,SCR2,and SCR3);Red tringles denote Quaternary volcanoes.

含有洋脊的洋-陆俯冲与一般的相比最显著的特征是俯冲洋壳年龄小、厚度薄、俯冲区域温度梯度高(Lagabrielle et al.,1999)，在靠近海沟轴线附近伴随有岩浆活动(Bourgois et al.,1996).现代观测数据表明，智利三联点南部区域呈现明显的热流正异常现象(Hamza et al.,2005)，明显高于南美洲的平均热流值52.7 mW·m-2(Sclater et al.,1980).

板片窗指板片在俯冲过程中由洋脊扩张、板片撕裂或板片断离形成的间隙.三联点附近地质构造类型存在明显差异，其中典型的案例就是南美板块中形成的巴塔哥尼亚板片窗(Bourgois and Michaud,2002).在智利巴塔哥尼亚区域和阿根廷巴塔哥尼亚弧后地区发现的始新世高原玄武岩(Demant and Morata,1996;Parada et al.,2001)已经被证明是洋脊-海沟俯冲的产物(Kay et al.,2002;Kay and Mpodozis,2002).研究洋脊-海沟俯冲机制对于了解三联点附近地质构造类型及板片窗的形成具有重要意义.

本文通过建立二维洋脊-海沟俯冲模型，采用有限差分法和粒子追踪技术研究岩石层尺度内地球流体介质在热、重力及构造作用下的流动、变形及温度场变化，定量模拟智利洋脊俯冲至南美板块的动力学过程.针对智利三联点附近区域洋脊俯冲的角度和速率不同，设置对照实验，研究不同的俯冲角度和速率条件下，洋脊俯冲对上覆岩石层温度结构及地质演化过程的影响.研究洋脊俯冲区域地表热异常机制对于了解地壳与上地幔的热结构及评定区域地热资源潜力等均具有重要意义.

1 数值模拟方法

1.1 控制方程

本文采用地球动力学程序“I2VIS”进行数值模拟(Gerya and Yuen,2003a)，控制方程分别为：斯托克斯流体动力学方程、不可压缩质量守恒方程、温度演化的能量守恒方程，可用于板块俯冲、碰撞、拆沉、折返与造山等复杂地质过程的时空演化(李忠海和许志琴，2015).

有效黏度ηeff存在两种形式，一种是韧性黏度，一种是塑性黏度.其中，韧性黏度计算如下：

(1)

塑性部分采用Drucker-Prager屈服准则 (Ranalli,1995)：

σyield=C0+Psin(φeff),

(2)

其中，

sin(φeff)=sin(φ)(1-λ).

(3)

定义塑性黏度如下：

(4)

式中，σyield为屈服应力；C0是P=0条件下的岩石剩余强度；φ是内摩擦角；λ是孔隙流体系数；φeff是有效内摩擦角，与孔隙流体系数和内摩擦角相关.

选取ηductile和ηplastic的较小者作为实际的有效黏度参与计算(Ranalli,1995;Gerya and Yuen,2003a):

ηeff=min (ηductile,ηplastic).

(5)

模型计算中考虑了地壳岩石的部分熔融行为(Gerya and Yuen,2003b;Burg and Gerya,2005;Gerya and Burg,2007)，根据实验岩石学的约束条件，将温度与岩石部分熔融体积比例之间的关系看作近似的线性关系.

2 洋脊俯冲模型的模型设置

根据智利三联点北部贝尼奥夫地震带的分布，智利洋脊俯冲至南美板块的倾角约为30°(Cahill and Isacks,1992)，而智利三联点南部没有观测到贝尼奥夫地震带，智利三联点南部的俯冲角度大致为10°～13°，但缺乏相关观测证据(Breitsprecher and Thorkelson,2009).本文以智利三联点南部的浅俯冲地质背景为基础，分别建立了初始俯冲角度分别为15°和30°的数值模型.同时，由于俯冲区域的速度边界条件是控制大洋俯冲模式的重要影响因素之一(Jarrard,1986;Lallemand et al.,2005).因此我们设计两组不同速率的模拟实验(3.5 cm·a-1，6.5 cm·a-1)讨论俯冲速率的影响，各模型的俯冲角度和速率如表1所示.

表1 模型对应的俯冲角度和速率Table 1 The subduction angle and speed used in the model

4个模型的初始宽度为1000 km，深度为300 km(如图2)，使用不规则的有限差分网格对模型进行离散化，俯冲带区域加密，网格设计为1 km×1 km单元，模型其余区域采用10 km×10 km单元.大陆板块顶部长度为400 km (蓝色区域)，大洋板块为600 km (红色区域)，俯冲带倾角为30°.大陆板块总厚度为120 km，分为三部分，包括1 km厚的沉积层，20 km 厚的上地壳，15 km厚的下地壳和余下的岩石层地幔部分.大洋板块总厚度为75 km，其中包括2 km 厚的沉积层，7 km厚的洋壳和余下的岩石层地幔部分.模型宽度500 km处为洋脊构造，洋脊的位置和几何参考前人结果(Hamblin and Christiansen,2003).不同的岩石类型和性质见表2.针对俯冲带区域流变性质相关参数的影响，我们主要考虑上覆大陆岩石层的流变性质，尤其是陆壳的流变强度，以及岩石圈地幔的流变性能，流变参数主要参考前人的文献(Ranalli,1995；李忠海和许志琴，2015).上覆大陆地壳的初始黏度与软流层黏度降低均导致两板块俯冲界面倾角增大(皇甫鹏鹏，2016).

表2 模型采用的相关参数(参考Ranalli,1995；Li et al.,2010;李忠海和许志琴，2015)Table 2 Related parameters used in the model (Refer to Ranalli,1995;Li et al.,2010;Li and Xu,2015)

图2 初始模型物质图.模型中的颜色代表不同的岩石组成.本模型主要分为8个分区：1.伪黏性空气层与沉积层；2.上地壳；3.下地壳；4.洋壳；5.洋脊板块；6.岩石层地幔；7.软流层地幔；8.软弱带.水化的和部分熔融的岩石类型在该初始模型中不存在，随着模型的演化而产生.Fig.2 Initial model material diagram.The colors in the model denotes different rock compositionsThis model is mainly divided into 8 partitions:1.pseudo-sticky air and sediment;2.upper continental crust;3.lower continental crust;4.oceanic crust;5.ridge plate;6.lithospheric mantle;7.asthenospheric mantle;8.weak zone.The hydrated and partially molten rocks are not shown in the initial model,but appear during the evolution of the model.

模型底部边界为渗透性边界 (Burg and Gerya,2005;Ueda et al.,2008;Li et al.,2010).模型在俯冲板块的远端施加俯冲速率(Vx)边界条件，在计算过程中保持恒定，其他边界采用自由滑动的边界条件.相当于静止参照物的绝对运动速率.模型自初始条件开始,为了更真实反映俯冲碰撞后的地表变化,在地壳表面和模型顶界面之间设计了一层相对高黏度的伪空气层，其与上地壳的接触面被用来模拟模型计算过程中的地貌起伏变化，且该地貌起伏面模拟包含近似的地表剥蚀和沉积作用(Gerya and Yuen,2003a;Burg and Gerya,2005).伪空气层的初始厚度在大陆岩石层之上为 7 km，在大洋地壳之上为8 km，大陆岩石层的初始高程比大洋高500 m.由于沉积层太薄，不便于单独的颜色区分，因此将其与伪空气层标注相同颜色.

模型顶部温度固定为0 ℃，莫霍面温度为650 ℃，岩石层中采用线性温度梯度0.5 ℃/km，岩石层底部温度为1400 ℃(Tetreault and Buiter,2012)，岩石层温度采用线性温度分布，如图3a所示.两侧边界的水平方向温度梯度为零(即零热流).底部边界采用的是外部边界固定温度条件，即在模型底边界下假设一个固定的地幔温度(李忠海和许志琴，2015；皇甫鹏鹏，2016).初始模型的密度设置如图3b.

图3 初始模型示意图(a)初始模型温度图；(b)初始模型密度图.Fig.3 Schematic diagram of the initial model(a)Initial model temperature map;(b)Initial model density map.

3 模拟结果

根据海洋板块俯冲速率和初始俯冲角度不同设计了四个模型，进行对比分析.

3.1 低俯冲角度低俯冲速率模型(模型1)

模型1的初始俯冲角度为15°，初始俯冲速率为3.5 cm·a-1(图4).模型运行至6 Ma时，洋脊中心俯冲至地幔60 km左右深度，并且向大陆方向运移200 km(如图4c).温度模拟结果显示，含有洋脊的海洋板块的俯冲使得俯冲区域上覆岩石层温度升高约 400～800 ℃，上覆岩层温度明显升高区域主要集中在模型水平方向220～350 km区域，由于模型设置的含洋脊的大洋板块俯冲开始时，洋脊中心位置位于模型水平方向500 km处，海沟位置位于模型水平方向400 km处，即洋脊俯冲 6 Ma后，上覆岩层温度异常区域位于220～350 km，距离海沟50～180 km，温度异常区水平长度约为130 km(如图5a—5c所示).而通过模拟洋脊俯冲进入海沟2 Ma后，在地下20 km深处温度可达到600 ℃(DeLong et al.,1979).模型运行至4 Ma时洋脊才开始俯冲到海沟之下，俯冲至6 Ma时地表下20 km深处温度为630 ℃，与DeLong等(1979)的研究结果相近，如图5c所示.含洋脊板块俯冲过程中受到上覆板块的刮擦效应及大陆板块逆冲作用的影响，岩石层之间剪切热大大增加，同时伴随着温度较高的洋脊俯冲至大陆板块之下与上覆岩层产生热交换，俯冲区域上覆岩层产生明显的热物质侵蚀上涌现象，对应俯冲区域密度升高200～300 kg·m-3(如图6所示).

图4 初始俯冲角度15°，俯冲速率3.5 cm·a-1模型物质演化结果Fig.4 Material evolution results of the model with initial subduction angle of 15° and subduction rate of 3.5 cm·a-1

图5 初始俯冲角度15°，俯冲速率3.5 cm·a-1模型温度演化结果Fig.5 Temperature evolution results of the model with initial subduction angle of 15° and subduction rate of 3.5 cm·a-1

图6 初始俯冲角度15°，俯冲速率3.5 cm·a-1模型密度演化结果Fig.6 Density evolution results of the model with initial subduction angle of 15° and subduction rate of 3.5 cm·a-1

3.2 低俯冲角度高俯冲速率模型(模型2)

模型2的初始俯冲角度为15°，初始俯冲速率为6.5 cm·a-1.模型结果如图7所示，与模型1相比，在相同的初始俯冲角度条件下(15°)，相同的模型演化时间内(6 Ma)，随着俯冲速率的增大，洋脊中心俯冲至地幔110 km左右深度，比模型1深50 km，并且向大陆方向运移350 km，比模型1远150 km.这种现象主要是由于初始俯冲速率较大时俯冲板块在俯冲区域浅层停留时间相对较短，弱化熔融作用较弱，表现为板块强度较大，更易克服板块前端挤压应力，相对更容易俯冲至地幔深处，向大陆方向运移的距离也相对较长.同时，对应的温度场显示(如图7b)，含洋脊的板块俯冲使得俯冲区域上覆岩石层温度结构异常区域位于模型210～330 km范围，距离海沟70～190 km，温度异常区水平长度约120 km.对应俯冲区域密度升高200～300 kg·m-3左右(如图7c).

图7 初始俯冲角度15°，俯冲速率6.5 cm·a-1模型6 Ma演化结果(a)模型物质图；(b)模型温度图；(c)模型密度图.Fig.7 Model evolution results with initial subduction angle of 15° and subduction rate of 6.5 cm·a-1 at 6 Ma(a)Model material map;(b)Model temperature map;(c)Model density map.

3.3 高俯冲角度低俯冲速率模型(模型3)

模型3的初始俯冲角度为30°，初始俯冲速率为3.5 cm·a-1，模型结果如图8所示.相同的俯冲速率条件下，随着俯冲角度的增大，含洋脊的海洋板块俯冲过程中，在相同的模型演化时间内(6 Ma)，洋脊中心俯冲的深度比模型1深20 km，同时，洋脊中心向大陆方向运移的距离减小10 km，俯冲区域上覆岩层的温度异常区域离海沟越近，其水平长度缩小70 km.具体而言，当俯冲速率为3.5 cm·a-1，初始俯冲角度为30°的洋脊板块俯冲至6 Ma时，洋脊中心俯冲至地幔80 km左右深度，并且向大陆方向运移190 km.对应的温度场显示(如图8b)，含洋脊的板块在6 Ma的俯冲过程中使得俯冲区域上覆岩石层温度结构异常区域范围位于模型水平方向 260～320 km，即距离海沟 80～140 km区域范围，温度异常区水平长度约60 km.

图8 初始俯冲角度30°，俯冲速率3.5 cm·a-1模型6 Ma演化结果(a)模型物质图；(b)模型温度图；(c)模型密度图.Fig.8 Model evolution results with initial subduction angle of 30° and subduction rate of 3.5 cm·a-1 at 6 Ma(a)Model material map;(b)Model temperature map;(c)Model density map.

3.4 高俯冲角度高俯冲速率模型(模型4)

模型4的初始俯冲角度为30°，初始俯冲速率为6.5 cm·a-1，当洋脊板块俯冲至6 Ma时，洋脊中心俯冲至地幔 120 km左右深度，并且向大陆方向运移290 km(图9a).对应的温度场如图9b，含洋脊的板块在6 Ma的俯冲过程中使得俯冲区域上覆岩石层温度结构异常区域位于模型水平方向280～310 km，距离海沟 90～120 km.在相同速率(6.5 cm·a-1)条件下，相比低角度(15°)俯冲，温度异常区域中心向海沟方向迁移了25 km左右，温度异常区水平长度缩减90 km.对应俯冲区域密度升高200～300 kg·m-3左右(如图9c).

图9 初始俯冲角度30°，俯冲速率6.5 cm·a-1模型6 Ma演化结果(a)模型物质图；(b)模型温度图；(c)模型密度图.Fig.9 Model evolution results with initial subduction angle of 30° and subduction rate of 6.5 cm·a-1 at 6 Ma(a)Model material map;(b)Model temperature map;(c)Model density map.

四个模型显示上覆岩层不同的温度异常区和洋脊中心俯冲深度(如表3).含洋脊的海洋板块俯冲过程中，在低初始俯冲角度条件下(15°)，俯冲速率越大，洋脊中心在相同的时间俯冲至地幔的深度越深，洋脊中心向大陆方向运移的距离越大，俯冲区域上覆岩层的温度异常区范围越大，温度异常区域中心位置距离海沟越远(模型1和2).在高初始俯冲角度条件下(30°)，俯冲速率越大，洋脊中心在相同的时间俯冲至地幔的深度越深，洋脊中心向大陆方向运移的距离就越大，俯冲区域上覆岩层的温度异常区范围越大，温度异常区域中心位置与海沟的距离越近(模型3和4).

表3 不同俯冲角度和速率的洋脊俯冲模型结果(单位：km)Table 3 Results of ocean ridge subduction models with different subduction angles and speeds (unit：km)

3.5 含洋脊的海洋板块俯冲区域上覆岩层近地表热流值变化

大地热流是表征地球内部热活动的重要物理参量，研究地表热流值的变化对于了解含洋脊的海洋板块俯冲区域上覆岩层热结构及地球深部地质变化具有重要意义.在一维稳态条件下，热流值可通过地温梯度与岩石热导率的乘积来表述(胡圣标等，2001).由于岩石热导率具有明显的各向异性，并且受压力、温度、岩石成分和结构影响很大，本文取上地壳岩石平均热导率2.5 W·(m·K)-1(Furlong and Chapman,2013)进行计算，由于模型设计的上地壳温度梯度为30 ℃/km，因此热流值高于62.5 mW·m-2的区域在模型中属于热流正异常区域，由于洋脊俯冲区域整体热流值偏高，设定热流值超过70 mW·m-2为热流明显异常区域，研究不同的初始俯冲角度，不同俯冲速率条件下，热流最大值和明显异常区域分布范围.在上述洋脊俯冲模型中水平方向0～400 km范围，地表下1 km深处水平方向每间隔20 km设置一个监测点位，通过记录含洋脊的海洋板块俯冲至6 Ma时地表下1 km处的温度值，和大陆地表温度对比得到温度梯度，与岩石热导率的乘积就是地表下1 km处的热流值，结果如图10所示.

图10 不同俯冲角度和俯冲速率模型在模型运行至 6 Ma时对应的近地表热流值Fig.10 Comparison of the near-surface heat flow values corresponding to different subduction angles and subduction rate models at 6 Ma

图10结果表明：含洋脊的大洋板块俯冲至6 Ma时，初始俯冲角度对近地表热流明显异常区域的范围以及最大值的大小和位置均具有重要的影响；而俯冲速率仅在低初始俯冲角度(15°)俯冲时，对近地表热流明显异常区域的范围以及最大值的大小和位置有较大影响，在高初始俯冲角度(30°)俯冲时，影响不大.具体而言，在低初始俯冲角度(15°)、高速俯冲时(6.5 cm·a-1)，近地表热流值明显异常区域范围最大且热流值也最高，最大热流值达到 81.8 mW·m-2,热流明显异常区域位于距离海沟 70～190 km范围(模型 210～330 km)，热流最大值位于距离海沟120 km(模型280 km)处，热流异常区域的均值约为75.3 mW·m-2；低初始俯冲角度(15°)、低速俯冲时(3.5 cm·a-1)，近地表热流值明显异常区域范围和热流最大值次之，最大热流值达到76.5 mW·m-2，热流明显异常区域范围位于距离海沟 50～180 km范围，热流最大值位于距离海沟120 km处，热流异常区域的均值约为 72.4 mW·m-2；而当初始俯冲角度为30°时，无论是高速还是低速俯冲，热流值明显异常的范围和最大值均大幅减小，且俯冲速率对于热流值影响不大.原因可能是俯冲角度增大后，洋脊中心随着洋脊板块俯冲的深度加大，从而对于近地表的温度结构扰动减弱，导致近地表的热流值明显异常区域范围相比低角度俯冲减小，最大热流值和热流异常区域的均值也相应减小.

4 讨论

4.1 模拟结果与南美洲实测热流值对比

在智利三联点向北迁移中，巴塔哥尼亚下方的软流层窗口打开，引起区域地幔对流扰动，造成局部地区的地壳减薄、火山作用加强(Kay et al.,1993;Gorring et al.,1997,2003;D′Orazio et al.,2000).因此，智利三联点附近地表的热流值高于南美板块其他地区，表现为正异常，说明该区域地幔深部构造活动强烈.受人为活动因素和地形环境的影响，南美地表热流数据测量点位分布密度差异巨大，数据测量方法种类多样，数据来源不一(Birch,1954；Carvalho and Vacquier,1977；Swanberg and Morgan,1978；Santos et al.,1986；Burkhardt et al.,1989).位于南纬46°30′的智利三联点以南的热流测量点位稀疏，难以和数值模拟结果直接对比，故本文与Hamza等(2005)采用最小二乘法获得插值多项式进行拟合的数据进行了对比.Hamza等通过数据分析结果表明，巴塔哥尼亚高原的平均热流值74±13 mW·m-2，高于南美洲大陆的平均热流.巴塔哥尼亚位于智利洋脊段俯冲形成的板片窗上方，对应洋脊俯冲模型的热流值异常区域.智利三联点南部区域属于低角度低速率俯冲，图10分析结果表明，热流明显异常区域位于海沟前方50～200 km内.当洋脊低角度俯冲时(15°)，热流异常区域平均热流为72.4 mW·m-2(3.5cm·a-1),75.3 mW·m-2(6.5 cm·a-1)，与Hamza等的分析结果较为吻合.

4.2 智利洋脊俯冲的地貌响应

随着智利三联点向北迁移，巴塔哥尼亚南部以下板片窗的开放，抵消了俯冲带上方大陆板块的向下偏转，导致部分地区隆起(Guillaume et al.,2009,2010).图11为俯冲角度为15°、俯冲速率为3.5 cm·a-1条件下地形随时间的演化规律.在洋脊俯冲过程中，海沟前方400 km内大陆岩石层高程均有所增大.初始模型中，大陆岩石层的高程比大洋岩石层高500 m；当俯冲进行至1 Ma左右，地表高程剧烈抬升；当俯冲进行至6 Ma，海沟前方400 km内大陆岩石层的高程为1300～2800 m，与初始高程相比抬升了800～2300 m，最大值位于海沟前方320 km处.图1中，SCR-1段洋脊所在位置高程在1000～2000 m左右，据此我们推断，位于南美洲智利三联点以南距离海岸线数百千米范围内的高山隆起及智利冰川的形成均与洋脊俯冲有很大关系.

图11 俯冲角度15°，俯冲速率3.5 cm·a-1模型对应的地貌图Fig.11 The geomorphological map corresponding to the model with the subduction angle of 15° and subduction speed of 3.5 cm·a-1