雁列构造是一种野外常见的非连续地质构造.拉张型和挤压型是两类常见且重要的雁列构造.地震界、地质界等领域的众多科技人员开展了大量的与雁列构造有关的研究工作(Segall and Pollard， 1980；Sibson，1985；Ma et al., 1986，2010；Aydin and Schultz， 1990；Swanson，1990；Thomas and Pollard， 1993；Zachariasen and Sieh， 1995；Talwani，1999； de Joussineau et al., 2003；Barker，2005；Smith and Yamauchi， 2005；Soliva et al., 2010).一些关于断层相互作用、地震机理、前兆及过程的研究也常以雁列构造为例进行分析和讨论(Sibson，1985；刘力强等，1986；Ma et al., 1986，2010；马文涛等，1995；马胜利等， 1995a，1995b，2004，2008；马瑾等， 1999，2002，2007，2008；蒋海昆等，2002；刘培洵等，2007；任雅琼等，2013).目前，关于雁列构造的研究主要集中在理论分析、实验研究、野外观测及统计等方面.相比之下，数值模拟方面的研究成果少见(Aydin and Schultz， 1990；Olson and Pollard， 1991；Du and Aydin， 1995；Harris and Day， 1999； de Joussineau et al., 2003；Ohlmacher and Berendsen，2005；Palusznya and Matth ib，2009).随着计算机软、硬件水平的不断发展，数值方法在上述问题研究中可望能发挥越来越重要的作用.采用合适的数值模型，开展实验室尺度模型的研究工作，是实验研究的有力补充，并能积极促进或推动有关的研究向前发展，极大地增进人们对于断层相互作用及失稳前兆的认识和理解.在数值模拟之中，一些费时费力甚至难于在当前技术条件下实现的实验方案，可能方便地实现；一些困扰实验研究的难题可以方便地被克服或避免，例如，试验机压头之间的干涉，多种观测元件之间的干扰，高围压、水压的施加，同批次标本之间或多或少存在差异，声发射监测中的“死时间”问题，观测区域受限或观测变量受限等问题.

过去针对雁列构造的数值模拟研究，多针对应力场、应变场、位移场及动态断裂过程等进行分析，考察的力学量有限.在一些数值模型之中，不允许单元破坏，因而难于研究雁列区的贯通过程及引起的多种力学量的变化规律；少见非均质性被引入.事实上，材料或几何非均质性的重要性怎么强调也不过分.在连续介质模型之中引入某种程度上的非均质性，即可模拟地震的时空复杂形为(Miyatake，1992；Beroza and Mikumo， 1996；Bizzarri et al., 2001；Hillers et al., 2006).

近年来，作者开展了典型断层构造破裂失稳前兆及粘滑过程的数值模拟研究(Wang et al., 2010，2012，2013a，2013b；王学滨等， 2013，2014).考虑到上述问题研究的难度，将上述问题分别对待，采用两种模型分别进行分析和讨论.

在第1种模型中，分别考虑了断层和岩石的非均质性，针对雁列构造及Z字型构造开展研究.这种模型适于模拟一个粘滑周期内断层系统的复杂力学行为.研究发现，挤压雁列区贯通过程中b₀值(事件的频次-能量释放关系的斜率的绝对值)表现为下降(Wang et al., 2010)，而对于拉张雁列构造则观察不到b₀值的明显变化(Wang et al., 2010，2012).在挤压雁列区贯通过程中，与剪切应变降有关的量的表现比较突出，这是由于1个单元的破坏，能引起周围众多单元剪切应变的降低，即弹性卸载(Wang et al., 2013b).

在第2种模型中，未考虑非均质性，断层系统由断层和之外的破坏区(针对雁列构造而言)以及弹性岩石构成.对于断层单元，在应变软化之后考虑了应变硬化，以模拟随着断层滑动距离的增加，断层面的强度愈合现象.采用这种模型，模拟了单一、交叉、拐折及雁列构造的粘滑过程.针对拐折构造的计算发现，大粘滑事件与断层不同段落的同步失稳有关，而小粘滑事件与不同区段的非同步失稳有关(Wang et al., 2013a).

为了探索事件的迁移及前兆规律，本文以拉张型雁列构造为例进行数值分析和讨论.数值计算在FLAC-3D软件中实现，采用了提出的非均质应变软化模型(Wang et al., 2010，2012，2013b).为了研究两条断层的相互影响和作用，在两条断层的内、外侧岩石上布置了若干被监测单元，监测它们的剪切应变的演变规律.在标本全场内对大于某一临界值的剪切应变降进行了捕捉，并对与此有关的3种量进行了统计.文中还阐明了雁列区内的破坏区不启动于断层端部现象的原因.1 计算模型及主要参数

