2013年4月20日四川省雅安市芦山县发生了M_S7.0地震，这是继汶川地震以来，在龙门山断裂带上的又一次强震，引起了严重的山体崩塌、滑坡、泥石流、砂土液化等次生灾害，造成了重大的人员伤亡和财产损失(Zhang et al., 2013).震后野外应急科考和浅层人工地震勘探调查发现，芦山地震极震区和附近断裂带沿线均没有出现构造成因的地震地表破裂带，推断此次地震发震断层属于隐伏断层(李传友等, 2013; 徐锡伟等, 2013; 雷生学等, 2014).不同研究者采用近震、远震数据反演芦山地震震源机制解，结果表明，主震发震断层走向约210°，滑动角近90°，可视为纯逆冲型地震；由于使用的地震数据或滤波频段不同，断层倾角和发震深度结果有较大差异，整体而言，断层倾角约为42°，主震深度约15 km(刘杰等, 2013; 吕坚等, 2013；谢祖军等, 2013; 许力生等，2013; 曾祥方等, 2013).芦山地震余震重定位结果显示，余震主要分布在10~20 km深度上，浅部余震稀少，在垂直于龙门山断裂带的剖面上，余震集中在倾向相反、相互交叉的两个条带上，构成“ Y ”字型分布(Fang et al., 2013; Long et al., 2015; 赵荣涛等, 2015; Lu et al., 2017).此外，深地震反射剖面进一步显示了芦山地震区复杂的地壳结构和断裂构造特征，震区地壳厚度约40 km，上部地壳分层特征明显，褶皱和断裂构造清晰.由于受到北西-南东向挤压应力的作用，上部地壳形成了一系列倾角上陡下缓的铲形断层.结合余震重定位和震源机制解等研究结果，推断芦山地震发震断层位于新开店断层和大邑断层之间，伴有反冲断层，断层系统呈“ Y ”字型分布，深部归并到约16 km左右的滑脱面上，断层具有一定埋深(王夫运等, 2015; 冯杨洋等, 2016).综合地表地质调查、余震定位、震源机制解以及深地震反射剖面等结果，推断芦山地震可能是一次发生在复杂断层系统中的盲逆破裂事件(徐锡伟等, 2013; 王夫运等, 2015).

这次地震的震源破裂过程是受到广泛关注的问题.不同的研究者利用地震波、GPS和InSAR等观测资料，采用有限断层模型对芦山地震的破裂过程和同震滑动分布进行了反演研究(Hao et al., 2013; Jiang et al., 2013; Liu et al., 2013; 王卫民等, 2013; 张勇等, 2013; Zhao et al., 2013; Zhang et al., 2014; 刘琦等, 2016; 郑绪君等, 2018).结果显示，此次地震以单个破裂事件为主，破裂方向性不明显；破裂持续约10 s；同震滑动主要集中在震源附近，分布在断层走向方向约20 km范围内，最大滑动量约为1.3 m，断层浅部没有明显滑动量分布，破裂没有突破地表；整体破裂方向接近90°，此次地震可视为纯逆冲型地震；地震震级约M_W6.6；部分结果显示震源东北方向最终滑动量较小(Hao et al., 2013; Zhang et al., 2014; 郑绪君等, 2018).而GPS、InSAR等反演手段可能受发震断层模型影响较大，得到的平均最大滑动量约为0.8 m(Jiang et al., 2013; 谭凯等, 2015; 刘琦等, 2016).

目前，针对芦山地震破裂过程和同震滑动分布的研究多通过波形拟合等反演手段进行，没有考虑导致破裂过程出现的原因——应力因素.为了从力学角度解释此次地震震源破裂过程的原因，需要进行震源动力学研究.本文基于芦山地震野外地质调查、震源机制解反演、余震定位、深地震反射剖面等结果，构建芦山地震发震断层模型，在此基础上，结合岩石力学实验得到的结果，以震源运动学反演结果作为约束，通过试错法调整震源动力学计算参数，对芦山地震进行大量动力学模拟研究，得到芦山地震破裂传播的可能情况，进而分析讨论不同动力学计算参数对芦山地震破裂过程和同震滑动分布的影响.

