共轭断层系统的非线性有限元模拟与震群模型讨论

郭婷婷 徐锡伟 邢会林 于贵华

1)山东省工程地震研究中心，济南 250021

2)山东省地震局，济南 250014

3)中国地震局地质研究所，活动构造与火山重点实验室，北京 100029

4)Earth Systems Science Computational Centre(ESSCC)，University of Queensland，St Lucia QLD 4072，Australia

0 引言

野外地震现场观测和岩石力学实验表明，地震的孕育和发生与活动断层的滑动习性密切相关，不同类型活动断层黏滑错动引起断层两盘块体弹性应变的快速释放是地震成核和发生的主要机制(Scholz，1990)。断层面上凹凸体或障碍体的几何形态和规模大小等几何结构和断层本身的摩擦非均匀性等均会对断层滑动失稳过程产生影响(马瑾，1999)。另外，共轭断层是最常见的断层分布样式之一，中国大陆许多强震的震源断裂表现为共轭断层，许多地震区存在相互交叉的共轭发震断层(张四昌，1991)，特别是华北平原、川滇菱形块体等地区常常发生与共轭活动断层错动相关的双震或震群型地震;例如华北平原1966年邢台地震群(M6.8，6.7和7.2)、1976年松潘－平武地震群(M7.2，6.7和7.2)、1988年澜沧－耿马地震(M7.6和7.2)等(Jones et al.，1984;俞维贤等，1994;徐锡伟等，2000)。川滇菱形块体与华南克拉通之间鲜水河－小江断裂系东侧，2014年8月3日发生的鲁甸M6.5地震在地表沿NW向包谷垴－小河断裂形成的2.2km长的地表破裂带表明，发震断层为NW向左旋走滑断层(徐锡伟等，2014);但重新定位后的余震呈NNW和NEE向“L”状分布，显示出共轭断层之间存在着破裂触发作用，类似于1987年美国加利福尼亚州Superstition Hills地震序列在NE向和NW向2组共轭断层上同时出现地震地表破裂带的现象，说明发育程度不同的2组交叉断层或共轭断层在孕震或发震过程中存在着相互作用或触发作用(Hudnut et al.，1989;王未来等，2014)。曾有学者用NE向与NW向2组共轭断层及其围限块体的转动与交替活动模型来解释首都圈及其邻区的现今构造变形与地震活动性:地震平静期的构造变形主要受易于错动的NW向断层控制，而大地震的发生则是由NE向断层的快速错动引起，或者出现断层黏滑(强地震)与断层蠕滑(弱地震)交替错动的现象，给出了交叉或共轭断层交替滑动现象的部分理论解释(马瑾等，2000;马胜利等，2003;马瑾，2009)。但在这种共轭断层系统中，有关NW向断层在区域孕震构造中的作用及其与NE向断层关系的问题，尚需作进一步的力学解释和论证(Scholz，1990;马瑾等，2000;徐锡伟等，2002;马胜利等，2003;马瑾，2009)。因此，研究共轭断层这种典型地震构造的相互作用及块体运动具有十分重要的科学意义。

有限单元方法是数值模拟的一种常用方法，广泛应用于地球动力学研究与断层系统的模拟中，并取得大量研究成果(Jungels et al.，1973;汪素云等，1980;王仁等，1980，1982;Melosh et al.，1989;陈连旺等，1999;Cai et al.，2000;朱守彪等，2004，2008;Xing et al.，2006;周仕勇，2008;姚琪等，2012a，b)。因此，本文针对共轭断层系统这一典型的地震构造模型，采用断层摩擦接触的非线性数值模拟方法，模拟共轭断层上交替发生的失稳事件，以及NE向与NW向2组共轭断层系在孕育与发生破坏性地震过程中的作用，用以解释中国大陆地区典型共轭断层系从应变能积累的闭锁状态到突然释放的解锁过程。这一研究对进一步深刻理解华北平原、川滇地区震群型地震的空间分布规律、不同构造部位应力－应变响应机理、间震期和同震构造变形的协调性问题，乃至强地震发生前可能的前兆异常的识别等均具有十分重要的理论和实际意义。

