弯折断层系统上超剪切破裂传播的产生条件

郑玲珑, 钱峰, 张海明

北京大学地球与空间科学学院地球物理系, 北京 100871

0 引言

地震发生过程中断层的破裂速度是控制地震波形态的重要因素，对于近场的强地面运动有显著的影响.随着观测记录的增多，在越来越多的地震中发现了超剪切传播的证据.已被发现出现超剪切破裂的地震包括1979年Imperial Valley地震(Archuleta, 1984)，1999年Duzce地震(Bouchon et al., 2001)，1999年Kocaeli地震(Bouchon et al., 2001)，2001年昆仑山口西地震(Bouchon and Vallée, 2003)，2002年Denali地震(Dunham and Archuleta, 2004)，2010年玉树地震(Wang and Mori, 2012)，以及2013年Craig地震(Yue et al., 2013)等.在超剪切破裂过程中，断层不同时刻辐射的剪切波互相干涉、叠加，形成马赫波，使得地震波振幅更大，辐射范围更广，往往造成更严重的地震灾害(Bernard and Baumont, 2005; Dunham and Archuleta, 2005; Dunham and Bhat, 2008).多数破坏性的地震往往都发生在几何形态复杂的断层系统上，由于断层空间延展尺度较大，断层的几何形状将会对破裂过程产生重要的影响(马瑾等，1996).研究复杂断层系统上的超剪切破裂转化的条件，对于震源过程的科学认识和防震减灾工作都具有重要的意义.

B-A机制是针对二维模型提出的，认为无量纲比值S临界值为1.77，只有满足一定的背景应力和摩擦条件时，才能够产生超剪切破裂(Andrews, 1985).针对三维均匀介质模型，Madariaga和Olsen(2000)提出了另一种无量纲化参数κ作为确定破裂速度的指标，定义为破裂应力释放能与断裂能之比，包含了更多破裂信息，适用于高应力场断层破裂情况；而S则需要小于更小的临界值1.19，才能触发超剪切破裂(Dunham, 2007).考虑包含地表的均匀半空间介质，S值大于1.19也产生超剪切破裂(Bizzarri, 2010; Bizzarri and Das, 2012; Zhang and Chen, 2006).Olsen等(1997)、Kaneko和Lapusta(2010)认为自由界面SV-P衍射波的相变转换所产生的附加应力场是超剪切破裂形成的主要原因，而走滑断层自由表面所诱发的超剪切破裂也可能会返回至次瑞利速度(Hu et al., 2019).即使不考虑自由界面的效应，三维全空间断层介质的非均匀性使得在不满足B-A机制所预测的背景应力水平下仍然有可能出现超剪切破裂(Passelègue et al., 2013)，但也可能使得当背景预应力超过所预测的临界水平时仍然以亚瑞利波速稳定传播，无法出现超剪切破裂(Fukuyama and Madariaga, 2000).

影响超剪切破裂出现的复杂因素还包括，非成核区介质摩擦性质和背景预应力分布的不均匀性，以及断层几何形状的复杂性等.一般以包含障碍体(barrier)和凹凸体(asperity)的断层面的粗糙度和准静态裂缝表示非成核区介质摩擦性质和背景预应力的非均匀性.Fukuyama和Olesn(2002)研究了不同凹凸体宽度和不同背景预应力对超剪切破裂的影响，发现在一定条件下，含有凹凸体的非均匀介质相对于均匀介质，更容易激发超剪切破裂；对于宽度较小的凹凸体需要更高的背景预应力来触发超剪切破裂，当凹凸体宽度大于某一临界值后，所需要的背景预应力不再降低，保持稳定后仍然可以产生超剪切破裂.当非成核区存在准静态二次扩张裂缝时，根据B-A机制，二次裂缝将会受到主破裂前S波和P波应力场的影响开始破裂，当主破裂前锋与二次裂缝前端相遇时速度会突然超过剪切波速，转换成超剪切破裂(Liu and Lapusta, 2008).针对高剪切强度的障碍体，Dunham等(2003)提出在障碍体与成核区之间具有一定距离时，破裂前锋绕过障碍体，其产生的凹形聚焦效应能够使破裂传播速度迅速增加到超剪切状态.之后的研究发现障碍体诱导形成超剪切破裂需要满足一定的条件(Duan, 2012; Yang et al., 2012, 2013; Weng et al., 2015)，障碍体埋深越浅、尺寸越大和强度越高, 越容易产生超剪切破裂，当障碍体位于自由界面，与成核区在自由界面上的距离大于一定长度时，能够扩大超剪切破裂的传播距离(Xu et al., 2016).当断层面本身大面积凹凸不平，粗糙度越大越容易产生局部超剪切破裂，但持续的超剪切往往出现在粗糙度相对中等的位置，过大的粗糙度会损耗能量，使超剪切破裂无法持续(Bruhat et al., 2016).

