离散元方法（Discrete element method）是1971年由Cundall建立理论基础并付诸实践的一种离散元数值模拟技术^[1]，已广泛应用于构造地质研究，如蔡申阳等进行了岩块的应变分析和裂缝预测^[2]、Iwashita等研究了剪切带问题^[3]、Strayer等诸多学者研究了断层及断层相关褶皱问题^[4-8]、张洁等探究了主动底辟盐构造问题^[9]、Steffen等对比分析了用离散元技术和物理模拟方法对继承性断层发育特征 ^[10-13]的研究等。

俯冲作用是汇聚板块边缘的相互地质作用过程，表现为一个上覆板块斜插入另一个板块之下，由于洋壳密度大于陆壳密度，典型的稳定俯冲作用带往往是一块洋壳俯冲到陆壳之下或者另一块洋壳之下，在大陆岩石圈由于碰撞过程而产生地壳变形，而在大洋岩石圈则表现为板块的消亡和被动大陆边缘的形成^[14]。多种因素影响着板块俯冲作用变形，如板块间耦合程度，板块年龄，地幔流体作用等^[15-16]，而俯冲带下切角度的变化显著影响着板块俯冲作用变形特征，且日益引起国内外研究者的关注。近年利用离散元模拟技术研究板块俯冲变形也取得了一定进展：Hardy等用离散元方法研究了具有不同基底摩擦系数对由于俯冲作用所形成的俯冲楔的影响，并与物理实验模拟结果进行了对比分析^[17]。Burbidge等用离散元研究了增生楔和褶皱-逆冲带的发育过程，且设置了不同的颗粒间摩擦系数和颗粒与基底间摩擦系数，并得出了颗粒与基底间的摩擦系数是影响地层变形特征的主要因素^[18]。Naylor等用离散元方法研究了不同边界条件对俯冲变形构造的影响，并与物理砂箱模拟进行了对比分析^[19]。然而对以不同角度进行俯冲所产生的构造变形的变化还没有较好的数值模拟研究。

本文以离散元原理为基础，使用Itasca公司旗下的PFC2D颗粒流数值模拟软件构建二维离散元模型，模拟板块俯冲变形过程，解析不同俯冲角度下地层变形过程中断裂系统的发育和演化特征，揭示俯冲角度对断层发育的影响，为深入了解板块俯冲变形演化历史解释提供参考。

1.   模型构建

本实验着重研究洋壳地层以不同角度向陆壳俯冲时的板块俯冲变形过程，以对比不同变形样式。本实验设置俯冲角度分别为30°、15°、5°的实验模型（图1），以离散元颗粒模拟主体地层，赋予不同微观接触属性模拟洋壳与陆壳地层，两侧及底部以刚性的墙体作为边界，通过底部活动墙体模拟洋壳向陆壳的俯冲作用。

图  1  离散元数值模拟俯冲变形实验装置示意图

Figure  1.  Schematic diagram of subduction deformation experimental device for discrete element numerical simulation

下载: 全尺寸图片 幻灯片

为了保证俯冲作用开始时洋壳相对于陆壳的位置不变，且俯冲带两侧洋壳和陆壳的长度相同，因此，随着洋壳俯冲角度的不同，设置3组实验模型总长度不同，随着俯冲角度变缓模型总长度随之增加，第1组实验模型俯冲角度为30°，总长度为100 m，第2组实验模型俯冲角度为15°，总长度为120 m，第3组实验模型俯冲角度为5°，总长度为200 m。除此之外，3组实验模型具有相同的其他参数设置，最终俯冲量都为25%，陆壳地层厚度都为10 m，洋壳地层厚度都为5 m，都以15%孔隙度、2 600 kg/m³密度随机充填粒径0.125～0.075 m的二维圆盘颗粒，不同粒径颗粒数量服从高斯分布。右侧洋壳基底设有活动墙体，与右侧墙体都以0.5 m/s的速度向左侧陆壳俯冲运动，实验模型左侧陆壳地层均分为10层，每层厚度1 m，中部及右侧洋壳地层均分为5层，每层厚度1 m，以不同颜色加以区分，方便观测俯冲过程中地层构造变形特征。

