2019, Vol. 62
2. 北京大学地球与空间科学学院, 北京 100871
2. School of Earth and Space Sciences, Peking University, Beijing 100871, China
影响断层破裂过程的因素有很多，复杂的几何形态是其中的一种，它可能导致破裂启动或停止(King and Nábělek, 1985)，从而影响地震波场.分叉形态是断层多种复杂的几何形态中较为常见的一种.过去人们通常将地震断层近似为单个平面，一方面是由于缺乏对断层几何形态的确切描述，另一方面也是由于技术手段的限制.但是，随着观测数据的增多，越来越多的证据表明分叉断层广泛地存在于世界各地，例如1992年美国的Landers地震和2002年美国的Denali地震，就是发震断层包含分叉特征的典型例证.因此，讨论分叉特征对于断层自发破裂过程的影响，可以更好地帮助人们了解断层复杂破裂过程的影响因素，从而对地震震源导致的波场及其造成的危险性有更好的把握.
尽管分叉断层普遍存在，但由于技术手段的限制，直到19世纪80年代中期才出现在震源动力学的研究中(例如King and Nábělek, 1985).技术的发展和观测记录的丰富，以及若干大地震的发生，促进了有关分叉断层自发破裂传播问题的研究.在本世纪初，有关分叉断层的自发破裂传播问题成为震源动力学研究的热点之一.断层的动力学破裂过程数值模拟结果显示，分叉断层的破裂传播，不但与分叉角度和预应力状态有关(Aochi et al., 2000)，也受到破裂前锋速度的影响(Poliakov et al., 2001)，断层的连续性也尤为重要(Aochi et al., 2002).其中，预应力状态决定破裂的最优方向，较高的破裂前锋速度使得非最优方向上的破裂成为可能，较小的分叉角度则降低主断层和分叉断层同时破裂的可能性(Kame et al., 2003; Bhat et al., 2004).
前人的研究，除理论分析之外，大多使用边界积分方程方法进行数值模拟.也有一些研究使用有限元方法(例如，Oglesby et al., 2003; Dreger et al., 2004; DeDontney et al., 2012).但有限元方法不易处理分叉断层三个表面的交线(以下称“三联点”)，甚至可能因为对三联点的不同处理方式得到不同结果(DeDontney et al., 2012).而边界积分方程方法，则不需要对三联点进行额外的处理.并且由于边界积分方程方法积分核计算和破裂过程计算分离的特点，计算量相对较小，节约存储空间.
目前已有的研究有关分叉角度对于断层破裂影响的讨论尚不充分，而且主要限于二维模型.此外，多数研究集中于研究一种特殊的分叉断层：对于两个断层分叉的情况，其中一个为原来主断层的自然延伸，而缺乏对于一般的Y型分叉断层的讨论.Y型分叉断层在自然界中并不罕见，图 1的方框中所示为Y型分叉断层的实例(Hérail et al., 1996).
|
图 1 Y型分叉断层实例(引自Hérail et al., 1996) 这是玻利维亚南部Nazareno盆地至San Vicente断层之间的剖面图，图中方框所示的断层具有Y型的分叉结构. Fig. 1 A true example of Y-shaped branched faults (Hérail et al., 1996) This is a cross section from the Nazareno basin to the San Vicente fault in southern Bolivia. A Y-shaped branched fault is shown in the rectangle. |
本文基于一般的Y型三维分叉断层模型，针对破裂前锋以超剪切速度达到三联点的情况，计算了不同分叉角度的情况，并对结果进行分析.
1 模型、参数及正确性检验 1.1 模型及参数考虑在无限大的三维各向同性均匀弹性介质中的分叉断层.这个分叉断层由共用一条边的三个矩形表面组成，如图 2a所示，共用的这一条边为OO′，其中破裂起始区域(成核区)所在的表面称为第一表面，其余两个表面分别称为第二表面和第三表面.第二表面和第三表面的排名不分先后，在必要的时候完全可以互换它们的名称.将第一表面从OO′一侧延长，第二表面与第一表面延长面所成的夹角为ϕ2，第三表面与第一表面延长面所成的夹角为ϕ3.约定从第一表面的延长面出发，逆时针旋转的角度为正，顺时针旋转的角度为负.
|
图 2 分叉断层示例及其离散化 (a)分叉断层示例：OO′A′A为第一表面；OO′B′B为第二表面；OO′C′C为第三表面；OO′为三联点；ϕ2为第二表面与第一表面延长面的夹角；ϕ3为第三表面与第一表面延长面的夹角.(b)分叉断层离散化示例. Fig. 2 An example of branched faults and corresponding discretized ones (a) An example of branched faults. The first surface: OO′A′A; the second surface: OO′B′B; the third surface: OO′C′C; triple junction: OO′; the angle between the second surface and the extension of the first surface: ϕ2; the angle between the third surface and the extension of the first surface: ϕ3. (b) Discretized branched faults in (a). |
为了进行数值模拟，首先将断层的三个表面进行离散化，分割成若干个大小差距不大的小三角形，如图 2b所示.我们使用相同的时间步长去模拟破裂过程，同时近似地认为，在一个时间步长内，同一个三角形单元内的物理量是相等的.某一时空单元当中沿断层平面的剪应力，在去除了初始应力的因素之后，是历史时间内所有时空单元的滑动速率线性叠加的结果(Kame and Yamashita, 1999a, 1999b, 2003)：
|
(1) |
其中Δτln表示编号为l的三角形在第n个时间步结束时的剪应力变化量，Vik表示编号为i的三角形在第k个时间步结束时的滑动速率，Kτln:ik为积分核，表示第i个单元在k时刻的滑动对第l个单元在n时刻的应力作用系数.
在(1)式中，第n时间步的各三角形的滑动速率同时出现在等式右端，彼此耦合，不便求解.当满足一定条件时，例如Δt≤Δs/(2cP)，该方程可以解耦：
|
(2) |
其中Δt表示时间步长，Δs表示所谓“空间步长”，对于三角形网格来说，则是这些三角形最小内切圆的直径，cP表示P波速度，cS表示S波速度，μ表示剪切模量.方程右边第一项表示共位单元在当前时刻内的应力响应，(2)式是对滑动速率怎样决定剪应力做出的描述.
自发破裂过程的求解还需联立摩擦准则，本文采用如下的滑动弱化准则(Ida, 1972; Palmer and Rice, 1973)：
|
(3) |
其中τ表示剪应力，τr表示剩余剪应力，τp表示剪应力强度，亦即发生破裂所需的最小剪应力，Δu表示滑动量，Dc表示临界滑动弱化位移，H(·)表示阶跃函数.(3)式是对滑动量怎样决定剪应力做出的描述.滑动速率的解，应当使(2)式求得的剪应力与(3)式求得的剪应力保持一致.
联立(2)、(3)两式，可以求出滑动速率、滑动量、剪应力这三个物理量.
由于本文的重点在于研究分叉角度的影响，因此我们采用了无限空间的介质模型.在对分叉断层系统进行非结构化的三角形划分(钱峰等，2019)之后，基于无限空间Green函数计算三角形单元的积分核(Feng and Zhang, 2017)，最后联立摩擦准则和边界积分方程求解自发破裂过程.对于存在地表的介质模型，当断层面与地表的距离超过1 km时，地表的影响可以忽略不计(Oglesby, 2001; Zhang and Chen, 2006)，本文的讨论仍然成立.
本文所涉及参数的取值，列在表 1中.其中，介质密度和波速的取值，位于地壳最常见的花岗岩、片麻岩、石灰岩的参数范围之内；满足
|
表 1 本文所涉及的参数及其取值 Table 1 Parameters involved and their values |
Aochi等(2000)同样使用三维边界积分方程方法研究分叉断层，其方法与本文基本一致，可作为本文计算程序的正确性检验标准.为使参数与Aochi等(2000)保持一致，在1.1节中陈述的一些参数在本小节中做了如下的修改：
断层的几何形状：调整为两个断层表面，第一表面为XY平面上0 km≤X≤6 km，0 km≤Y≤4 km的矩形，第二表面与第一表面在X=4 km处形成分叉，与X轴正向成10°角，向外延伸2 km.
网格剖分：略微增加了网格的大小，使得三角形单元的最小内切圆半径为0.049 km.网格的剖分方式不应显著地影响结果，但因断层面积增加导致三角形单元的数量增加，从而使计算耗时和占用内存增加，因此不得不调整网格大小.由此，本文模型的时间步长与正确性检验标准有一定差别，标准当中给出10倍、30倍、50倍、70倍、90倍时间步长的快照，在本文模型中分别对应于10倍、31倍、51倍、72倍、93倍时间步长.在调整之后，三角形单元的个数为1798个.
成核区位置：仍在第一表面的正中间，X=3 km, Y=2 km.
物理参数：成核区以外介质的初始剪应力值与剩余剪应力的差值修改为τe－τr=75 MPa；临界滑动弱化位移修改为Dc=0.80 m.其他物理参数没有调整.
如图 3，在做出以上调整之后，定性地看，所得结果与正确性检验的标准基本一致；定量地看，存在一些差别，这是因为在Aochi等(2000)中，没有给出参数的具体数值，而是只给出了一些同量纲的物理量的比值，因此以上的参数未必与其数值模拟的真实参数相一致.
|
图 3 分叉断层上的滑动速率快照与Aochi等(2000)结果的比较 下图：Aochi等(2000)计算得到的分叉断层上的滑动速率快照，标出了无量纲时刻，左侧和右侧的矩形分别为第一表面和第二表面；上图：本文对同一模型的数值模拟结果. Fig. 3 Slip rate snapshots on a branched fault compared with Aochi et al. (2000) The lower figure shows the slip rate snapshots on a branched fault computed by Aochi et al. (2000), with the dimensionless time marked, the rectangles on the left representing primary plane and the ones on the right representing off branch. The upper figure shows the numerical simulation result of the same model with our method. |
由于第一表面的破裂情况几乎不受断层的分叉特征影响，以下我们不讨论第一表面的破裂情况；由于第二表面和第三表面的顺序具有对称性，规定我们所讨论的表面是第二表面，并用ϕ2/ϕ3表示.例如，30°/-45°表示，第二表面与第一表面延长面的夹角为30°，第三表面与第一表面延长面的夹角为-45°，以逆时针为正、顺时针为负.
本文的数值算例共计306例，其中有17例ϕ2=－ϕ3的情况，它们的第二表面和第三表面是完全对称的.不考虑第一表面，仍有595个表面；如果把单纯的弯折当作ϕ2=ϕ3的情况，则有612个表面.限于本文的篇幅，不能对这些表面的破裂情况一一叙述，以下先根据模拟结果不同方面的特征，进行分类讨论，再对两个典型的分叉断层序列进行详细的说明.
2.1 分类：是否发生超剪切按表 1所述参数进行模拟，当破裂传播到三联点时，破裂已经进入超剪切状态，破裂前锋速度约为1.19cS.可以预期，在第二表面和第三表面上，超剪切破裂应当是一般情况，而亚剪切破裂状态、破裂终止等其他现象，则为特殊情况.数值模拟结果显示，这种猜测大体上是正确的.612个表面中，超过半数从始至终都是超剪切破裂，图 4a和4b是一个例子，0°/20°始终为超剪切破裂.图中左侧是其断层中间线上滑动速率随时间和位置变化的图像，横坐标为与三联点之间的距离，纵坐标为时间；右侧是断层中间线上破裂前锋慢度随位置变化的图像，图中四条横线分别代表P波和S波慢度的正值和负值.在破裂前锋的慢度图像中，除了左右两端以外(由于网格分布不均匀，慢度的计算不准确)，破裂前锋的慢度总是介于P波慢度和S波慢度之间，说明始终为超剪切破裂.
|
图 4 第二表面破裂过程的多种情况 左列：颜色所示物理量为第二表面中间线上的滑动速率(单位：km·s-1)，横坐标为与三联点之间的距离，纵坐标为时间.右列：折线为第二表面中间线上的破裂前锋慢度与破裂前锋位置的关系，四条水平直线从上到下分别代表S波慢度的正值、P波慢度的正值、P波慢度的负值、S波慢度的负值.(a)和(b) 0°/20°情况，始终为超剪切破裂；(c)和(d)35°/-15°情况，退出超剪切破裂状态后又重新回到超剪切破裂状态; (e)和(f) 40°/45°情况，退出超剪切破裂状态后，没有重新回到超剪切破裂状态；(g)和(h) 30°/20°情况，破裂自然终止；(i)和(j) 85°/65°情况，破裂不从三联点开始. Fig. 4 Different kinds of rupture process on the second surfaces Left column: The slip rate in km·s-1 on the mid-line of the second surface is shown with color. The abscissa shows the distance from triple junction. The ordinate shows time. Right column: The fold line shows how the slowness of rupture front on the mid-line of the second surface varies with the location of rupture front. Four horizontal lines in this figure represent the positive slowness of S wave, the positive slowness of P wave, the negative slowness of P wave, the negative slowness of S wave, respectively. (a) and (b) The rupture is supershear throughout the surface in the case 0°/20°. (c) and (d) The rupture goes subshear and then supershear again in the case 35°/-15°. (e) and (f) The rupture goes subshear until the maximum time step in the case 40°/45°. (g) and (h) The rupture ends by itself in the case 30°/20°. (i) and (j) The rupture does not start at triple junction in the case 85°/65°. |
除此以外，还有以下四种情况：
以图 4c和4d的35°/-15°为例，以超剪切破裂开始，在破裂传播一段距离后退出超剪切破裂状态、进入亚剪切破裂状态，再传播一段距离后重新进入超剪切破裂状态.在破裂前锋的慢度图像中，可以看到，这类情况下破裂前锋的慢度在一定区域内超过了S波慢度，而最终以低于S波慢度结束.这一类情况，常出现于中等角度，即ϕ2≈55°，或第二平面与第三平面的角度之差满足ϕ2－ϕ3≈60°的时候.
以图 4e和4f的40°/45°为例，在破裂传播一段距离后退出超剪切破裂状态、进入亚剪切破裂状态，但是在数值模拟的时间范围内，破裂没有重新进入超剪切破裂状态.在破裂前锋的慢度图像中，这类情况下代表破裂前锋慢度的折线在最右端仍然高于S波慢度，并且由于破裂传播得更慢，常常不能在200个时间步长内抵达断层的末端，因此折线偏短.
以图 4g和4h的30°/20°为例，破裂自然终止.在破裂前锋的慢度图像中，这类情况的破裂前锋慢度折线极短，可能只有几个数据点.
还有极个别情况，某一表面上没有任何破裂发生，整个表面的滑动量都为0.
以上五种情况在612个表面中的分布，展示在图 5中.本文没有考虑正应力的作用.因此，将模型沿xy平面翻转，即同时改变ϕ2和ϕ3的符号，破裂过程没有任何差别.图中把ϕ2 < 0的情况也画了出来，并且加上了ϕ2=ϕ3=0的情况，因此共有1225个点.
|
图 5 不同分叉角度情况的破裂行为 蓝色：始终超剪切；绿色：退出超剪切状态后又恢复到超剪切状态；黄色：退出超剪切状态后维持在亚剪切状态；红色：破裂终止；青色：不发生破裂. Fig. 5 This figure shows how the rupture develop for different branch angles Blue dots mean the rupture is supershear all the time. Green dots mean the rupture goes subshear and then supershear again. Yellow dots mean the rupture goes subshear until the maximum time step. Red dots mean the rupture ends by itself. Cyan dots mean no rupture happens on this surface. |
在612个分叉面当中，超过半数的分叉面是从X=0 km的三角形单元开始破裂的，即从三联点开始，2.1节中的所有例子都是这种情况，不再举例；极个别分叉面，是2.1节所述的没有破裂发生的情况，无法讨论其破裂起始位置.其余分叉面的破裂不是从X=0 km的三角形单元开始、单侧传播，而是从0 km < X < 0.2 km的某个三角形单元开始，然后双侧传播.图 4i和4j所示85°/65°，可以作为这类情况的一个代表.在破裂前锋的慢度图像中，可以看到，这类分叉面的特点是慢度存在负值，这是双侧传播的标志.
图 6展示的是破裂开始点的情况分布.从图中可以看到两类情况的分界线，在上下两端是直线ϕ2=55°，而中间一段较为复杂，无规律可循.
|
图 6 不同分叉角度情况下的破裂开始位置 蓝色：破裂从三联点开始；黄色：破裂不从三联点开始；青色：不发生破裂. Fig. 6 This figure shows where the rupture starts for different branch angles Blue dots mean the rupture starts at triple junction. Yellow dots mean the rupture does not start at triple junction. Cyan dots mean no rupture happens on this surface. |
为了增加图像的对比度、避免一些细节看不清楚，本小节和2.4节的图像都没有调整颜色条、使得各图的颜色条保持一致.
ϕ3=0°这种分叉断层最常见，在既有的研究中也得到了较为充分的讨论，具有特殊意义.因此，这里首先分析这一分叉断层序列.
由图 5和图 6可知，ϕ3=0°时：对于5°≤ϕ2≤50°的情况，破裂自发终止；对于55°≤ϕ2≤70°的情况，破裂在退出超剪切破裂状态、进入亚剪切破裂状态一段时间后，又重新回到超剪切破裂状态；对于75°≤ϕ2≤85°的情况，破裂从X≠0的某处开始，并且其破裂始终在超剪切状态进行.
从断层中间线上滑动速率的图像，如图 7，可以得到更多的信息.在ϕ2较小时，随着ϕ2的增加，发生破裂的区域收缩，滑动速率整体减小，如图 7a—c.但是当ϕ2=40°，第一次破裂终止后又发生了新的滑动，两次滑动在时间和空间上不连续，这是此前的图像中不曾见到的，如图 7d；当ϕ2继续增大，到55°为止，两次滑动彼此连在一起，能够进入超剪切状态而不终止，滑动速率的最大值也有提高，如图 7e—f.随着ϕ2的增加，破裂行为原本的变化趋势发生了改变.
|
图 7 ϕ3=0°情况下破裂过程随ϕ2的变化 Fig. 7 The rupture process varies with ϕ2 when ϕ3=0° |
ϕ2=－ϕ3这一序列的特征是其对称性，因此第二表面和第三表面的破裂行为是完全一致的.由图 5和图 6可知，在这一序列当中：ϕ2≤30°时，破裂直到时间步长最大值为止，仍处于亚剪切状态；35°≤ϕ2≤45°时，破裂始终为超剪切状态；50°≤ϕ2≤65°时，破裂初为超剪切状态，后进入亚剪切状态，最后又回到超剪切状态，并且从ϕ2=55°起，不再由三联点开始破裂；最后，70°≤ϕ2≤85°时，破裂又是始终为超剪切状态.
从断层中间线上滑动速率的图像上看，当ϕ2=5°，有一个自行终止的超剪切破裂过程和一个没有终止的亚剪切破裂，如图 8a；随着ϕ2增大，至20°为止，超剪切破裂过程不断缩短至不可见，而亚剪切破裂则略微提速，如图 8b和8c.ϕ2继续增大，在亚剪切破裂过程之前出现了一段超剪切破裂过程，且随着ϕ2的增大不断延长，直至ϕ2=35°时完全取代亚剪切破裂，如图 8d—f.而当50°≤ϕ2≤65°时，破裂则在中途退出超剪切状态.
|
图 8 ϕ2=－ϕ3情况下破裂过程随ϕ2的变化 Fig. 8 The rupture process varies with ϕ2 when ϕ2=－ϕ3 |
为了解释分叉角度对于破裂过程的影响，我们在计算破裂过程的程序中做了这样的修改：
首先，在判断破裂是否发生的时候，严格禁止第二表面的破裂.这种禁止不是通过增加τp实现的，而是通过人为地让其不发生破裂来实现的，因为在不清楚应力可能出现的最大值的时候，前一种方法可能失效，但后一种方法不会.
而后，在计算第二表面的应力的时候，单独考虑第一表面或第三表面的影响.也就是人为地修改(2)式，有选择地进行求和.
最后把由此计算出来的第二表面的应力记录下来，利用与图 4左列相似的方法，取断层的中间线，记录下这条线上三角形单元应力值与剩余应力的差τ－τr随时间、位置变化的情况.
这种修改方法，有如下的作用：严格禁止第二表面的破裂，可以使得我们对应力所做的这些小动作不致在第二表面引起滑动；由(2)式，第二表面不发生滑动，也就不致影响到第一表面和第三表面的破裂过程.这样，就相当于是一个由第三表面和第一表面组成的弯折断层，而我们能看到第二表面处的应力变化，由此可以尝试对第二表面的破裂过程做出一些解释，尽管这种解释存在一个缺陷，即完全忽略了第二表面自身的破裂对它自身的应力变化产生的影响，但仍优于完全静态的分析.以下将分别对第一表面和第三表面对于第二表面应力变化的贡献进行讨论.
本节涉及到的图片的颜色条范围设为60 MPa≤τ－τr≤100 MPa.因为τ0－τr=80 MPa，80 MPa代表第一表面或第三表面的滑动既不增加第二表面的应力，也不减少第二表面的应力，在图中表示为绿色；红色和黄色的，表示第一表面或第三表面的滑动增加了第二表面的应力；蓝色则表示第一表面或第三表面的滑动减少了第二表面的应力.
3.1 第一表面对第二表面应力变化的贡献第一表面的情况，相对比较简单，因为ϕ3的变化几乎对第一表面没有影响，因此第一表面的破裂过程基本没有差别.尽管如此，我们在挑选图片的时候，仍然控制ϕ3的取值保持一致，取为85°，如图 9.
|
图 9 第一表面对第二表面应力变化贡献的图像 Fig. 9 The contribution to the change of shear stress on the second surface, made by the first surface |
从图中可以看出，当ϕ2=0°，在图像中出现了两条“嵴”，中间夹着一条“谷”：对于断层中间线上的大部分三角形单元，随着时间推移，由第一表面造成的应力变化，由正变负，再由负变正，而负值区域较小.
当ϕ2开始逐渐变大，图像也就随之变化，并且变化逐渐加快.在ϕ2不是特别大的时候，第一表面造成的应力变化的绝对值变小，在图形上表现为“嵴”变低而“谷”变浅.尽管如此，对于ϕ2≤35°的情况，仍然可以看到两条“嵴”和一条“谷”.在图 8a中有一个超剪切破裂区和一个亚剪切破裂区，如果和图 9b对照，是和红黄色区域的“嵴”，也就是受第一表面促进作用的区域一致的，并且把受第一表面抑制作用的区域(即“谷”)夹在中间.这两张图大体上可以相互印证，但还需要3.2节所讨论的两个断层分叉之间的相互作用作为补充.
当ϕ2增大到40°，在X=0 km、t≈80Δt处，应力值开始持续减少，直到ϕ2=85°时，已经可以在这个位置看到比较深的蓝色.当ϕ2≥55°时，破裂基本不再从三联点开始，而且随ϕ2增加不断右移，与此有一定关系.
当ϕ2继续增加，达到约50°，下方红色的“嵴”开始上移，几乎移到了最初蓝色的“谷”所在的位置.这一过程在ϕ2=60°时基本完成.此后，“嵴”开始向右扩展，夹在两条“嵴”之间的“谷”也开始加深.
综合以上图片，直观地来说，对于不同角度的断层分叉，大角度和小角度都是占有优势的，而55°附近的中等角度则占有劣势.如果我们把简单的弯折断层看作是有“单个断层分叉”，则这一点是能够得到验证的.图 5中，ϕ2=ϕ3的点就是这种情况：当40°≤ϕ2≤70°，破裂会在中途退出超剪切状态；对于ϕ2=60°，破裂最终没有回到超剪切状态.ϕ3=0°的情况下破裂的规律也可由此得到一定的解释.
对照正确性检验当中的模拟结果，这一结论的成立是有条件的，与破裂前锋的速度有关，因为用来做正确性检验的0°/10°的情况下，破裂前锋到达三联点处的速度不同，10°表面更占优势.因此，上面这个结论是在当前的超剪切破裂的条件下成立的.
只分析第一表面对第二表面应力变化的贡献，只能部分地解释第2节的数值模拟结果，这种解释是片面的.要完全地解释数值模拟结果，仍然需要分析第三表面对第二表面应力变化的贡献.
3.2 第三表面对第二表面应力变化的贡献第三表面对第二表面应力变化的贡献，明显地受到第三表面破裂情况的影响，在此仅讨论第三表面小角度的情况，用ϕ3=0°代表，如图 10.
|
图 10 小角度第三表面对第二表面应力变化贡献的图像 Fig. 10 The contribution to the change of shear stress on the second surface, made by the third surface when ϕ3 is narrow |
在5°/0°的图像中，可以看到一条直线，直线下方是第三表面对第二表面应力变化没有贡献的绿色区域，上方是使第二表面的应力减小的深蓝色区域.这条直线的方向，与第三表面的破裂过程有关.而绿色区域与深蓝色区域之间的过渡带极其狭窄，这意味着两个断层分叉之间的抑制作用极强，可以在几个时间步长之内从0迅速增加到超过20 MPa.受此影响，在ϕ2－ϕ3=5°的情况下，第二表面和第三表面很难同时发生完全的破裂，甚至在ϕ2=50°, ϕ3=45°和ϕ2=55°, ϕ3=50°这两例算例当中，两个断层分叉受到的第一表面的促进作用相近，因而势均力敌地相互抑制，以致双双破裂终止.
当ϕ2－ϕ3开始增加，这种抑制作用开始减弱.可以看到，当ϕ2－ϕ3增大到25°，原来的深蓝色区域中间，出现浅蓝色；随着ϕ2－ϕ3的继续增加，浅蓝色区域扩大，并且转为绿色，同时原本的深蓝色区域与绿色区域的界限也开始变得模糊不清.ϕ2－ϕ3增大到约40°的时候，新出现的绿色区域当中又出现了红黄色代表促进作用的区域，而深蓝色区域已经压缩到了两条“走廊区域”之内.ϕ2－ϕ3继续增加，达到55°，已经可以认为两个断层分叉之间的相互作用以促进作用为主了，这个抑制作用减弱、促进作用增强的过程，可以对ϕ2=－ϕ3的小角度情况变化规律做出解释.此后，代表促进作用的红黄色区域开始向左移动，蓝色区域开始变浅，而在X=0 km, t=80Δt附近，也出现了一定的促进作用，可以在一定程度上抵消3.1节所提到的第一表面对此处的抑制作用.一些情况下第二表面破裂起始点的转化需要高于55°的ϕ2，可以由此得到解释.
当ϕ2－ϕ3增大到70°时，红色区域的右下边界开始出现一些锯齿.如果与第一表面对第二表面应力变化的贡献进行叠加，则此处会出现一些竖条纹，这可能是数值误差导致的.
对于ϕ3=0°，ϕ2－ϕ3最大只能达到85°；其他情况下，ϕ2－ϕ3可以继续增加，第三表面对第二表面应力变化贡献的变化规律没有变化，仍然是促进作用增强，抑制作用减弱，甚至当ϕ2－ϕ3>160°，断层分叉之间只有相互促进，而不再有相互抑制.
本小节所叙述的“ϕ2－ϕ3很小的时候，会产生强烈的相互抑制作用”，与Kame等(2003)的分析是一致的；而这一夹角较大时可能的相互促进作用，则鲜有提及.
4 结论本文利用边界积分方程方法，对612个不同角度的分叉面的自发破裂传播过程进行了模拟.模拟结果显示，在三联点处的破裂前锋速度超过剪切波速度的情况下，分叉面的破裂过程复杂多样，按有无超剪切破裂和破裂终止可以对破裂过程进行分类，破裂起始点的位置也有差异.通过单独地研究主断层对断层的一个分叉的作用、断层的两个分叉之间的相互作用，数值模拟结果可以得到如下的解释：
分叉断层的一个断层分叉的破裂情况，取决于主断层和另一分叉对其应力变化的贡献.其中，前者与该分叉与主断层延长面的夹角有关，而与另一个分叉关系不大；后者主要与两个分叉之间的夹角有关，但同时也要考虑另一分叉的破裂状况.对于本文所讨论的破裂以超剪切状态到达三联点的情况，前者对于大角度和小角度的促进作用较强，而对中等角度的促进作用较弱，不足以使破裂维持在超剪切状态；后者随着两个分叉之间夹角的增加，对破裂的作用由强烈的抑制转为促进.在两个分叉之间的抑制作用较强的情况下，受主断层促进作用较弱的分叉更易发生破裂终止，这与Kame等(2003)、Aochi等(2000)的结论一致.另外，主断层的破裂到达三联点处的时候，会在较大角度的分叉表面靠近三联点的位置产生一定的抑制作用，这可以对较大角度的断层分叉的破裂大多不从三联点开始这一现象做出解释.
在本文当中，尚没有考虑正应力的作用，没有考虑除了分叉角度以外的其他因素的影响，这将在未来的工作中考虑.
致谢 钱峰、吴葆宁、冯禧三人合作编写了本文所使用的数值模拟程序，钱峰对论文提出了修改建议，在此表示感谢！
Aochi H, Fukuyama E, Matsu'ura M. 2000. Selectivity of spontaneous rupture propagation on a branched fault. Geophysical Research Letters, 27(22): 3635-3638. DOI:10.1029/2000GL011560 |
Aochi H, Madariaga R, Fukuyama E. 2002. Effect of normal stress during rupture propagation along nonplanar faults. Journal of Geophysical Research:Solid Earth, 107(B2): 2038. DOI:10.1029/2001JB000500 |
Bhat H S, Dmowska R, Rice J R, et al. 2004. Dynamic slip transfer from the Denali to Totschunda faults, Alaska:Testing theory for fault branching. Bulletin of the Seismological Society of America, 94(6B): S202-S213. DOI:10.1785/0120040601 |
DeDontney N, Rice J R, Dmowska R. 2012. Finite element modeling of branched ruptures including off-fault plasticity. Bulletin of the Seismological Society of America, 102(2): 541-562. DOI:10.1785/0120110134 |
Dreger D S, Oglesby D D, Harris R, et al. 2004. Kinematic and dynamic rupture models of the November 3, 2002 MW7.9 Denali, Alaska, earthquake. Geophysical Research Letters, 31(4): L04605. DOI:10.1029/2003GL018333 |
Feng X, Zhang H M. 2017. Equivalent formulae of stress Green's functions for a constant slip rate on a triangular fault. Earthquake Science, 30(3): 115-123. DOI:10.1007/s11589-017-0186-3 |
Hérail G, Oller J, Baby P, et al. 1996. Strike-slip faulting, thrusting and related basins in the Cenozoic evolution of the southern branch of the Bolivian Orocline. Tectonophysics, 259(1-3): 201-212. DOI:10.1016/0040-1951(95)00108-5 |
Ida Y. 1972. Cohesive force across the tip of a longitudinal-shear crack and Griffith's specific surface energy. Journal of Geophysical Research, 77(20): 3796-3805. DOI:10.1029/JB077i020p03796 |
Kame N, Yamashita T. 1999a. A new light on arresting mechanism of dynamic earthquake faulting. Geophysical Research Letters, 26(13): 1997-2000. DOI:10.1029/1999GL900410 |
Kame N, Yamashita T. 1999b. Simulation of the spontaneous growth of a dynamic crack without constraints on the crack tip path. Geophysical Journal International, 139(2): 345-358. DOI:10.1046/j.1365-246x.1999.00940.x |
Kame N, Yamashita T. 2003. Dynamic branching, arresting of rupture and the seismic wave radiation in self-chosen crack path modelling. Geophysical Journal International, 155(3): 1042-1050. DOI:10.1111/j.1365-246X.2003.02113.x |
Kame N, Rice J R, Dmowska R. 2003. Effects of prestress state and rupture velocity on dynamic fault branching. Journal of Geophysical Research:Solid Earth, 108(B5): 2265. DOI:10.1029/2002JB002189 |
King G, Nábělek J. 1985. Role of fault bends in the initiation and termination of earthquake rupture. Science, 228(4702): 984-987. DOI:10.1126/science.228.4702.984 |
Oglesby D D, Day S M. 2001. Fault geometry and the dynamics of the 1999 Chi-Chi (Taiwan) earthquake. Bulletin of the Seismological Society of America, 91(5): 1099-1111. |
Oglesby D D, Day S M, Li Y G, et al. 2003. The 1999 Hector Mine earthquake:the dynamics of a branched fault system. Bulletin of the Seismological Society of America, 93(6): 2459-2476. DOI:10.1785/0120030026 |
Palmer A C, Rice J R. 1973. The growth of slip surfaces in the progressive failure of over-consolidated clay. Proceedings of the Royal Society A:Mathematical, Physical and Engineering Sciences, 332(1591): 527-548. DOI:10.1098/rspa.1973.0040 |
Poliakov A N B, Dmowska R, Rice J R. 2001. Dynamic shear rupture interactions with fault bends and off-axis secondary faulting. Journal of Geophysical Research:Solid Earth, 107(B11): 2295. DOI:10.1029/2001JB000572 |
Qian F, Wu B N, Feng X, et al. 2019. 3D numerical simulation of dynamic ruptures on complex fault systems by BIEM with unstructured meshes. Chinese Journal of Geophysics (in Chinese), 62(9): 3421-3431. DOI:10.6038/cjg2019M0642 |
Zhang H M, Chen X F. 2006. Dynamic rupture on a planar fault in three-dimensional halfspace-1.Theory. Geophysical Journal International, 164(3): 633-652. DOI:10.1111/j.1365-246X.2006.02887.x |
钱峰, 吴葆宁, 冯禧, 等. 2019. 基于非结构化网格的边界积分方程方法的断层自发破裂模拟. 地球物理学报, 62(9): 3421-3431. DOI:10.6038/cjg2019M0642 |