目前对于芦山地震破裂过程和同震滑动分布的研究，大多将发震断层简化为具有固定倾角的单一平面(Hao et al., 2013; Liu et al., 2013; 王卫民等, 2013; Zhang et al., 2014; 郑绪君等, 2018).但是芦山地震余震定位、深地震反射剖面等结果表明，此次地震实际的发震断层可能是一个上陡下缓的铲型断层，浅部连接近水平的滑脱带，几何结构比较复杂(Fang et al., 2013; Long et al., 2015; 王夫运等, 2015; 冯杨洋等, 2016; Lu et al., 2017).由于断层的几何结构对地震破裂传播有着重要影响(Aochi et al., 2000a, b; Aochi and Fukuyama, 2002)，因此本文结合芦山地震余震定位、深地震反射剖面等结果，构建芦山地震三维铲型断层模型，该模型包括NW倾向的主断层及其连接的深度约16.6 km的基底滑脱带，如图 1所示.为了更好地刻画铲型，本研究采用非结构化的三角形网格离散方案(钱峰等, 2019)，在计算量允许的情况下对模型进行离散，如图 1所示.断层模型走向方向长30 km，倾向方向长22 km，共1640个三角形单元，三角形网格的尺度Δs=0.245 km.

图 1

Fig. 1

图 1 芦山地震断层模型 x表示断层倾向方向，y表示断层走向方向，z表示深度. Fig. 1 The fault model of the Lushan earthquake x is the direction of fault inclination, y is the direction of fault strike, and z is the depth.

本研究采用能够灵活处理非平面断层破裂问题的边界积分方程法(BIEM)，对芦山地震破裂过程进行数值模拟.由于芦山地震发震断层上断点埋深大于1 km，自由表面对破裂过程的影响可以忽略不计(Zhang and Chen, 2006)，因此可以采用基于全空间模型的BIEM.这种方法利用基本解建立相应的边界积分方程，将域问题转化为边界问题，因此不仅降低了求解维度，而且便于求解具有复杂几何形状的断层破裂问题，比如弯折断层、铲型断层等.由于仅与断层几何形状有关的积分核计算与破裂过程的计算是分离的，因此非常适合几何形态固定，需要通过试错法探究动力学参数对破裂过程影响的研究.

BIEM基于格林函数求解弹性波动方程，建立断层上的牵引力与滑动速率函数之间的关系(Fukuyama and Madariaga, 1998; Tada et al., 2000)：

(1)

其中，T(x, t)表示接收点x在t时刻的剪切牵引力，T⁰表示初始剪切牵引力，表示源点ξ在τ时刻的滑动速率函数，表示源点ξ处在初始时刻一个单位速率的滑动在x处t-τ时刻产生的应力响应. 是包含格林函数导数的积分核函数，反映了单元与单元间的相互影响的程度，只与介质参数和单元几何信息有关，不随破裂过程而变化.

为了求解自发破裂问题，需要结合一定的摩擦本构关系.摩擦本构关系决定了破裂行为，控制了破裂过程.在同震破裂计算问题中，通常采用滑动弱化准则(Ida, 1972)，其表示形式为：

(2)

其中，T(D)表示滑动量为D时的剪切应力，T_u=σ_u-σ_f表示应力降，σ_u表示破裂强度，σ_f表示剩余应力，D_c表示临界滑动弱化位移，H(D_c-D)表示阶跃函数，当D_c>D时值为1.对于断层面上的某一点，在构造应力场作用下，应力逐渐积累，直到剪应力值超过σ_u发生滑动.随后剪应力随着滑动量的增加而线性地减小；当D超过D_c后，剪应力保持为恒定的常数σ_f.

结合边界积分方程和滑动弱化准则，定义初始应力分布、临界滑动位移分布、成核区分布、破裂强度分布等动力学参数后就可以进行破裂过程模拟.

对动力学破裂过程计算的结果依赖于多个动力学参数的共同作用，因此有必要谨慎设置每一个计算参数.