图 1中岩石标本包含两条具有一定重叠量的平行断层，在重叠区域构成雁列构造.断层的重叠量和间距相等，均为0.0184 m.在纸面内，标本的形状为正方形，边长为0.3 m，被剖分成9万个尺寸相同的单元.两条断层与水平方向的夹角均为45°.

在平面应变条件下进行计算，仅允许单元发生小变形.计算过程包括两个阶段：

首先，在静水压力(σ₁=σ₃=-2 MPa)条件下，将标本加载至静力平衡状态.这一阶段消耗2万个时步.若标本中节点的最大不平衡力已经足够小则意味着标本已经达到静力平衡状态.一个时步是指一次计算循环，从运动方程经由几何方程，到本构方程，再到运动方程.循环一次，则时步数目增加1.

其次，在标本左、右两端进行压缩的位移控制加载，加载速度v=1×10^-9 m/时步，而标本上、下两端的压应力与第1阶段相同.通过计算标本左、右两端上各节点不平衡力的合力与这两端的几何尺寸的商得到这两端的压应力.计算发现，该压应力的值大于标本上、下两端的压应力的值(2 MPa)，所以应为σ₁.由于两断层端部连线方向与σ₁垂直，所以图 1中的雁列构造为拉张型.

标本中包含两种单元：岩石单元和断层单元.这两种单元的本构模型相同，但有些参数的取值不同.在弹性阶段，两种单元均服从各向同性的线弹性模型；在破坏阶段，两种单元均服从带拉伸截断的应变软化莫尔-库仑模型.

图 1

Fig. 1

图 1 拉张雁列构造的计算模型(a)静水压力加载阶段；(b)水平方向位移控制加载阶段.Fig. 1 A computational model for the extensional echelon fault(a)The initial loading stage；(b)Horizontal displacement-controlled loading stage.

为了延缓单元的破坏过程，增加发现前兆的可能性，我们对传统的均质模型进行了改进(Wang et al., 2010，2012，2013b)，假定一些材料参数在空间上呈不均匀分布，即弹性模量、初始粘聚力和抗拉强度在标本内部服从Weibull分布，非均质性参数取值适中，即m=9，且它们彼此之间不相关.其余参数被认为是均匀的，例如泊松比、初始内摩擦角等.

对于岩石单元和断层单元，弹性模量的均值分别取为55 GPa和5.1 GPa；泊松比分别取为0.25和0.2；初始粘聚力的均值分别取为37.5 MPa和5 MPa；抗拉强度的均值分别取为24 MPa和1.2 MPa；初始内摩擦角分别取为50°和10°，扩容角均取为0°.

为了模拟应变软化，除了需要规定强度参数的初始值，还需规定软化阶段结束时的残余强度参数和所对应的塑性剪切或拉伸应变的值.我们假定，随着单元塑性变形的增加，单元的强度参数倾向于相同.也就是说，具有不同初始强度参数的单元的残余强度参数及所对应的塑性应变相同.这也意味着，强度参数越高的单元的峰后脆性越强.

上述参数的具体取值及依据参见Wang等(2010，2012，2013b). 2 计算结果分析

2.1 破坏区分布

图 2给出了雁列区贯通过程及之前不同时步数时标本中破坏单元的分布规律.黑色单元代表已发生剪切或拉伸破坏的单元.

由图 2可以发现，随着时步数的增加(或σ₁方向线应变的增加)，破坏单元由远离雁列区的断层上逐渐向雁列区迁移，直至雁列区贯通.

在4.8万个时步数之前(图 2a~c)，雁列区内尚未出现破坏单元.当时步数达到4.9万个时(图 2d)，在右断层某一位置(过左断层端部且垂直于两条断层的直线与右断层的交点处)扩展出了破坏区.当时步数达到5.0万个时(图 2e)，扩展出的破坏区已接近左断层端部.当时步数达到5.1万个时(图 2f)，雁列区已被上述破坏区贯通，与此同时，在左断层某一位置(过右断层端部且垂直于两条断层的直线与左断层的交点处)也扩展出了破坏区.随后，这条破坏区也将贯通雁列区.

图 2

Fig. 2