据前人资料和离散元原理可知，PFC2D软件是基于颗粒单元之间的接触关系而进行模型参数设定的，不同于对应的宏观属性，例如在设置实验参数弹性模量时，其实验设定值小于宏观的弹性模量，实验设定的颗粒单元之间微观摩擦系数大于宏观的摩擦系数，模型的宏观抗张强度与抗剪强度则以设定微观颗粒正向与切向黏结强度来表现^{[11, 20-21]}。参照前人的模型设置及宏、微观参数对应关系，结合物理砂箱模拟实验材料石英砂的力学属性，本文将实验主体陆壳地层弹性模量设为2×10⁷ Pa，洋壳地层的弹性模量设为4×10⁸ Pa（表1），以显著表现洋壳地层刚度大不易变形的特点。此处的摩擦系数指的是在构成板块的离散元颗粒间的摩擦系数，是一个微观的设置参数，一般岩石的摩擦系数都为0.6，且本次实验研究的是俯冲角度对构造变形样式的影响，应该保持每组的摩擦系数不变，因此，地层颗粒间属性均设定为0.6的摩擦系数，使模拟材料形成更加符合库仑剪破裂准则的一系列断层组合形式。

表  1  离散元数值模型参数

Table  1.  Parameters of the discrete element numerical simulation

实验数据	弹性模量/Pa	摩擦系数	法向与切向黏结强度/Pa	速度/（m·s−1）	俯冲角/（°）
实验1	陆壳地层2×107 洋壳地层4×108	0.6	3.6×105	0.5	30
实验2	陆壳地层2×107 洋壳地层4×108	0.6	3.6×105	0.5	15
实验3	陆壳地层2×107 洋壳地层4×108	0.6	3.6×105	0.5	5

下载: 导出CSV 

| 显示表格

2.   实验结果

2.1   30°俯冲角

第1组离散元实验模型俯冲角为30°，左侧和右侧洋壳地层呈水平展布，中部洋壳地层呈倾斜展布，洋壳地层及其边界共同以0.5 m/s的水平速度向左侧陆壳方向匀速运动，陆壳地层受到挤压作用产生相应形变，并跟据断层出现的先后顺序编制其序号。在俯冲量达到5%的早期俯冲变形阶段，随着洋壳地层向左侧运动，带动上覆陆壳地层产生断裂变形，倾向陆壳的断层F1、F2、F3和F4相继形成，随着俯冲量的增大，逐步形成倾向洋壳的反冲断层F5、F6。整体倾向陆壳的断层发育早于F5、F6（图2），且发育方式为前展式，因其右侧的洋壳地层限制了断层的继续发育，所以其前展式特征不明显，随后发育反冲断层，发育方式也为前展式，逐步向陆壳方向发展，且其上发育小型分支断层F7，断距小，主要起到调节陆壳地层挤压量的作用。在俯冲量达到10%的早期挤压变形阶段，反冲断层数量增多，断层以前展式向陆壳方向扩展，形成断层F8，并与先前形成的反冲断层呈近平行关系，随后形成小型分支断层；倾向洋壳断层数量无明显变化，且随着洋壳俯冲活动性较强的影响，断层向深部延伸，且其断距明显增大；整体上俯冲楔宽度及高度增大。在俯冲量达到15%的中期挤压变形阶段，反冲断层继续发育，断层向陆壳方向扩展，形成断层F10，与先前形成的反冲断层呈近平行关系；倾向洋壳断层数量无明显变化，其断距继续增大；整体上其形成俯冲楔宽度及高度继续增大。在俯冲量达到20%的中后期挤压变形阶段，断层数量明显增多，反冲断层继续发育，相继形成断层F10、F11、F12，且这3条断层控制了俯冲楔左侧地层的抬升；倾向洋壳断层数目无明显变化，其断距继续增大；整体上俯冲楔宽度变大，左侧地层继续抬升，其形态发生改变，从之前的左低右高演变为左侧比右侧略高的楔形。在俯冲量达到25%的晚期挤压变形阶段，断层数量无明显变化，反冲断层F10、F11、F12断距增大，这3条断层控制了俯冲楔左侧地层的继续抬升；倾向洋壳断层数目无明显变化，其断距继续增大；整体上俯冲楔宽度变大，左侧地层继续抬升，演变为左侧比右侧高的楔形。俯冲量从0逐渐增加到25%，整体上来讲，在陆壳与洋壳接触部位，先以前展式依次发育倾向陆壳断层，断层断距不断增大，随后以前展式相继发育反冲断层，断层逐渐向陆壳内部发展，小型分支断层发育，当俯冲量达到13%时，第1条控制俯冲楔整体形态的断层F10形成，逐渐增大俯冲量，当俯冲量达到16%时，第2条控制俯冲楔整体形态的断层F11形成，俯冲作用继续进行，当俯冲量达到17%时，第3条控制俯冲楔整体形态的断层F12形成，俯冲楔高度不断增高，宽度不断变大，且其形态逐渐由左低右高演变为左侧比右侧高的楔形。