首先是区域初始应力场.由于芦山地震可以视为纯逆冲型地震，因此本文不考虑沿断层走向方向的应力.虽然很少有研究直接给出地下深部应力情况，但是由于温度和围压的变化，应力在很大程度上与深度有关(如, Aochi and Madariaga, 2003; Zhang et al., 2019)，因此本文考虑随深度变化的应力设置.为了使破裂能够在震源附近成核传播，需要优先考虑震源所在深度的应力分布情况.假设参考平面为同一深度下的水平面，其受到竖直方向应力σ_V和水平方向应力σ_H、τ_H，水平方向应力垂直于断层走向方向，沿着x方向为正方向.根据子断层单元倾角的大小，通过应力变换可以计算出子断层单元上的初始正应力σ与剪应力τ.由于芦山地震动力学模拟受运动学反演结果约束，因此当进行试错实验时得到的最大滑动量与运动学反演结果基本符合时，即可得到合适的区域初始应力设置.本文选取的应力设置如图 2a所示.

图 2

Fig. 2

图 2 动力学计算参数随深度的变化 (a)区域初始应力:σV是竖直方向应力，σH、τH为水平方向应力; (b)临界滑动弱化位移Dc. Fig. 2 Variation of dynamic parameters with depth (a) Regional initial stress:σV is vertical stress, σH and τH are horizontal stress; (b) Critical slip weakening displacement Dc.

断层的破裂强度一般是正应力与静摩擦系数的乘积.当正应力已知，设置静摩擦系数即可计算得到破裂强度，本文中静摩擦系数μ_s取0.6，动摩擦系数μ_d取0.1(如, Aochi and Madariaga, 2003; Zhang et al., 2019).当剪应力超过破裂强度时，断层开始发生滑动.

对于临界滑动弱化位移D_c，大量的实验室实验以及野外观测资料分析表明，在不同地震中D_c取值自由度较大.通常在脆性岩层中D_c为一常数，在地壳浅部和韧性岩层中D_c较大(Ohnaka，1992；张丽芬和姚运生, 2013).在实际的地震动力学模拟中，还应根据具体情况不断调整D_c取值，直至破裂与运动学反演结果大致符合.本文假设D_c随深度变化(如, Aochi and Fukuyama, 2002; Zhang et al., 2019)，满足公式(3) :

(3)

其中h₁为4 km，h₂为16 km，d₀为0.3 m，见图 2b.h₁、h₂的取值是参考余震重定位时采用的速度模型设置的(Fang et al., 2013; 赵荣涛等, 2015).

为了使破裂在开始时刻被触发，在当前网格条件下，成核区半径小于3 km无法孕育破裂，因此设置一个半径为R=3 km的成核区，对成核区内部初始剪应力T_i重新赋值，使其略大于破裂强度T_u，也就是假定成核区已发生破裂T_i=1.01T_u.

在动力学破裂模拟中，若没有其他障碍体的情况下，破裂永不停止.但是根据动力学反演结果来看(如, Hao et al., 2013; Liu et al., 2013; Zhang et al., 2014; 郑绪君等, 2018)，滑动量主要分布在震源附近约20 km范围内，即破裂被限制在一定区域内.因此本研究假设震源20 km范围以外存在一近圆形的障碍体，使得断层的破裂强度是原来的2倍，当破裂接触到障碍体后，会自动减慢甚至停止，以达到限制破裂区域的目的.

由于本研究主要分析初始应力分布、临界滑动位移分布、成核区分布、破裂强度分布对破裂过程的影响，其余参数均采用固定值，为降低复杂度，本文采用均匀介质模型，如表 1所示，其中V_P, V_S, ρ分别为介质的P波、S波速度以及密度, CFL为稳定性系数(Fukuyama and Madariaga, 1998)，Δs为计算步长，Δt为时间步长，Step为计算步数，H为震源深度.

表 1 (Table 1)

表 1 介质参数与断层模型参数设置 Table 1 Parameters of medium and fault model α/(km·s-1) β/(km·s-1) ρ/(kg·m-3) CFL 5.85 3.45 2700 0.98 Δ s/km Δ t/s Step H/km 0.245 0.041 300 15