断层的几何形状是影响超剪切破裂产生的另一个重要因素.目前涉及到断层几何形状的研究主要集中于几种端元的复杂断层模型：阶跃断层、弯折断层和分叉断层.关于阶跃断层上超剪切破裂的影响因素研究较为完善，当满足一定范围内的断层间隔和足够长度的平行重叠距离后，阶跃位置的S值会发生变化降低到临界值之下，积蓄能量在次级断层上产生超剪切破裂(Kase and Kuge, 1998; Ryan and Oglesby, 2014; Hu et al., 2016; 袁杰和朱守彪，2017).因为弯折和分叉断层系统的复杂性，在数值模拟中相对较难实现，故而对其产生超剪切破裂的研究更加困难.Aochi等(2000a,b)基于边界积分方程法(BIEM)对简单弯折和分叉断层进行了数值模拟，发现当光滑弯曲断层面的最终角度为π/18时，破裂能够以超剪切速度传播，初始剪切应力和摩擦参数之间的关系对于简单分叉断层系统动态破裂传播的断层选择是最重要的，法向应力的变化通常较小.不同的弯折角度对破裂的传播和停止也有着重要影响(张丽芬等，2016).Kame等(2003)通过边界积分方法研究了二维分叉断层破裂速度影响条件，发现只有当分叉断层的S值从3突降到临界值以下的0.8时，才能在断层上产生超剪切破裂.因此对于弯折和分叉断层来说，B-A机制所提出的影响因素仍占据主导地位，针对不同形状的弯曲断层可能有着不同的临界状态.对于分叉断层系统，在一定条件下分叉间的相互作用也有可能使其中一支分叉断层产生超剪切破裂(Bhat et al., 2004).

目前有关复杂断层系统上的超剪切破裂的研究较少.本文将基于三维空间中非结构化的BIEM，针对不同弯折角度的弯折断层进行大量高效的数值模拟，在参数空间中通过破裂前锋传播图和破裂速度相图研究弯折断层的破裂速度变化，分析弯折断层的破裂行为和超剪切破裂的产生条件.

1 自发破裂模拟方法及破裂相图算法

1.1 边界积分方程方法(BIEM)

本研究采用的自发破裂数值模拟方法为BIEM，它建立了断层应力与滑动速率的关系(Das and Aki, 1977; Aochi et al., 2000a,b)：

(1)

(2)

对于断层自发破裂过程而言，必须结合控制破裂行为的摩擦准则.常用的摩擦准则有滑动弱化摩擦准则(Ida, 1972)和速率-状态摩擦准则(Dieterich, 1979).对于同震破裂过程的模拟，通常采用滑动弱化摩擦准则：

(3)

其中，τ是当前时刻单元的剪切应力值，τu和τf分别是破裂强度和残余应力，Dc为临界滑动弱化位移，D为当前时刻的单元累计滑动量.当断层单元上的剪切应力累积到了破裂强度时，单元开始破裂，其剪切应力随着滑动量线性减小，直至达到稳定滑动后保持不变的摩擦力.将BIEM和滑动弱化准则联立，给定初始应力后就可计算得到破裂过程.

1.2 破裂速度自动提取

完整破裂相图的参数空间包含几千个破裂结果，需要逐个判断每个算例的破裂速度.为了提高效率，本文采用了适用于全空间均匀走滑断层模型的破裂速度自动提取方法，能够在破裂传播计算过程中得到算例的破裂速度，步骤如下：