图  2  30°俯冲角离散元数值模拟实验变形过程解释图

Figure  2.  Interpretation of the deformation process in the discrete element simulation experiment at 30° subduction angle

下载: 全尺寸图片 幻灯片

2.2   15°俯冲角

减小洋壳俯冲角为15°，进行第2组离散元俯冲变形模拟实验，陆壳总长度为80 m，洋壳总长度为60 m，其他条件及参数设置与实验模型1相同（表1），右侧及中部洋壳地层边界共同以0.5 m/s的水平速度自右向左匀速运动。俯冲量从0逐渐增加到25%，其变形过程与30°俯冲角下变形过程相似（图3），但其变形特征有些差异。在陆壳与洋壳接触部位，以前展式依次发育倾向陆壳断层，断层断距随着俯冲作用的进行而不断增大，随后以前展式相继发育反冲断层，断层逐渐向陆壳内部发展，且小型分支断层发育，当俯冲量达到17%时，第1条控制俯冲楔整体形态的断层F14形成，逐渐增大俯冲量，当俯冲量达到22%时，第2条控制俯冲楔整体形态的断层F15形成，继续增大俯冲量，当俯冲量达到25%时，第3条控制俯冲楔整体形态的断层F16形成，由于俯冲作用形成的俯冲楔宽度不断增大，高度不断增高，且其形态逐渐由左低右高演变为左右两侧高度近乎相等的较为规则的楔形。

图  3  15°俯冲角数值模拟实验变形过程解释图

Figure  3.  Interpretation of the deformation process in the discrete element simulation experiment at 15° subduction angle

下载: 全尺寸图片 幻灯片

2.3   5°俯冲角

继续减小洋壳俯冲角为5°，进行第3组离散元俯冲变形模拟实验，陆壳总长度为160 m，洋壳总长度为100 m，其他条件及参数设置与前两组实验模型相同，右侧及中部洋壳地层边界共同以0.5 m/s的水平速度自右向左匀速运动。在俯冲量达到5%的早期俯冲变形阶段，倾向陆壳的F1、F2、F3、F4、F5断层相继形成（图4），随后 F1、F3上发育小型倾向洋壳分支断层F6、F7。在俯冲量达到10%的早中期俯冲变形阶段，断层数量明显增多，距F1、F2、F3、F4、F5很近的下方分别发育有与其平行的断层，这些新发育的断层与其相近的先存断层具有相似性质，所以分别用代号F1、F2、F3、F4、F5表示，相距很近的断层F1、F2、F3、F4、F5分别共同控制的地层断距明显增大，与此同时，反冲断层大量增多，发育断层F8、F9、F10、F11、F12，俯冲楔高度增大。随着俯冲作用的不断进行，在俯冲量达到15%的中期俯冲变形阶段，F1、F2、F3、F4、F5控制的地层断距继续增大，发育倾向洋壳分支断层 F14，俯冲楔高度明显增大。在俯冲量达到20%的中后期挤压变形阶段，反冲断层数量增多，相继形成断层F15、F16，且这两条断层控制了俯冲楔左侧地层的抬升；倾向陆壳断层断距继续增大，俯冲楔左侧地层继续抬升。在俯冲量达到25%的晚期俯冲变形阶段，断层数量无明显变化，反冲断层F15、F16断距明显增大，这两条断层控制了俯冲楔左侧地层的继续抬升，倾向陆壳断层断距继续增大，且断层倾角变陡，整体上断层之间排列更加紧密，俯冲楔左侧地层继续抬升，演变为左右坡度近于相等的楔形。