图 2 不同时步数时破坏单元的分布(a~f)时步数分别为4.1万、4.6万、4.8万、4.9万、5.0万及5.1万.Fig. 2 Distributions of failed elements for different timestep numbers(a~f)The timestep number is 4.1×104，4.6×104，4.8×104，4.9×104，5.0×104 and 5.1×104，respectively..

上述破坏区不启动于断层端部，而启动于断层某一位置的现象似乎有些反常.显然，对于标本中仅有一条断层的情形而言，这种现象不应出现.在实验室及野外，也能观察到上述反常现象(图 3).在图 3a中，拉张雁列区内部的裂纹不启动于两条断层的任一端；在图 3b中，大量声发射事件集中于下方断层的某一位置，预示着这一位置的开裂；在图 3c中，两条断层之间的新断裂启动于下方断层的某一位置.据此推测，上述反常现象的出现应与两条断层之间的相互影响和作用有关.

图 3

Fig. 3

图 3 实验室及野外发现的反常裂纹启裂位置(a)陈顺云等(2005)；(b)马胜利等(2008)；(c)Sibson(1985).Fig. 3 Abnormal crack initiation positions observed in experiments and in situ(a)Chen et al.(2005)；(b)Ma et al.(2008)；(c)Sibson(1985).

为了揭示上述反常现象的原因，我们给出了一条反常裂纹启动前后雁列区附近的速度场(图 4).标本在水平方向压缩条件下，右断层下盘的岩石主动地向左运动，而上盘岩石被动地斜向上(或垂直于右断层方向)运动.随着时步数的增加，向上运动的岩石尺寸越来越大，其前锋越来越向雁列区逼近.由于拉张雁列区之内的应力是以拉应力为主，当其前锋抵达拉张雁列区的边缘时，过高的拉应力将使前锋处的岩石单元发生拉破坏，进而诱发出垂直于右断层方向的拉裂纹.

图 4

Fig. 4

图 4 反常裂纹启裂前后雁列区附近的速度场(a~d)时步数分别为4.8万、4.85万、4.9万和4.95万；圆圈位置表示斜向上运动的岩石的前锋. Fig. 4 Velocity fields prior to the abnormal crack initiation and beyond(a~f)The timestep number is 4.8×104，4.85×104，4.9×104 and 4.95×104，respectively.

在FLAC-3D中，剪切应变增量定义为

为了获得标本中特定位置剪切应变的演变规律，在两条断层的内、外侧均布置了一些被监测单元(图 5).在某条断层的任一侧，被监测单元的连线平行于该断层，且被监测单元等间隔布置.被监测单元均布置在岩石单元上.左断层内侧、外侧、右断层内侧、外侧的被监测单元编号分别为3~32、33~62、63~92及93~122.值得指出的是，有一些被监测单元布置在断层端部的前方，例如，被监测单元3~8等.

图 5

Fig. 5

图 5 部分被监测单元布置Fig. 5 Arrangement of a part of monitored elements

图 6给出了左断层内侧(包括断层端部前方)的多个被监测单元的剪切应变随时步数的演变规律，同时，也给出了标本左、右端面的σ₁的值随时步数的演变规律(标记为1的曲线).红色箭头指明了雁列区贯通的时刻.由此可以发现：

(1)在雁列区贯通之前及之后，σ₁的值总体上呈现上升的态势，但前、后上升速度不同.在大约4.5万个时步数之前，σ₁的值上升较快，但此后，σ₁的值上升放慢.这与马胜利等(2008)的实验结果(差应力随时间增加的速度明显减缓)基本一致，雁列区的贯通难于在σ₁的演变规律中找到明显的标志，但却引起了某些位置被监测单元剪切应变的快速下降，例如被监测单元9、12及15，它们均位于雁列区内部，这与这些单元附近发生破坏有关.马胜利等(2008)通过实验也发现，靠近雁列区的应变释放极为明显.反之，对于被监测单元3及7，由于附近没有破坏单元，所以在雁列区贯通过程中，没有观察到剪切应变快速下降现象.

图 6

Fig. 6

图 6 最小主应力(σ1)及左断层内侧被监测单元的剪切应变随时步数的演变规律Fig. 6 Evolution of the minor principal stress(σ1) and shear strains of monitored elements at the inner side of the left fault with timestep number

