板块的绝对或相对运动控制着板块边界的断裂带形成和应力场特征.印度洋板块与欧亚板块的碰撞-挤压作用，形成了喜马拉雅弧形造山带，其东部弧顶沿N—NE方向插入青藏高原的东北缘(Molnar and Tapponnier, 1975；England，1986；England and Houseman, 1988；Harrison et al., 1992；付碧宏等，2008)，造成了巴颜喀拉块体和龙门山断裂系深、浅部构造强烈活动和变形，而高原地壳、地幔物质以大型走滑断裂为通道边界向E—ES方向运移(滕吉文等，2014)，导致青藏高原持续隆升，致使周缘发生激烈的地壳形变和频繁的大地震活动等构造运动，并控制了中国西部乃至其以北的广大地区(Van Der Rob et al., 1997；Borwn et al., 2002)，如图 1所示.根据火山岩化学全分析资料和GPS监测计算研究现代板块边界的运动速度，印支、燕山与四川期板内平均缩短速度均在5.4~5.6 cm·a^-1，喜马拉雅期在中国西部平均缩短速度达5.2 cm·a^-1，新构造期仅为2.5~2.8 cm·a^-1，喜马拉雅地区表面运动速度都呈现N—NE方向运动(Sugisaki，1976；万天丰，1993；万天丰和任之鹤，1999；叶正仁和王建，2004；张培震，2008).

图 1

Fig. 1

图 1 龙门山活动断裂与地震断层空间分布图(据付碧宏等, 2008) Fig. 1 Spatial distribution map of the active faults and seismic fault in the Longmenshan fault zone (Fu et al., 2008)

板块之间的相互运动，必然造成板块岩体应力的变化，通过“世界应力图”发现全球大部分地区的最大水平主应力方向与板块绝对运动迹线保持较好的一致性，反映出构造应力与板块运动的关系密切(谢富仁等，2003).沿着喜马拉雅地区的印度洋板块和欧亚板块的碰撞运动产生强烈的挤压构造应力(徐纪人等，2008)，当岩体应力超过其强度极限时必然进入塑性状态(徐芝纶，1980；沈明荣和陈建峰，2006).以上研究表明，地壳板块之间存在相互运动，且运动方向与最大主应力方向具有较好的相关性，岩体的主应力状态决定着岩体的弹塑性状态.但板块塑性区发育的几何特征与主应力之间有何关系？塑性区的发育与断层的形成有何关系？这些问题仍需要我们深入研究.龙门山断裂带是在一定深处收敛成为一条NE方向的深部剪切断裂带(滕吉文等，2014)，本文拟以板块运动为研究的基本条件，依据龙门山附近区域的构造运动特性以及深部岩体力学特性，采用FLAC模拟软件模拟龙门山断裂形成过程和塑性区分布的几何特征，研究区域构造运动对于龙门山断裂带形成及地应力状态的影响.

龙门山地处青藏高原东缘与四川盆地的交接部位, 山-盆构造地貌反差强烈, 是青藏高原周边地形梯度变化最大的地区, 在不到50 km的范围内地形高差达约4000 m(邓起东等，1994; 苏生瑞等，2002；付碧宏等，2008；颜照坤等，2014)，如图 2所示.龙门山断裂带主要由大致平行的3条逆冲、逆冲兼走滑性质的断裂组成，分别为：汶川—茂县断裂F1、映秀—北川断裂F2和灌县—安县断裂F3，该断裂带长500多公里，宽30~50 km，大致呈北东—南西向延伸构造(邓起东等, 1994; Burchfiel et al, 1995; 陈国光等, 2007).

图 2

Fig. 2

图 2 龙门山断裂地形地貌剖面图(颜照坤等，2014) Fig. 2 Topographic profile of the Longmenshan fault zone (Yan et al., 2014)

青藏高原Moho面的深度为60~70 km, 四川盆地Moho面的深度为40 km左右，高原内部物质在推挤压力与重力的作用下向四周运移, 出现了高原物质东流现象(曾融生等, 1992; Royden et al., 1997; 朱守彪和张培震，2009)，具有流变性质的下地壳速度V₂大于上地壳速度V₁(朱守彪和张培震，2009)，如图 3所示; 根据GPS观测结果显示, 龙门山断裂带附近地区位移速度差约为2 mm·a^-1(赖锡安等，2004；王连捷等，2009).

图 3

Fig. 3

图 3 龙门山断裂带运动模型示意(朱守彪等，2009) Fig. 3 Motion model of the Longmenshan fault zone (Zhu et al., 2009)

本文选用FLAC^3D5.0数值模拟软件，它是基于拉格朗日差分法的一种显式有限差分程序(Fast Lagrangian Analysis of Continua)，内置了摩尔-库仑准则，适宜于处理大尺度、大变形工程和地质问题(杨铭键等，2012；王涛等，2015).