图  4  5°俯冲角离散元模拟实验变形过程解释图

Figure  4.  Interpretation of the deformation process in the discrete element simulation experiment at 5° subduction angle

下载: 全尺寸图片 幻灯片

3.   模拟结果

通过对比3组实验俯冲量达到25%时俯冲变形结果可以看出，在经历了25%的俯冲量后，实验1、实验2、实验3均发育了板块俯冲变形断裂系统。将3组实验进行对比发现,随着俯冲量的增加，断层数量不断增多，整体断距不断增加；在不同俯冲角度下，俯冲角度越低，形成断层数量越多，地层变形区域越大，形成俯冲楔的高度越高，宽度越大，且不同俯冲角度下形成俯冲楔的形态不同。

根据实验模型中各离散元颗粒最终位移量及位移方向，在PFC软件中可提取累积位移量图（图5），每幅图中图a表示离散元颗粒最终位移量，图b通过矢量箭头显示出伸展变形过程中模型内部颗粒的瞬时位移方向，在PFC中对累积位移矢量图视图窗口进行放大，可清晰观测到陆壳在此部分为箭头向左的细小箭头形式，因颗粒按照箭头的形式时密度太大，因此箭头间有叠加，在图片上显示箭头方向并不明显，所以特意以红色箭头形式标明。色度差异反映了最终累积位移量的大小，累积位移矢量图直观且定量化地表现了实验模型的最终形变特征。

图  5  离散元模拟实验累计位移量图

Figure  5.  Cumulative displacement from the discrete element simulation experiments

下载: 全尺寸图片 幻灯片

通过3组实验的累积位移矢量图对比可以看出，洋壳地层不断向陆壳方向运动因而位移量最大，其次，在洋壳与陆壳接触处位移量较大，向左侧位移量逐渐减小，整体上位移量具有明显分带特征。实验1中在洋壳地层与陆壳地层相接处和左侧陆壳地层俯冲楔边缘处具有明显的位移量突变界面，前者是由于洋壳地层带动着陆壳地层跟随洋壳向俯冲方向运动地层突变，后者突变界面是由于陆壳地层受到挤压向上抬升地层突变，模型总体以两处位移突变界面划分为3个区带，分别为左侧的静止带、中央的位移渐变带和右侧的洋壳大量位移带。实验2与实验3的累积位移矢量图特征基本一致，与实验1相比，也具有两处位移突变界面和3个位移区带，但实验2和实验3的中央位移渐变区域较大，且实验3的中央位移渐变区域最大；两处位移突变界面整体倾角比实验1有逐渐减小的趋势。整体来讲，两处位移量突变界面组合成楔状，且中间具有位移渐变区域，这与实验结果中所反映陆壳地层受洋壳地层俯冲作用所形成的俯冲楔变形断裂系统的特征相符合。

通过检测边界断层两盘对应地层垂直距离获得铅直地层断距，断距随俯冲量的变化能够反映边界断层在俯冲变形过程中活动性的改变。在实验1、实验2和实验3中，F1断层出现的时间最早，随着俯冲量的变化断距也更加明显，在实验过程的各个阶段都有F1断层的发育，相较于其他断层，如F2、F3断层，其出现时间比F1略晚，且其变化特征较弱，因此，F1断距更能反映俯冲楔整体的变化特征。统计F1断距随俯冲量的变化并绘制折线图，可见在3组实验中，F1断距均随俯冲量的增大而不断加大（图6），断距增长速率呈现先快后慢的趋势，断层活动性逐渐变小，且在俯冲作用进行到最后阶段，俯冲量达到25%时，洋壳以30°俯冲角俯冲时形成的F1断距大于以5°俯冲角俯冲时形成的F1断距，15°俯冲角向下俯冲时形成的F1断距最大。据分析可知，断层F1是由洋壳向陆壳俯冲时，带动陆壳地层运动，陆壳地层受到挤压而发生破裂变形所产生。