表 1 介质参数与断层模型参数设置 Table 1 Parameters of medium and fault model

为了分析初始应力分布、临界滑动位移分布、成核区分布、破裂强度分布对破裂过程和最终滑动量分布影响，本文采用控制变量法(张丽芬等，2016)进行大量数值模拟实验，通过对比模拟的结果探究不同震源动力学参数对芦山地震破裂过程和最终滑动量分布的影响.

为了研究断层在不同初始应力状态下的破裂情况，在不改变其他动力学参数的情况下，选取三组不同的初始应力设置，以震源深度15 km为例，图 3显示了A、B和C三种应力状态在莫尔圆中的位置，横坐标代表断层面上正应力，纵坐标代表断层面上剪应力，τ_p=μ_sσ代表破裂线.其中，A状态代表图 2a区域应力设置中震源位置的受力情况.A状态下震源处所受正应力为7.53 MPa，剪应力为3.2 MPa，B状态下震源处所受正应力为7.71 MPa，剪应力为3.59 MPa，C状态下震源处所受正应力为7.13 MPa，剪应力为2.89 MPa.从图 3可以看出，C点距离破裂线最远，B点距离破裂线最近.

图 3

Fig. 3

图 3 初始应力分布 圆表示应力莫尔圆，τp=μsσ和τr=μdσ分别表示破裂强度线和残余应力线，A、B和C表示三种不同应力状态，σ、τ分别表示正应力和剪应力，μs、μd分别表示静摩擦系数和动摩擦系数. Fig. 3 Initial stress distribution Circle is the Mohr circle of stress. τp=μsσ and τr=μdσ are the breaking strength line and the residual stress line, respectively. A, B and C represent three different stress states. σ and τ are the normal and shear stresses, μs and μd are the static and kinematic friction coefficients, respectively.

图 4展示了不同初始应力状态下的滑动速率快照(a-c)和最终滑动量分布图(d-f).可以看到，在A、B应力状态下，破裂在成核区孕育，随着剪应力的逐渐积累，破裂沿断层走向方向呈椭圆形向四周传播，在接触到障碍体后，破裂逐渐停止；而在C应力状态下，剪应力的积累没能克服断层的破裂强度，破裂在此应力状态下无法稳定地向外传播.根据滑动弱化摩擦准则，当剪应力超过断层的破裂强度时，断层两侧才开始相对移动，由此可知，初始应力是决定断层是否发生错动的关键因素之一.

图 4

Fig. 4

图 4 不同初始应力状态下的滑动速率快照图(a-c)和最终滑动量分布图(d-f) Fig. 4 Snapshots of slip rate (a-c) and final slip distribution (d-f) under different initial stress states

相比较而言，在A应力状态下破裂传播较慢，沿走向方向约8 s到达障碍体产生停止震相，最大滑动速率约0.36 m·s^-1，沿断层倾向方向破裂未能传播到断层顶部，最终滑动量分布范围比较小，最大滑动量约为1.1 m；而在B应力状态下，破裂快速向外传播，沿走向方向约6 s到达障碍体产生停止震相，最大滑动速率约0.68 m·s^-1，沿断层倾向方向破裂传播接近断层顶部，最终滑动量分布范围较大，最大滑动量约为1.37 m.这意味着，初始应力设置越接近破裂线，破裂传播越快，最大滑动速率越大；同时，破裂沿着倾向方向的“爬坡能力”更强，最终滑动量分布范围越大，最大滑动量数值越大.

根据Hao等(2013)对芦山地震运动学反演结果，破裂持续时间约10 s左右，最大滑动量约1.2 m，浅于8 km的断层面上滑动量分布接近于0，破裂并未表现明显的“爬坡能力”，因此在当前参数设置下，A状态下的应力取值较为合适.

临界滑动弱化位移D_c是摩擦准则中的重要参数.由于本研究采用分段函数表示D_c随深度变化分布，如图 2b，因此设置不同的d₀即可代表D_c对破裂情况的影响，本文选取d₀=0.3 m、d₀=0.2 m和d₀=0.4 m三组情况值进行破裂模拟.