1 摩擦接触非线性有限元模拟方法

采用摩擦接触非线性方法处理断层问题，可有效地解决先存断层间的相互作用问题，并利用库仑摩擦准则解释从地震孕育到发生的复杂力学过程:断层从黏着状态到滑移状态，进而产生地震的行为。本文数值计算所用程序Esys-crustal是基于摩擦接触的有限元平衡方程，其中，采用与速率相关的摩尔库仑摩擦准则来描述断层接触面的黏着与滑移状态，通过基于R最小策略的显示算法，进行积分求解(Xing et al.，2002)。

1.1 有限单元平衡方程

采用更新的拉格朗日表达式描述非线性摩擦接触问题，根据增量虚功率原理(Xing et al.，2002)，可以将平衡方程表述如下:

式(1)中，V和S分别表示t时刻整个变形体的体积与表面积;SF是外载力作用在S表面的边界;SC是接触表面;δ˙u是在速度边界上满足的虚速度场是Cauchy应力的Jaumann率;L是速度梯度张量，且是张量L的对称部分;˙f是接触表面SC上接触力的速率;是接触对之间的相对滑动速度。

1.2 摩擦接触的本构关系

大多数构造地震起因于先存断层或板块边界的突然摩擦滑动，Brace等(1966)认为地震是岩石的黏滑失稳过程(Dieterich，1979;Ruina，1983)。在岩石的摩擦滑动中，断层接触面上具有规则的黏着(sticking)和滑动(slipping)现象，故地震是沿着先存断层面上的滑动不稳定现象。

1.2.1 法向接触力的计算

将先存的2组共轭断层面理解为变形体之间的接触，其中1组为主接触面(master)，另1组为次接触面(slaver)。运用罚函数方法处理法向约束，以确保断层面接触时法向接触力的传递以及2接触面间没有穿透;对于次接触面，其法向接触力为

式(2)中，En是罚因子，用以防止接触面发生穿透，gn是法向穿透距离，且gn=n·(xs－xc)。

1.2.2 接触摩擦力的计算

采用摩尔库仑摩擦准则来描述断层接触面上的黏着(sticking)与滑动(slipping)状态，相当于理想弹塑性材料变形的屈服准则，将黏着状态类比为弹性，滑动状态类比为塑性。为了避免状态变量计算导致的计算困难，采用速率相关的Coulomb摩擦准则，使用将摩擦中黏着与滑移两部分分开的办法处理，而不考虑状态量变化的效应，若同时忽略温度影响，则断层接触面上的摩擦力描述为

式(3)、(4)中，μ是摩擦系数，它取决于法向接触力fn，等效切向速度状态变量φ，即，Dieterich，1979;Ruina，1983;Scholz，1998);Et是沿接触面切向方向的一个常量;是临界摩擦力，，其中的初始值;的滑动部分(不可逆)，则是的闭锁部分(可逆)。

1.3 模拟方法的优点

本文采用摩擦接触非线性有限元程序Esys－crustal，对共轭断层系统这类地震构造模式进行数值模拟、对比和分析，进而解释这类构造模式在实际地质中的构造意义。其主要优点如下:

(1)采用有限元数值模拟方法，通过数值模拟岩石断层的破裂过程与机理，探究相似发震构造的实际地壳结构，具有很好的可行性，能更真实地反映物性参数与断层接触关系。

(2)以往的数值模拟常将断层处理为一软弱带，通过物性参数的选取，降低其强度。但这样处理断层无法模拟地壳介质中存在的不连续间断面，无法模拟断层面上存在的滑动摩擦机制。采用接触摩擦分析单元处理断层，能较好地模拟断层面上存在的摩擦滑动机制。

(3)采用显式的时间积分算法(Xing et al.，2002，2006)，通过时间步长增量的比例控制系数(Rmin＜1)来控制荷载增量的大小。即选取大小合适的时间步长，使得每个增量步内的单元力学性状与界面上的接触状态平稳变化，从而保证整个计算过程稳定可靠。

2 地震构造模型与加载方式

2.1 地震构造模型