计算机数值模拟分析模型的主要平面与龙门山断裂带走向垂直，几何尺寸如图 4所示.走向长度为160 km，包含马尔康至四川盆地的部分区域；高度为40 km，与龙门山断裂带下部的Moho面大致相同；厚度取1 km.

图 4

Fig. 4

图 4 龙门山断裂带数值计算模型 Fig. 4 Numerical calculation model of Longmenshan fault

模型的边界条件和加载方式如图 4所示.上部为自由边界，地形的起伏简化为折线；底部大体与Moho面的深度相同，简化为位移边界条件，在x方向可以运动，z方向为固定支座，即w=0；右侧是板块运动相对静止的边界，简化为在z方向可以运动，其他方向为固定铰支座，即u=v=0.在模型左侧边界施加位移前，对整个模型范围内的岩体施加重力应力场，使模型具备初始应力环境.

模型左侧的移动量按加载步逐步施加，即每个加载步施加1a的移动量.模型的左侧边界相对于右侧可以移动，移动速度取两个边界实际移动速度的差值；考虑到地表与Moho面处的移动速度的不同，在左侧边界上的水平速度根据深度的增加成线性增长，模型上部加载的正向水平速度2 mm·a^-1，底部加载的正向水平速度4 mm·a^-1.

为研究模型不同位置的主应力变化情况，在模型的水平位置中部(80 km)处，地表下1 km、11 km、21 km、31 km位置处设置最大最小主应力追踪测点.

岩石破坏选用摩尔-库仑准则，其基本的物理力学参数包括弹性模量、泊松比、内聚力、内摩擦角、剪胀角和抗拉强度.

迄今为止，对地壳深部岩石力学性质的研究不是很多，仅限于几公里范围内(苏生瑞等，2002).一般认为，随着深度的增加，岩体更加致密，岩体的密度和弹性模量等都会随深度的增大而增大(蔡美峰和王双红，1997；秦向辉等，2012；柳畅等，2014; 熊魂等，2015).本文数值模拟分析模型的岩体弹性模量参考花岗岩取值，为4.0×10¹⁰~10.6×10¹⁰Pa(沈明荣和陈建峰, 2006)，对整个模型按均匀梯度对弹性模量和密度进行赋值，其他参数为定值，物理力学参数如表 1所示.

表 1 Table 1

表 1 模型岩体物理力学参数 Table 1 Physical parameters of rock for the model 力学参数 弹性模量/GPa 抗拉强度/MPa 内聚力/MPa 摩擦角/(°) 泊松比 密度/(kg·m-3) 重力加速度/(m·s-2) 地表 底部 地表 底部 数值 40 106 12 16 35 0.286 2575 2809 9.8

表 1 模型岩体物理力学参数 Table 1 Physical parameters of rock for the model

图 5描述了龙门山断裂带3条主要断裂的形成过程.由该图可见，区域构造运动6万年之前，龙门山地表出现塑性破坏区，地壳下方岩体保持完好，处于弹性状态；区域构造运动6.39万年时，在马尔康下方深度约32.7 km、水平位置为29 km处开始有塑性区产生，此处为F1断层发育初始位置；到6.72万年时，其右上方的塑性区发育较为明显，F1断层初露端倪；区域构造运动7万年时，F1断层发育至地表，在莫霍面附近呈现不对称共轭塑性区；到10.07万年时，地表下方22.0 km、水平57 km处出现塑性区，此处为F2断层发育初始位置；到10.25万年时，F2断层由中部向两侧发育；到11万年时，F2断层发育至地表，F2的底部同时发育不对称共轭断层；到14.51万年时，在深度为24.2 km、水平69.9 km位置处出现塑性区，此处为F3断层发育初始位置；到15万年时，塑性区继续发育；到16万年时，F3断层基本发育至地表；到17万年时，其底部同时发育了不对称共轭断层.

图 5

Fig. 5

图 5 龙门山断裂形成过程模拟图 Fig. 5 Simulation of formation process for the Longmenshan fault zone

因此可认为，F1、F2、F3断层是从中部向上下两侧开始发育，上侧逐渐发展至地表，而且F1、F2和F3断层依次形成断层.这些塑性区的发育方向与水平面的夹角大致呈40°左右，并出现了各自的共轭断层，与实际的龙门山断裂大致吻合(胡幸平等，2012).

由图 5还可以看出，经过长期的区域构造运动作用，在区块岩体内形成了纵横交错的塑性区，它们随着区域构造运动而不断发生变化，新的塑性区和断层不断产生.这种断层不断产生的过程也许会引发不同的地震，当区域构造运动活跃时，产生的地震会增多，反之亦然；同时，随着区域构造运动，地表的破坏深度也具有很大差别，破坏深度大的区域或许形成谷地，而破坏深度小的地域成为山峰.