图 5显示了d₀=0.3 m和d₀=0.2 m的破裂滑动速率快照图(a, b)和最终滑动量分布图(c, d).结果显示，当d₀=0.3 m时，破裂传播较慢，最大滑动速率约0.36 m·s^-1，破裂前锋比较宽，破裂传播稳定，沿断层倾向方向破裂未能传播到断层顶部，最终滑动量分布范围比较小，最大滑动量约为1.1 m；当d₀=0.2 m时，破裂传播较快，最大滑动速率约1.22 m·s^-1，破裂前锋非常窄，沿断层倾向方向破裂传播接近断层顶部，最终滑动量分布范围较大，最大滑动量约为1.45 m.然而，在当前动力学参数设置下，当d₀=0.4 m时，破裂在成核区孕育失败，无法向外传播.这意味着，D_c对破裂滑动速率有着很大的影响.D_c越小，破裂越容易发生，破裂传播越快，最大滑动速率越大；同时，破裂沿着倾向方向的“爬坡能力”更强，最终滑动量分布范围越大，最大滑动量数值越大；但是相对的破裂前锋比较窄，在数值模拟中容易积累误差.

图 5

Fig. 5

图 5 不同d0的滑动速率快照图(a, b)和最终滑动量分布图(c, d) Fig. 5 Snapshots of slip rate (a, b) and final slip distribution (c, d) under different d0

当给定断层面应力情况后，在运动学反演最大滑动量结果的约束下，D_c就只能在较窄的范围内取值.根据郑绪君等(2018)反演得到的芦山地震滑动速率最大为0.48 m·s^-1可知，芦山地震破裂滑动速率并未超过0.5 m·s^-1，因此在当前参数设置下，取d₀=0.3 m较为合适.

由于破裂主要发生在震源附近20 km范围内，因此在设置成核区的时候不宜太大.为了研究成核区的大小对破裂情况的影响，本文选取成核区半径R=3 km、R=3.5 km和R=2.5 km三组情况值进行破裂模拟.

图 6显示了R=3 km和R=3.5 km的破裂滑动速率快照图(a, b)和最终滑动量分布图(c, d).结果显示，当R=3 km时，破裂传播较慢，在6 s时还未到达障碍体，最大滑动速率约0.36 m·s^-1，最大滑动量约为1.1 m；当R=3.5 km时，破裂传播较快，最大滑动速率约0.5 m·s^-1，最大滑动量约为1.45 m；两者的最终滑动量分布范围相差不大.然而，在当前动力学参数设置和网格条件下，R=2.5 km时破裂在成核区孕育失败，无法向外传播.这意味着，在固定的网格划分情况下，成核区需要足够大，蕴含足够多的能量，才能够使破裂成功向外传播；同时，成核区半径大小更多影响的是破裂成核的快慢.成核区半径越大，应力积累的越快，裂纹向外传播的越早，成核区最终滑动量越大.

图 6

Fig. 6

图 6 不同成核区半径的滑动速率快照图(a, b)和最终滑动量分布图(c, d) Fig. 6 Snapshots of slip rate (a, b) and final slip distribution (c, d) under different radius of nucleation zone

根据Zhao等(2013)对芦山地震运动学反演结果，芦山地震在前2 s并未见明显破裂，直到4 s左右破裂才明显向外传播，因此在当前参数设置下，取R=3 km较为合适.

为了使破裂能被触发，在设置成核区初始剪应力T_i时，通常使其略高于断层的破裂强度T_u.为了研究成核区的初始剪应力大小对破裂情况的影响，本文选取T_i=1.01T_u、T_i=1.001T_u两组情况值进行破裂模拟.