根据华北平原和川滇地区发育交叉或共轭断层地震区已有震源机制解、区域构造应力场和发震构造或孕震构造研究成果可知，发生震群型地震的地区区域主压应力近于水平，发震断层基本上为走滑型断层，一组为左旋走滑断层，与其垂直的另一组为右旋走滑断层(邓起东等，1976;李钦祖等，1982;Jones et al.，1984;许忠淮等，1989;张诚等，1990;张四昌，1991;俞维贤等，1994;徐锡伟等，2000，2002;马文涛等，2004;崔效锋等，2005;周仕勇，2008)，其中一组为主断层，另一组为次断层。例如，华北平原北部地区 N30°E右旋走滑断层为主断层，而N60°W左旋走滑断层为次断层。因此，可将发震构造模型简化为单条N30°E断层模型以及N30°E和 N60°W向共轭断层系统模型(简称交叉断层模型)，且均为走滑型;另外，最大主压应力轴优势方位取近EW向，最小主压应力优势方位为SN向，以保证NE向断层以右旋走滑为主，NW向断层以左旋走滑为主。通过上述2种不同的地震构造模型的数值模拟与比较，探究NW向次断层在共轭断层系统强震孕育中的作用。

2.2 断层面与介质性质

在断层相互作用的构造力学分析中，关键与难点之一是对断层的正确处理。本文将断层面两侧块体考虑为弹性体，断层接触面两侧为摩擦非线性接触关系。断层接触面上的接触摩擦分析可真实反映出断层存在所引起的地壳不连续面，断面上建立主、次接触面，其滑动遵从库仑摩擦准则，即一旦断层面上的剪应力超过极限摩擦强度，断层由闭锁状态进入黏滑状态，即发生地震。

2.3 计算模型与加载方式

2.3.1 计算模型简化

鉴于地震构造模型和选定的构造应力场方位、断层面与介质性质等基本特征，依据相似性原理，计算模型选取100mm×100mm×15mm的1个3维块体，其中断层取为竖向直立断层。断层为接触面，分别定义为主面与次面。在有限元计算中，利用摩擦接触单元模拟摩擦接触，其接触面上的摩擦依速度而定，服从下列公式:

式(5)中，μ为摩擦系数，V与Vref分别是滑动速度和参考速度，本文取Vref=10－4m/s。

计算模型中的断层模型简化为以下2种模式:

(1)N30°E发震断层的单断层模式(图1):本模型中只考虑单条N30°E发震断层，计算在构造应力场作用下时间步长的相对速度变化，并分析速度突变造成的一次地震事件的发生以及其与计算时间步长的关系。

(2)N60°W与N30°E两组共轭断层系统的交叉断层模型(图2):本模型中不仅考虑到N30°E向主断层，同时还考虑了N60°W向次断层。将强震孕育、发生的2组共轭断层看作1个断层系统，计算该断层系统的应力－应变关系，寻找断层面上的速度突变，分析计算时间步长与速度突变大小的关系，以及孕震构造在断层系统强震孕育中的作用。

图1 N30°E单断层模型Fig.1 Single-fault of N30°E.

图2 N30°E与 N60°W的交叉断层模型Fig.2 Conjugate faults of N30°E and N60°W.

2种计算模型中，除选取断层系统不同外，其他网格划分以及相关断层设置参数均相同。计算时不考虑重力作用，块体材料假定为线弹性体，杨氏模量E取为70GPa，泊松比为0.25，密度取为2 700kg/m3。2种模型的有限元网格划分均为30 618个节点，25 248个单元(图1，2)。

2.3.2 加载设置

2种断层计算模型加载条件均相同，即最大主压应力轴为EW向，最小主压应力轴为SN向。在计算荷载设置中，分2步加载:1)在模型的西边界面与北边界面上分别加上1个小的应力边界，可使得整个模型得到1个初步的应力均衡;模型东边界的各节点EW向(x方向)自由，SN向(y方向)固定;模型南边界的各节点SN向(y方向)自由，EW向(x方向)固定。2)在模型的西边界面上加载1个x方向的Vref/10大小的相对速度，其他边界条件与第1步中的设置相同。