提取测点1~测点4的最大最小主应力数值，处理为随时间变化的曲线图. 图 6a为测点1的最大最小主应力曲线，区块运动4万年之前, 最大最小主应力基本随时间的增大而增大，即随着区域构造运动的持续，主应力随之增大，符合岩体应力随着应变量的增大而增大的规律.当区域板块运动到一定程度时，测点的主应力数值达到相对稳定的状态，但达到稳定状态的时间不同，如图 6b—图 6d所示，测点1~测点4的主应力达到稳定状态的时间大致为4万年、10万年、12万年、17万年.地壳岩体水平挤压会导致水平应力的增大，而地壳岩体变形、塑性破坏等又会导致地应力的释放，因此，地壳岩体内部因区域构造运动而增加的主应力与释放的主应力可能形成相对平衡的状态，该平衡机制类似于自组织临界性，即一个具有持续能量供给且由许多基本单元组成的系统，当单元之间具有相互作用时，会自发地演化到某个临界状态(Bak et al., 1987；易桂喜等，2013).

图 6

Fig. 6

图 6 不同测点最大(σ1)、最小(σ3)主应力随时间变化曲线图 (a)深度为1 km; (b)深度为11 km; (c)深度为21 km; (d)深度为31 km. Fig. 6 Curves of maximum (σ1) and minimum (σ3) principal stresses varying with time at different sites (a) Depth=1 km; (b) Depth=11 km; (c) Depth=21 km; (d) Depth=31 km.

为研究断层形成前、F1断层贯通、F2断层贯通和F3断层贯通等几个关键时间点的主应力变化规律，截取6万年、7万年、11万年和16万年时的4个测点的最大最小主应力数值，计算其比值如表 2所示.处于不同深度的各测点最大最小主应力虽然数值变化较大，但比值变化较小，测点1在6万年~16万年的最大最小主应力比值介于4.14~4.39，平均为4.3.深部测点也大致呈同样的规律，测点2、测点3和测点4的最大最小主应力比值变化幅度同样较小，平均比值分别为3.41、3.43和3.39.

表 2 Table 2

表 2 测点主应力及其比值汇总表 Table 2 Summary table on principal stress and its ratio of measuring point 时间/万年 测点1 测点2 测点3 测点4 σ3/MPa σ1/MPa σ1/σ3 σ3/MPa σ1/MPa σ1/σ3 σ3/MPa σ1/MPa σ1/σ3 σ3/MPa σ1/MPa σ1/σ3 6 86.35 378.84 4.39 252.83 859.42 3.40 356.22 1189.91 3.34 469.37 1554.43 3.31 7 86.06 377.37 4.38 296.81 1003.76 3.38 405.36 1357.02 3.35 526.87 1756.16 3.33 11 87.59 376.57 4.30 385.65 1347.32 3.49 529.03 1801.95 3.41 679.87 2320.11 3.41 16 89.47 370.19 4.14 391.00 1320.74 3.38 556.65 2021.33 3.63 770.84 2710.34 3.52

表 2 测点主应力及其比值汇总表 Table 2 Summary table on principal stress and its ratio of measuring point

根据蝶形塑性区理论(Zhao et al., 2013；马念杰等，2015)，塑性区的发生发展主要是由偏应力所决定，偏应力能控制岩体的破坏，当最大最小主应力比值η达到某一临界条件(K₁)后，圆孔(或地质变异体，或最先出现的塑性体)周边的塑性区呈“X”形共轭分布，塑性区半径R_max增长趋势会发生急剧扩展现象.从本模型计算的结果来看，模型的最大最小主应力比值η在3.31~4.39之间，表明持续的水平运动，可能会致使最大最小主应力的持续变化，当其比值达到临界状态时，塑性区的发展将发生无限扩展，这也是断层从深部发育至地表的原因之所在.

(1) 板块水平运动是地壳主应力变化的主要原因，岩体内应力增大至岩体强度极限时，必将产生塑形破坏.在一定应力环境中岩体的塑形区发育具备2个特性：一是方向特性，即断层的发育延伸具有一定的方向性，本计算模型的塑性区发育倾角约为40°左右；二是共轭特性，在大塑性区附近会形成近似垂直相交的不对称共轭塑性区.

(2) F1、F2、F3断层的发育是从中部开始，逐渐发展至地表，F1断层首先出现，F2和F3断层依次出现.模拟形成的断层几何特征与实际的龙门山3条断层大致吻合.

(3) 区块内主应力变化分为2个阶段，一是在区块挤压初期，主应力随着区块内应变量的增大而增大；二是主应力相对稳定阶段，当区块挤压持续进行时，区块岩体内的主应力保持相对稳定的状态，表明区块岩体内部存在应力积聚-应力释放的平衡机制.

(4) 不同深度岩层中的最大最小主应力比值均存在临界状态，浅部较大而深部略小，较高的偏应力环境可能是断层形成的重要原因.

感谢本文的评审专家提出的中肯建议，使得文章内容更加完善.