图  6  3组实验中F1断层断距随俯冲量变化图

Figure  6.  Fault displacement variation of F1 with subduction in three experiments

下载: 全尺寸图片 幻灯片

在俯冲作用前期，F1断距增长速度较快，俯冲作用后期阶段产生了其他破裂变形，分担了F1断层的活动性，因此，F1断距增长速率变缓。相较而言，当俯冲角为30°时，俯冲作用早期阶段断层F1活动性最大，其断距增长速率也最快，但到了俯冲作用晚期阶段，断层活动性减弱，其断距增长速率极为缓慢，虽然早期阶段F1断距在15°俯冲角条件下比30°俯冲角条件下要小，但因在15°条件下时F1断层活动性持续时间较长，因此，到了俯冲作用晚期阶段其累积断距大于在30°俯冲角条件下的累积断距；俯冲角为5°时，因为其俯冲角度较为平缓，洋壳向陆壳进行俯冲作用时，其带动上部陆壳地层产生破裂变形作用强，其他断层分担了F1断层的断距，所以其在折线图上表现出的断距及其增长速率最小，但这并不意味着洋壳以5°俯冲角向陆壳俯冲时其产生的俯冲变形作用最弱，相反，由上文的实验结果图可以看出，其产生的破裂变形作用最为明显，产生断层数量最多。

检测洋壳俯冲时带动上部陆壳地层变形所形成俯冲楔距地表的高度，俯冲楔高度能够反映陆壳地层整体变形情况。图7绘制出实验1、实验2和实验3的俯冲楔高度随俯冲量变化图，可见3组实验俯冲楔高度均随俯冲量的增大而不断增高，在30°和15°俯冲角条件下，俯冲楔高度增长速率较为平缓，而5°俯冲角条件下，俯冲楔高度增长速率较为迅速，且在俯冲作用进行到最后阶段，俯冲量达到25%时，洋壳在15°俯冲角条件下俯冲楔高度大于30°俯冲角条件下俯冲楔高度，洋壳以5°角向下俯冲时形成的俯冲楔高度最大，由此可知俯冲角越小，形成的俯冲楔高度越高。据分析可知，随着俯冲作用的不断进行，陆壳受洋壳挤压程度越来越大，地壳地层变形越高，所形成的俯冲楔高度也就越大，当洋壳以30°高角度进行俯冲时，陆壳地层所产生变形的区域较小，因此，其形成的俯冲楔高度较低；当洋壳以15°较低角度俯冲时，陆壳变形的区域较大，因此，其形成的俯冲楔高度较高；当洋壳以5°低角度俯冲时，陆壳地层产生较大幅度及范围的变形，形成的俯冲楔高度最高。

图  7  3组实验中俯冲楔高度随俯冲量变化图

Figure  7.  Variation of wedge height with subduction in three experiments

下载: 全尺寸图片 幻灯片

由此可见，俯冲角度的变化是影响板块俯冲变形的重要因素，俯冲角度不同，其最终形成的构造样式也不同。总体来讲，洋壳以低角度俯冲时，陆壳地层变形程度较大，其形成的俯冲楔高度较高，地层变形范围较大，形成的断层数量越多，且控制俯冲楔整体形态的反冲断层形成时间较晚；当洋壳以较高角度进行俯冲作用时，陆壳地层变形程度较小，其形成的俯冲楔高度较低，地层变形范围较小，形成的断层数量较少，且控制俯冲楔整体形态的反冲断层形成时间较早。

4.   讨论

本文参考前人研究成果，选取了菲律宾板块向欧亚板块俯冲部位的几条地震剖面图（图8）^[22-23]，并选择具有代表性的区域进行了相应解释，观察其断层分布特征。在靠近菲律宾板块俯冲部位发育多条逆冲断层，陆壳内侧分布有多条倾向洋壳的断层，陆壳经由俯冲作用而形成局部隆起，整体形态为左侧比右侧高的楔形。将其与离散元俯冲变形模拟实验对比（图9），发现30°俯冲角实验模拟结果与之相似程度最高，具有相同的俯冲变形特征，俯冲楔形态最为相似，参照前人地震资料统计结果，菲律宾板块向欧亚板块俯冲角度约为30°^[24-26]，这与本文实验模拟结果吻合，由此可知，本次实验模拟结果可以为探讨菲律宾板块俯冲演化史以及对东海盆地的成因分析提供参考。