(2)在雁列区贯通之前及之后，被监测单元的剪切应变的演变规律有重大的差别.在雁列区贯通之前，可以观察到一些被监测单元的剪切应变的缓慢上升和快速下降现象，这样的现象重复了多次(图 6中蓝色虚箭头)，就像粘滑一样；离雁列区越远，剪切应变下降现象出现得越早；在雁列区贯通之后，仅观察到了剪切应变的小幅度震荡现象.马胜利等(2008)通过实验发现，在雁列区贯通之前，断层内侧各点应变在随时间缓慢增加的背景上伴有准周期性的台阶式下降.这与本文的上述结果基本类似.

图 7a~b给出了左、右断层内侧被监测单元的剪切应变的演变规律.由此可以发现下列现象：

(1)在虚线1及1′的位置，剪切应变降仅发生在左断层内侧，例如，第23、25号被监测单元.

(2)在虚线2及2′的位置，剪切应变降同时发生在两条断层内侧，例如，第23、89号被监测单元.

(3)在虚线3及3′的位置，剪切应变降仅发生在右断层内侧，例如，第79号被监测单元.

(4)在虚线4及4′的位置，剪切应变降先发生在左断层内侧，后发生在右断层内侧，例如，第25、79号被监测单元.

(5)在虚线5及5′的位置，剪切应变先发生降低(在左断层内侧)，然后发生突增(在右断层内侧)，例如，第27号和第87号被监测单元.

至少，部分现象(左、右断层内侧剪切应变降同时发生等)应该根源于断层之间的相互影响和作用.两条断层附近的剪切应变的演变规律非常复杂，同步、不同步共存.马瑾等(2002)从剪应力的角度分析了断层相互作用引起的促震和减震效应.从本文的结果看，断层之间的相互作用对于有的部位产生了减震的作用(剪切应变的突然降低)，而对于有的部位产生了促震的作用(剪切应变的突然增加).

利用图 7a，c可以分析左断层内、外侧被监测单元的剪切应变的演变规律.由此可以发现，通常，这些被监测单元的剪切应变的演变都能表现出较好的同步性.例如，在虚线1及1″的位置，当左断层内侧的第23、25号被监测单元的剪切应变下降时，在左断层外侧的第55、57号被监测单元的剪切应变也下降，第25、55号被监测单元(编号相差30)的连线刚好垂直于左断层；在虚线2及2″的位置，第23、53号被监测单元的剪切应变同时表现为下降；在虚线4及4″的位置，第29、59号被监测单元的剪切应变同时表现为下降，等等.上述计算结果表明，通常，可以仅对断层一侧的被监测单元的剪切应变进行分析即可了解断层上事件的发生情况.上述结果与马胜利等(2008)的实验结果有所差异，他们发现，在雁列区贯通之前，断层内侧各点应变的释放引起外侧应变的台阶式上升，当雁列区贯通之后，断层内、外侧应变变化才趋于一致.

图 7

Fig. 7

图 7 被监测单元剪切应变随着时步数的演变 (a)左断层内侧；(b)右断层内侧；(c)左断层外侧；虚线1、1′及1″对应相同的时步数，依此类推.Fig. 7 Evolution of shear strains of monitored elements with timestep numbers(a)At the inner sides of the left fault；(b)At the inner sides of the right fault；(c)At the outer sides of the left fault；the dotted lines1，1′ and 1″correspond to the same timestep number， and so on.

我们重点关注那些较大的剪切应变降发生的地点及迁移规律.在图 8中，给出了超过5×10^-6的剪切应变降在不同时段内的空间分布规律.如果将上述阈值取得过小，则可能标本各处都有剪切应变降；如果取得过高，可能捕捉不到剪切应变降.根据经验，上述阈值比较合适.每个时段持续1千个时步.第1个时段的结果为4万至4.1万个时步数之间的结果，依此类推.图 9给出了对剪切应变降的一些统计结果，分别是剪切应变降的的累计(在任一时段之内，将所有单元的剪切应变降量求和)、最大值及发生剪切应变降的单元数目.在图 8中，圆圈的半径定量代表剪切应变降的大小，剪切应变降的比例尺均相同，图 8各子图中的最大剪切应变降易于从图 9中看出.

图 8

Fig. 8

图 8 不同时段时剪切应变降的空间分布 (a)时段1;(b)时段4;(c)时段5;(d)时段7;(e)时段8;(f)时段9;(g)时段10;(h)时段11;(i)时段13;(j)时段14;(k)时段15;(l)时段16. Fig. 8 Distributions of shear strain decreases for different timestep intervals

图 9

Fig. 9