图1 破裂速度拾取示意图 (a) 平面断层模型，黄色高亮单元为破裂前锋中心点传播路径的网格单元, 绿色单元为成核区单元. (b)平面断层滑动速率时空分布，图中VP和VS两条白线分别代表P波和S波的传播速度，红圈为破裂前锋随时间变化的x方向位置，绿 线为红圈的拟合线条，其斜率近似反映破裂的传播速度.Fig.1 Illustration of the rupture speed picking (a) The planar fault model. The highlighted yellow elements are the grids on the propagation path of the midpoint of the rupture front, the green element is a grid inside the nucleation zone. (b) The spatiotemporal evolution of the slip velocity of the planar fault model. The white lines denote the reference P and S wave velocity respectively. The red circles are the locations of rupture fronts along the fault varying with time, whereas the slope of the fitted connecting green line approximates the rupture propagation speed.

(1)图1a中显示为绿色单元的成核区(初始破裂区域)近似为圆心坐标为(22.5,7.5)，半径为3 km的圆，根据成核区的位置，标记破裂前锋的中心点在断层上的传播路径中黄线所示，对路径上单元按照x方向坐标从小到大排列；

(2)在算例计算过程的每个时间步循环内，对步骤1中已排序的传播路径单元的滑动速率值进行自动筛选，找到成核区沿传播方向上滑动速率降为0的最小x坐标单元，即破裂前锋对应单元，并记录时间及该单元x方向坐标，用图1b中的红圈标记；

(3)每隔一定时间间隔对步骤2中记录单元的坐标和时间进行x-t线性拟合，拟合结果可示意为图1b中红圈拟合的绿线，计算拟合线的斜率得到破裂传播速度，当破裂的时间循环计算结束时就可自动得到一系列实时变化的破裂速度.

1.3 平面断层的相图对比

为了检验破裂速度自动拾取方案的正确性，我们使用与Xu等(2015)近似相同的参数(表1)和相同的无量纲化公式(4)计算得到平面断层破裂相图：

(4)

图2显示了本文与前人的平面断层破裂相图，对比发现平面断层破裂相图整体趋势相同，验证了破裂过程自动拾取的正确性.

图2 平面断层破裂相图 左(本文)，右(Xu et al., 2015).0～0.2区间内过小，导致计算结果不可靠，故整体在破裂相图中以灰色表示， Supershear区域表示算例为超剪切破裂， SubRayleigh区域表示算例破裂速度为亚瑞利波速，Self-arresting区域表示算例破裂未持续传播自发停止.Fig.2 Phase diagram of rupture on a planar fault (Left) Result in this study, and (right) that in Xu et al. (2015). Since the small value in the range 0～0.2 would produce unreliable results, this part of parameter space is shaded by gray. Region in which super-shear ruptures occur is marked in “Supershear”. Rupture speed below Rayleigh wave speed, “SubRayleigh”, and rupture ceases spontaneously, “Self-arresting”.

表1 平面断层模型参数Table 1 Parameters of the planar fault models

2 弯折断层破裂模拟结果与讨论

利用前面所描述的方法，本文分别对弯折角度为15°、25° 和40° 的弯折断层模型进行破裂相图计算.弯折断层模型如图3所示，断层模型从x方向45 km处开始弯折，整个断层模型由平面部分15 km×45 km的矩形平面和弯折平面部分为15 km×15 km的正方形平面以弯折角α拼接而成，成核区位于矩形平面的中央.

图3 弯折角α=15°的弯折断层模型 黄色三角形单元为图6所取弯折面上靠近弯折交界线单元.Fig.3 The fault model with a bending angle α 15° The highlighted yellow triangle element is the grid on the bending plane chosen in Fig.6 that is close to the bending boundary.

(5)

其中：τa=σ33(μs-μd)，σ33为弯折断层平面部分初始正应力.

表2 弯折断层模型参数Table 2 Parameters of the fault models with a bending angle