3 模拟结果与分析

3.1 2种模式下的N30°E断层相对速度云图对比

通过单断层模型与交叉断层模型在相同主压应力作用下的应变速度云图对比，可以形象、直观地观察到先存断层应变速度逐渐增大进而发生突变(即地震孕育和发生的过程)，以及应变速度累积与时间步长的关系。

由单断层模式计算结果可见(图3)，在EW向主压应力作用下，NE向断层在较短的时间步长下很快发生速度突变(图3b)，断层面从闭锁进入到解锁状态(图3c—f)，NE向断层产生滑动。由于应力－应变积累时间短，这种断层模型的速度突变较小，故在此模型下一般只会产生一些中小强度的地震，断层滑动属蠕滑型或黏滑－蠕滑过渡型。

图3 单断层在不同时间步长的相对速度云图Fig.3 Relative velocity of single-fault at different time steps.红色线为断层，速度标尺单位为mm/s a 42.6s;b 43.6s;c 4.0s;d 44.4s;e 44.7s;f 46.5s

图4 交叉断层在不同时间步长的相对速度云图Fig.4 Relative velocity of conjugate faults at different time steps.红色线为断层，速度标尺单位为mm/s a 42.3s;b 43.4s;c 44.7s;d 46.5s;e 358s;f 370s;g 480s;h 594s

交叉断层模型属相互作用的共轭断层系统模型(图4)，同样在EW向主压应力作用下，由于NW向断层与主压应力夹角较小，其断面上正应力分量较小，所受到的摩擦强度较小，容易发生滑动。因此，在较短时间步长内NW向断层从闭锁状态进入解锁状态(图4b—d)，断层面上发生较小的速度突变，进而产生较小应变能的释放，很容易产生蠕滑，故一般发生中小地震;由于NE向断层与主压应力夹角较大，断层面上的正应力分量较大，其摩擦强度大，不易发生滑动，易使应力得到积累达到较高的量值再释放。因此，在一定时间步长的应力积累下，NE向断层从闭锁状态进入解锁状态，断层面上发生较大的速度突变，从而发生一次大的应变能释放，故为发生强震的断层(图4e—h)。因此，在包含NW与NE向2组共轭断层系统的交叉断层模型中，NW向断层经常蠕滑活动，发生一系列中小地震事件，同时为NE向断层应力－应变积累提供了力学条件。

3.2 交叉断层的结点滑动速度曲线对比

以交叉断层为例，在 N30°E和 N60°W断层master面与slaver面上各取6个结点(图5)，观测并讨论不同时刻结点的滑动速度变化。

在较小时间步长下，N60°W断层上结点(图5中的结点1，3，5，7，9，11)发生了较小的速度突变，出现滑动或1次较小的能量释放，即发生一次小震级的地震事件;但经过一定时间的应力积累，N30°E断层(图5中的结点2，4，6，8，10，12)才发生速度突变，且其速度突变较大，即产生了较大的能量释放，对应1次中强或强地震事件的孕育和发生。可见，在区域应力场的作用下，交叉断层模型中2条共轭断层具有交替活动的特征:N60°W断层以释放能量较小的中小地震频繁活动为主，而N30°E断层以释放能量相对较大的中强地震或大地震长间隔复发为主。即与最大主压应力轴夹角较小的N60°W断层面上正应力较小，处于易于错动的方向，失稳事件多;而与最大主压应力轴夹角较大的N30°E断层面上正应力大，处于不易错动的方向，失稳事件少但影响范围大。这一结果与岩石力学实验结果基本一致(马瑾等，2000;马瑾，2009)。

3.3 断层闭锁到解锁过程的初步解释

图5 交叉断层的结点选取Fig.5 Node selection of conjugate faults.

图6 交叉断层不同结点在不同时刻的速度曲线Fig.6 Velocity curves of different nodes of conjugate faults at different time steps.

图7 交叉断层各结点的摩擦剪切应力与极限摩擦强度的关系对比Fig.7 Contrast between shear stress and limiting friction strength of each node of conjugate faults.