图  8  东海陆架盆地及地震剖面地理位置平面图（a）及地震剖面反射图（b, c, d）b图左下方为P波速度剖面图（据文献[22-23]修改）

Figure  8.  Geographical map of East China Sea Shelf Basin (a) with seismic profile (b), (c) and (d) The lower left of (b) shows the P wave velocity profile (after references[22-23])

下载: 全尺寸图片 幻灯片

图  9  实验结果与实际地震剖面对比

a. 图8中经过构造解释的放大图，b. 30°俯冲角离散元数值模拟实验结果构造图（据文献[22]修改）

Figure  9.  Comparison of the experimental results with seismic profile

a. Enlarged view in Fig.8 after structural interpretation, b. Result of 30° subduction angle simulation experiment (after reference[22])

下载: 全尺寸图片 幻灯片

北印度洋的马克兰增生楔是阿拉伯板块俯冲到欧亚板块之下形成的，其低俯冲角形成了宽广的增生楔，马克兰增生楔主要由大量逆冲断层组成^[27-29]。现将其构造特征与本次实验进行对比可知（图10），5°俯冲角的实验与其具有较高的相似度，具有相似的俯冲变形特征，俯冲楔形态也较为相似，断层形态整体较缓，参照前人地震资料统计结果，阿拉伯板块向欧亚板块俯冲角度小于5°，这与本文实验模拟结果较为吻合，由此可知，本次实验模拟结果可以为探讨阿拉伯板块俯冲演化史提供参考。

图  10  实验结果与马克兰增生楔对比

a. 马克兰增生楔，b. 5°俯冲角离散元数值模拟实验结果构造图（据Grando等改, 2007）

Figure  10.  Comparison of the experimental results with Makran's accretion wedge

a. Maklan accretion wedge, b. Result of 30° subduction angle simulation experiment (from Grando et al., 2007)

下载: 全尺寸图片 幻灯片

东海大陆架盆地位于欧亚板块的西南部，与菲律宾板块相邻，是经过多期俯冲作用形成的（图11）。中生代时期太平洋板块向东亚大陆边缘发生低角度俯冲，使东海陆架盆地发生内陆拗陷，形成挤压拗陷型盆地；古近纪时期，由于菲律宾板块高角度俯冲作用，导致琉球岛弧前缘增生楔发育^[30-37]。经过以上的实例对比可知，本次实验结果具有一定的指导意义，因此，可用15°俯冲角的模拟结果与太平洋板块向东亚大陆边缘发生俯冲时形成的东海陆架盆地中生代增生楔部分做对比，以验证中生代时期太平洋板块的俯冲角度。本文因资料有限，未能截取合适的剖面以作比对，期待下一步有更多的资料开放以支持研究。

图  11  东海陆架盆地南部演化模式图（据文献[30]修改）

Figure  11.  Evolution model of the East China Sea Shelf Basin (after reference[30])

下载: 全尺寸图片 幻灯片

5.   结论

（1）板块俯冲变形特征与俯冲角度有关，俯冲角度不同，其最终形成的变形样式也不同。在不同俯冲角度下，其相同位置形成的断层断距也不同。

（2）随着俯冲角度由30°、15°到5°的减缓，断层数量不断增多，且断层扩展的水平距离不断增大。俯冲楔高度不断增加，陆壳地层变形幅度越大，最终形成俯冲楔形态类型不同。

（3）控制俯冲楔整体形态的反冲断层形成时间随着俯冲角度的减缓而变晚。

（4）实验模拟结果与实例进行对比发现，两者表现出相似的构造特征，具有良好的对应性，可依据不同的俯冲楔样式推测俯冲角度，对板块俯冲作用具有指导意义。