观察三个破裂相图的自发停止破裂、亚瑞利破裂和超剪切破裂之间的两条分界线I和II可以看出，自发停止破裂区域范围在三个不同角度的弯折断层相图中完全一致，不受弯折角度改变影响；而三个相图中分界线II的不同，则显示当破裂能够稳定传播时，弯折面的角度对破裂速度有明显影响.随着弯折角度增加，超剪切区域范围减小，亚瑞利破裂范围增大，15° 与25° 弯折断层破裂相图的差异较小，弯折角度为40° 时亚瑞利波破裂范围明显增加，且超剪切破裂边界线变化趋势与平面破裂相图较为相似.

图4 弯折断层破裂相图对比图 (a)、(b)、(c)分别为弯折角α=15°，25° 和40° 的结果.绿色虚线I为自发停止破裂和亚瑞利破裂之间的分界线，蓝色虚线II为亚瑞利破裂和超剪切破裂之间的分界线，Sub-Rayleigh 1代表算例破裂速度为亚瑞利波速且破裂未传播至弯折面，Sub-Rayleigh 2代表破裂速度 接近但不超过剪切波速的算例，图中A1—I6所圈参数的算例依次对应图6相同编号的子图.Fig.4 Comparison of phase diagrams of rupture on faults with different bending angles (a), (b) and (c) are results with the bending angle of 15°, 25° and 40°, respectively. The dashed green line I separates the Self-arresting and SubRayleigh. The dashed blue line II separates the SubRayleigh and Supershear. Examples marked in “SubRayleigh 1” represent ruptures whose speed are below the Rayleigh wave speed. Examples marked in “SubRayleigh 2” represent ruptures whose speed are approaching, but not beyond the shear wave velocity. The examples circled and annotated with A1—I6 correspond to each subplot in order in Fig.6.

图5a、b、c分别为平面初始正应力(σ33,σ11)设置为(-10, -26) MPa，(-10, -15) MPa, (-10, -10) MPa的不同弯折面应力莫尔圆示意图，每张子图展示了给定同一平面初始应力时，不同弯折角的断层面初始应力和S.可以看出初始状态下对于15°、25°和40°弯折断层面，α越小的断层面越靠近破裂线τp，参数S越小，意味着破裂越容易传播，也就越容易产生超剪切破裂.图5a为图4相图初始应力设置，与图4的弯折角越小超剪切破裂区域越大规律一致.三幅子图中图5a即本文所采用应力设置，初始状态下各个弯折面对应S值最小，作为对照组的图5b、c的应力设置使得平面和弯折断层在莫尔圆上位置整体左移，逐渐偏离破裂线.

图5 断层面不同初始正应力莫尔圆示意图 (a)、(b)、(c)分别为平面初始正应力(σ33,σ11)设置为(-10, -26) MPa，(-10, -15) MPa, (-10, -10) MPa的不同弯折面应力莫尔圆示意图，三个弯折模型除弯折角α外的其他参数相同，所取为0.85，S=(τp-τo)/(τo-τr)，τp为破裂强度，τo为初始剪应力，τr为剩余剪应力，σn为正应力.图中弯折角α=0°的灰色正方形所标记的横纵坐标展示平面断层的初始正应力和剪应力， 其他形状标记对应了 不同弯折角的弯折面初始正应力和剪应力.Fig.5 Diagram of Mohr′s circle under different initial normal stresses on the fault model The planar initial normal stress (σ33,σ11) for the Mohr′s circle diagram (a), (b) and (c) are set to be (-10, -26) MPa, (-10, -15) MPa and (-10, -10) MPa, respectively. All parameters but the bending angle α are the same for three models. is set to be 0.85 and the ratio S=(τp-τo)/(τo-τr), where τp is the rupture strength, τo is the initial shear stress, τr is the residual shear stress, and σn is the normal stress. The case α=0° is denoted by the gray square, whose horizontal and vertical coordinate represent the initial normal stress and shear stress on the planar part, respectively. Other symbols correspond to the initial normal stress and shear stress on the bending part， with different bending angles.