图8 交叉断层局部结点的摩擦剪切应力与极限摩擦强度的关系对比放大图Fig.8 Magnified figure of contrast between shear stress and limiting friction strength of local nodes of conjugate faults.

共轭断层系统的相互作用和单条断层闭锁到解锁的应变能释放过程，即地震孕育、发生的力学过程可以用库仑摩擦准则来解释。仍以 N30°E与 N60°W交叉断层为例，在 N30°E和 N60°W 断层 slaver面上各取6个结点(图5)，图7，8为断层面上各点的摩擦剪切应力与极限摩擦强度的关系对比图和局部放大图。例如，黑色方框实点为结点1各加载时刻的摩擦剪切应力f，黑色交叉点为各加载时刻的摩擦强度μFn。在区域应力作用下，断层面的摩擦剪切应力最初小于极限摩擦强度，断层是处于闭锁状态;随着应力作用的持续，断层面上摩擦力在一定时间后超过极限摩擦强度，断层则由闭锁状态突然进入解锁状态，发生速度突变，产生滑动，形成1次地震事件，与岩石力学实验结果一致(Scholz，1998)。例如，华北平原北部有NW向张家口－蓬莱断裂带，是NE向和NW向共轭断层发育区，也是历史地震强烈活动地带，1679年三河－平谷8级地震和1976年河北唐山7.8级地震的发震断层分别为NE向右旋走滑的夏垫新断裂和唐山断裂(表1)，即与主压应力夹角较大的 N30°E主断层在间震期处于闭锁状态，应力易于积累，积累到一定程度时可发生复发间隔较长但震级较大地震的黏滑;而沿张家口－蓬莱断裂带与主压应力夹角较小的N60°W断层易频繁发生中小地震的蠕滑;但NW向断层的存在，为NE向断层的应力积累提供了条件。同样，1966年邢台地震群(M6.8，6.7和7.2)、1976年松潘－平武地震群(M7.2，6.7和7.2)、1988年澜沧－耿马地震(M7.6和7.2)等均发生在共轭断层发育区(Jones et al.，1984;俞维贤等，1994;徐锡伟等，2000)，这些震群的发生无疑与共轭断层相互作用和破裂触发作用有关(Jones et al.，1984;闻学泽等，2006)。因此，采用库仑摩擦准则解释了地震孕育与发生的过程。

表1 华北平原北部主要历史地震Table 1 Major historical strong earthquakes in the northern North China Plain

4 结论

本文针对共轭断层系统这种典型的地震构造模型，采用断层摩擦接触的非线性数值模拟方法，对单断层与交叉断层2种断层模式分别进行了数值模拟，主要取得以下认识:

(1)采用接触摩擦分析单元处理先存断层，可以模拟出地壳介质中存在的不连续间断面;利用断面滑动摩擦机制解释了断层由闭锁状态进入发生地震的滑动状态的过程。

(2)数值模拟所得速度云图结果展示:与主压应力夹角较小的 N60°W断层易频繁发生中小地震的蠕滑，与主压应力夹角较大的 N30°E断层在间震期处于闭锁状态，应力易于积累，积累到一定程度时可发生复发间隔较长但震级较大地震的黏滑;NW向断层的存在，为NE向断层应力积累提供了条件，可用来解释华北地区北部或川滇地区共轭断层系统构成的发震构造模型的滑动习性及其发震危险性的差异。

(3)共轭断层系统的非线性有限元模拟结果，可为地震断层现今运动状态监测工作提供依据。例如，华北平原北部应加强对NE向和NW向共轭断层的监测，并在易于滑动的NW向次断层上捕捉可能的“滑动前兆”，为NE向主断层的大地震预测提供可能的参考依据。

(4)华北平原区和川滇地区地震地质环境复杂，活动断层众多，本文仅对其中共轭断层系统发育区简化了的地震构造模型进行了初步数值模拟探讨。今后应采用摩擦接触非线性有限元方法对整体断层系统做模型的合理简化，进行应力－应变3维数值模拟，尝试性地对大地震发生地点和震级上限进行预测研究，以便有效地减轻可能遭遇的地震地质灾害。