图 7显示了不同成核区初始剪应力的破裂滑动速率快照图(a, b)和最终滑动量分布图(c, d)，结果显示，当T_i=1.01T_u时，破裂在3.96 s时已经明显地向外传播，最大滑动速率约为0.36 m·s^-1, 最大滑动量为1.11 m；当T_i=1.001T_u时，3.96 s时破裂仍旧在成核区积累能量，在5.96 s时明显地向外传播，最大滑动速率约0.34 m·s^-1，滑动量分布范围相对较小，最大滑动量为0.96 m.这意味着，设置成核区初始剪应力略大于破裂强度，初步触发破裂，破裂起始后需要一定时间积聚能量才会显著地向外传播.当T_i越大，成核区初始蕴含的能量越大，需要积聚能量向外传播的时间越短，破裂越容易孕育并向外传播，最终滑动量分布范围就越大，最大滑动量值越大；T_i越小，并不会特别影响破裂的滑动速率.

图 7

Fig. 7

图 7 不同成核区初始剪应力的滑动速率快照图(a, b)和最终滑动量分布图(c, d) Fig. 7 Snapshots of slip rate (a, b) and final slip distribution (c, d) under different initial shear stress in nucleation zone

根据Hao等(2013)、Zhang等(2014)对芦山地震运动学反演结果，破裂持续时间约10 s左右，因此在当前参数设置下，T_i=1.01T_u的取值较为合适.

当断层局部地区存在破裂强度不均匀现象，可能导致破裂传播情况发生变化.本文假设震源东北方向有一破裂强度为20 MPa的局部障碍体进行破裂模拟.

图 8显示了加入局部不均匀破裂强度的破裂滑动速率快照图(a)和最终滑动量分布图(b).结果表明，在未传播到局部障碍体的时候，破裂传播一切正常；在5 s时破裂接触到局部障碍体后，由于其破裂强度过大，破裂无法突破，从而产生停止震相，最终破裂没能继续向震源东北方向传播；最终滑动量分布范围不再呈椭圆形分布，震源东北方向滑动量很小.

图 8

Fig. 8

图 8 加入局部不均匀破裂强度的滑动速率快照图(a)和最终滑动量分布图(b) Fig. 8 Snapshot of slip rate (a) and final slip distribution (b) with introducing local inhomogeneous fracture strength

Hao等(2013)、Zhang等(2014)和郑绪君等(2018)对芦山地震运动学反演结果表明，震源东北方向滑动量很小，也就是说破裂在这个区域基本不发生传播，断层最终滑动量分布呈不规则形状，这与图 8结果一致，因此可能是断层面上存在局部障碍体导致的.当破裂传播遇到无法突破的障碍体，就会导致破裂行为发生改变，最终影响断层面上滑动量分布情况.

本文基于芦山地震野外地质调查、余震定位、深地震反射剖面等信息，构建芦山地震铲型断层模型，在震源运动学反演结果的约束下，利用BIEM计算芦山地震动力学破裂过程，讨论不同动力学计算参数对破裂过程和同震滑动分布的影响.结果显示，初始应力是决定断层是否发生错动的关键，初始应力设置越接近破裂线，破裂越容易孕育传播，滑动速率越大，滑动量越大；D_c对破裂滑动速率有着很大的影响，D_c越小，破裂越容易发生，滑动速率越大，但是破裂前锋比较窄，在数值模拟中容易积累误差；成核区半径和成核区初始应力主要影响破裂成核的快慢，成核区半径越大，应力积累的越快，裂纹向外传播的越早，成核区最终滑动量越大；成核区初始应力越大，成核区初始蕴含的能量越大，需要积聚能量向外传播的时间越短，破裂越容易孕育并向外传播，最终滑动量分布范围就越大，最大滑动量值越大；局部不均匀破裂强度主要影响破裂行为和断层最终滑动量分布，当破裂传播遇到局部障碍体，会导致破裂行为发生改变，导致断层最终滑动量分布呈不规则形状.

本文的动力学模拟结果可以很好地再现芦山地震破裂特征，从震源动力学角度解释为何断层面上滑动量分布不规则，对于了解导致这次地震的力学原因有一定的参考价值.但是，由于实际地震的震源过程是多个因素共同作用的结果，在没有其余参数的确切取值的情况下，很难给出某一个物理参数的确切值.这还有赖于未来的综合多种手段的研究来进一步约束.