图6 15°、25°和40° 弯折断层分别取0.35, 0.70和0.93时的滑动速率时空分布组图 每个子图代表一种参数组合下结果，子图左下角所标序号对应着图4相应相图中所标的序号，红圈代表破裂前锋，为了方便判断破裂 前锋的速度，标出了P波和S波速度.破裂算例对应参数以标注在白色括号内，白色虚线为弯折交界线x=45 km的标记.Fig.6 The selected spatiotemporal evolution of slip velocity in 15°-, 25°- and 40°- bending fault when is 0.35, 0.70 and 0.93, respectively Each subplot represents the result under specific combination of parameters which is marked in white brackets as The sequential number at the lower left corner of each subplot is the same as that in Fig.4. Red circles denote the rupture front. P and S wave velocities are illustrated for reference. The white dashed line marks the location of bending boundary (x=45 km).

图7 弯折断层相同初始正应力的滑动速率时空分布对比图 (a)、(b)、(c)分别为平面初始正应力(σ33,σ11)设置为(-10, -15) MPa的弯折角α=15°，25° 和40° 的弯折断层滑动速率时空分布图，子图所取为(0.70, 0.85)，图中VP和VS两条白线分别代表P波和S波的传播速度，红圈示意为破裂前锋随时间的变化 趋势，绿线为红圈的拟合线条，其斜率近似反映破裂的传播速度.Fig.7 Comparison of the spatiotemporal evolution of slip velocity under the same initial normal stresses on the bending fault model The planar initial normal stress (σ33,σ11) for the spatiotemporal evolution of slip velocity on the model with a bending angle of 15°, 25° and 40° are set to be (-10,-15) MPa, respectively. The parameters for the subfigures are (0.70,0.85), respectively. The white lines denote the reference P and S wave velocity, respectively. The red circles are the locations of rupture fronts along the fault varying with time, whereas the slope of the fitted connecting green line approximates the rupture propagation speed.

图8 弯折面上交界线附近一个单元的不同弯折角度断层模型应力时间演化对比图 (a)、(b)、(c)分别为平面初始正应力(σ33,σ11)设置为(-10,-26) MPa，(-10,-15) MPa, (-10,-10) MPa的不同弯折角度断层模型应力时间演化对比图.每个子图中三个弯折模型算例除弯折角度外的其他参数相同，为(0.70, 0.85)，Te为 剪切应力，Tu为破裂强度，N为正应力.Fig.8 Comparison of the temporal variation of stress at one element grid close to the bending line on the fault model with different bending angles (a), (b) and (c) correspond to the comparison of the temporal variation of stress on the fault model with different bending angles when the planar initial normal stress (σ33,σ11) are set to be (-10,-26) MPa, (-10,-15) MPa and (-10,-10) MPa, respectively. All parameters but bending angles are the same for three subfigures. is (0.70,0.85), where Te,Tu and N are the shear stress, the fracture strength and the normal stress, respectively.

3 结论

对不同弯折角度的相图和选定参数组合下的滑动速率时空图分析，可以发现对于本文所设应力条件下弯折面处于压缩区(α>0)的弯折断层模型，有以下结论：

(1)弯折角度对于无法持续传播的自发停止破裂没有影响，相图中自发停止破裂区域所占参数空间完全相同；

(3)体现在相图上，15° 和25° 的超剪切区域明显比40° 的大，即破裂能够越过弯折交界线在弯折面上继续传播时，更容易形成超剪切破裂； 15° 的超剪切区域又略大于25°，即弯折角越小，越容易形成超剪切破裂；

(4)对初始剪应力较小的情况，超剪切破裂更容易在弯折面上形成，而当初始剪应力较大时，超剪切破裂在平面上就可形成；

根据滑动速率时空分布图，可以直观地分析破裂传播速度和滑动速率的变化，当破裂传播至弯折面上时，需要克服弯折交界线上的正应力，从而破裂前锋在交界线处出现错位延迟效应，同时在弯折面上的滑动速率相对于平面断层有了明显的增大.弯折角度越大，需要克服的时间越长，延迟效应越明显.

本文考虑的全空间弯折断层模型，包括平面部分和相对弯折部分，破裂成核区设在平面断层，弯折断层部分位于平面断层的压缩区，这只是比平面断层稍复杂一点的断层系统，本文进行了该类断层系统的破裂行为研究的初步尝试.对于其他类型的复杂断层，比如分叉断层，甚至更为复杂的断层系统，其破裂行为的控制因素将更为复杂，这将在进一步的研究中予以考虑.