图 9 与剪切应变降有关的3种量的演变规律 Fig. 9 Evolution of three kinds of quantities related to shear strain decreases

由图 8可以发现，随着σ₁的值的增加，剪切应变降由远及近迁移.也就是说，由远离雁列区的断层附近，向雁列区附近迁移；之后，剪切应变降云集在雁列区附近；最终，仅两条断层附近个别位置发生剪切应变降.

在4.4万个时步数之前(第1~4个时段，图 8a~b)，在雁列区附近，未观察到剪切应变降现象.这意味着雁列区及附近的断层尚未破坏.在远离雁列区的断层附近，剪切应变降密集分布，反映了在断层上发生了大量的破坏事件.剪切应变降的累计及发生剪切应变降(二者的演变规律一直类似)的单元数目在第4个时段内经历了最高值(图 9中左数第1条虚线).

当达到4.5万个时步数时(第5个时段，图 8c)，左断层附近的剪切应变降已逼近雁列区，在右断层端部附近，已观察到剪切应变降.

当时步数在4.6~4.7万之间时(第7个时段，图 8d)，尽管σ₁的值持续增加，但在雁列区附近，未能观察到剪切应变降，雁列区形成了一个空区.此时，在左断层附近出现一个全局(第1~16个时段)最大的剪切应变降，其值约1.7×10^-2(图 9).随后，当时步数在4.8~5.0万之间时(第9~10个时段，图 8f~g)，上述空区消失.此时，雁列区内已出现破坏单元.

当时步数在5.0~5.1万之间时(第11个时段数，雁列区发生贯通，图 8h)，剪切应变降云集在雁列区及附近的断层上，统计的与剪切应变降有关的3种量均不高(图 9中左数第3条虚线).

此后，剪切应变降累计及发生剪切应变降的单元数不断变少，在特定时段内，仅发生在两条断层的某些特定位置上，这些位置可能具有高的强度，因而有时剪切应变降的最大值并不小.

综合分析图 8可以发现，在4.1万个时步数之前，在特定时段之内，剪切应变降仅发生在某些特定位置上，它们彼此孤立.这应该是断层上少量的低强度单元破坏的反映.随后，随着σ₁的值的增加，剪切应变降明显增加，发生剪切应变降的区域连成片，密集地分布在远离雁列区的断层附近.这应该是断层上大量的中等强度单元破坏的反映.当雁列区贯通之后，两条断层上某些位置仍可以观察到剪切应变降.这应该是断层上少量的高强度单元破坏的反映.

总之，随着σ₁的值的增加，标本中的剪切应变降的表现先从平静变得活跃，雁列区曾一度成为空区，随后空区消失.当拉张雁列区贯通之后，剪切应变降的表现再一次回归平静.马胜利等(2008)通过实验发现，对于拉张雁列构造，雁列区破坏发生在前，而断层滑动发生在后；雁列区的破坏对断层滑动具有指示作用.在雁列区贯通之前的剪切应变降的活跃表现及雁列区之内空区的出现均具有一定的断层失稳前兆意义.由平静到活跃反映了标本应力或应变的增加而引起的发生在断层上的事件越来越多，空区的出现反映了标本应力或应变的增加过程中，拉张雁列区对断层错动一定程度的阻碍作用，尽管这种阻碍作用和挤压雁列的阻碍作用相比要小得多. 3 结 论

3.1 在拉张雁列区贯通过程中，观察到了启动于断层某一位置的破坏区，这根源于标本应力增加过程中两条断层的相互影响和作用.

3.2 在拉张雁列区贯通之前，观察到了一些被监测单元剪切应变的多次缓慢增加及之后的快速下降现象，这反映了标本应力增加过程中，断层上单元的不断破坏造成的卸载.

3.3 在拉张雁列区贯通之前，观察到了剪切应变降由远及近的迁移、雁列区成为空区及空区消失等现象，这些现象可视为断层失稳的前兆.

3.4 随着标本应力的增加，两条断层附近的剪切应变降依次表现为平静、活跃及平静.前后两个平静期与少量的低或高强度单元破坏有关，活跃期与大量的中等强度单元破坏有关.

3.5 在同一时刻，一条断层两侧相对位置的被监测单元的剪切应变演变规律比较类似，但不同断层附近的剪切应变演变规律非常复杂：有时突增，有时陡降，同步、非同步现象共存.这些复杂的表现根源于断层的非均质性及复杂的相互影响和作用，值得深入研究.