我们知道，断层两侧介质运动的形态通常是不同的(即：运动速度大小及方向是不同的)，这样经过漫长的地质年代，势必导致断层两侧介质物理性质的不同(bimaterial contrast)，具体表现为断层两侧介质的地震波速度不同(Le Pichon et al., 2005; Lewis et al., 2007; McGuire and Ben-Zion, 2005; Shapiro et al., 2005; Thurber et al., 2006; Wdowinski et al., 2007; Zhao et al., 2008, 2010).前人的大量研究结果表明，断层两侧介质物理性质差异会严重影响断层的自发破裂过程和破裂特征(Adams, 1995; Ampuero and Ben-Zion, 2008; Andrews and Ben-Zion, 1997; Ben-Zion and Andrews, 1998; Brietzke et al., 2009; Duan, 2008; Ranjith and Rice, 2001; Rubin and Ampuero, 2007; Weertman, 1980; Yuan and Zhu, 2016; Zhu et al., 2013; 朱守彪和袁杰, 2016).例如，Yuan and Zhu(2016)利用有限元计算方法考察了断层两侧介质物性差异对地震破裂过程中辐射的地震波在空间的分布特征，发现介质物性差异程度越高，强地面运动加速度分布在空间的不对称性就越强.朱守彪和袁杰(2016)通过计算发现，正是龙门山断裂带两侧的物性差异导致了2008年汶川大地震断层破裂结果的空间不对称性分布(主要呈现朝东北向的单侧破裂).可见断层两侧介质物理性质差异对地震震源过程及其地震灾害空间分布的影响是非常严重的.

一般情况下，断层两侧介质地震波速度的差别在5%~20%的范围内，但也有像美国San Andreas断层那样的个别地区，断层两侧介质的地震波速度差异能够达到30%(Ben-Zion, 2001; Brietzke et al., 2007; Ranjith et al., 2001).我们知道在断层两侧介质的材料属性完全相同的情况下，若断层面上的摩擦本构关系为单一的常系数，且当应力场在空间均匀分布时，则不会出现破裂失稳现象.但是，如果断层两侧介质物性不同，即使摩擦系数保持为常数，也会发生破裂失稳现象(Adams, 1995; Andrews and Ben-Zion, 1997; Ben-Zion, 2001; Weertman, 1980; 朱守彪和张培震，2009).

在断层两侧介质物理性质完全相同的情况下，破裂位错不会导致断层面上的正应力发生变化；然而，当断层两侧介质物性不同时，Ⅱ-型破裂中断层面上的正应力会随着位错发生变化，且其变化大小与位错的空间导数成正比(Adams, 1995; Ben-Zion, 2001; Ranjith and Rice, 2001; Weertman, 1980).这样的力学机制，就会造成在断层较软弱一侧位移方向(通常称为正方向)传播的破裂前端产生动态扩张；而在相反方向(称为负方向)传播的破裂前端出现动态压缩；最后导致断层破裂的不对称性空间分布(Yuan and Zhu, 2016; 朱守彪和袁杰, 2018).此外，破裂的不对称性程度随着破裂传播速度和断层两侧介质物性差异程度的增加而增加(Weertman, 1980; Ben-Zion and Andrews, 1998; Ben-Zion, 2001; Brietzke et al., 2007; Yuan and Zhu, 2016).

但是，前人关于断层两侧物性差异(bimaterial contrast)影响地震破裂过程的研究，都是假定断层两侧介质是均匀、各向同性的，而事实上地球内部介质的物理性质存在着不可忽略的不均匀性特别是各向异性(Cochran et al., 2003; Faccenda et al., 2008; Huang et al., 2015, 2017; Lewis et al., 2007; Peng and Ben-Zion, 2004; Zhao et al., 2008; 高原等, 2010; 滕吉文等, 2012; 王琼和高原, 2018; 朱守彪等, 2017).当断层两侧介质为各向异性材料，并且两侧介质的各向异性程度不同时，这时断层产生自发破裂过程中又可能出现哪种不同形式的力学表现？或者说，断层两侧介质各向异性及其异性程度的不同对断层自发破裂过程有何影响？到目前为止，对这个问题的研究国际上未见有文献报道，但断层两侧各向异性差异对地震破裂过程或许起着重要作用.因此，本文的目的就是要通过有限元计算方法，探讨在断层两侧介质出现各向异性以及各向异性程度存在差异时，断层自发破裂动力学行为的不同特征，以进一步加深我们对震源破裂过程的认识.

为研究断层两侧不同的各向异性材料对震源破裂动力学过程的影响，特别是为抓住主要的力学本质，我们采用二维模型，来研究Ⅱ-型破裂问题.模型中将断层面简化为一条直线.有限元模型几何、初始条件和边界条件见图 1.模型空间是一个100 km×100 km的正方形，内嵌一条50 km长的断层，位于模型的中央.为了便于描述，我们让断层走向沿着直角坐标系X轴的方向(如图 1所示).图 1中，线段AB代表断层的空间位置，红色的五角星表示震源位置(即成核中心).模型四周为吸收边界(ABAQUS中的无限单元)，其目的是为了消除地震波反射造成的影响.此外，我们在整个模型区域内施加统一的法向应力-(张力为正)和剪切应力作为初始应力场.

图 1

Fig. 1

图 1 模型几何形状，初始条件和边界条件.模型区域大小是100 km×100 km；图中线段AB表示断层位置，五角星表示成核区，M1和M2分别表示断层上侧和下侧的介质；模型周围的灰色区域是吸收边界(ABAQUS中的无限单元). Fig. 1 Model geometry, the initial and boundary conditions. The model domain is 100 km by 100 km; Line AB in the figure represents the location of the fault, the star marks the nucleation patch, and M1 and M2 denote the materials in upper and lower parts of the model, respectively; The gray region around the model is the absorbing boundary (infinite elements in ABAQUS).

在模型的有限单元网格剖分中，全部采用三角型网格(见图 1)，目的是消除沙漏现象的影响(Hibbitt et al., 2006).由于我们主要关注断层区域，因此在断层附近的网格密度大，向外网格逐渐变粗.在模型的中央区域，单元边长为100 m；而在模型外围单元边长则为500 m.整个模型有149，295个节点和295，178个单元组成.

需要特别指出的是，由于完全各向异性介质的物性参数其数学表达比较复杂，作为探索研究，文中只研究相对简单的正交各向异性情况.

摩擦本构关系决定着破裂的形式及其演化过程，文中采用前人普遍使用的速率和状态相关的摩擦本构关系(Ampuero and Ben-Zion, 2008; Cochard and Madariaga, 1996; Shaw and Rice, 2000).具体数学表达式如下：

(1)

其中μ_s是静摩擦系数，V是滑移速率，V_c是特征速度，α_f和β_f是两个表示演化效果的常数，θ是状态参量，t_c是特征时间尺度.

通常情况下，当断层两侧介质物性不同时，断层自发破裂会出现Adams不稳定现象(Adams, 1995; Cochard and Rice, 2000; Ranjith and Rice, 2001).这种现象在物理上表现为破裂产生的脉冲随着传播距离变得更加尖锐并且发散.从数学角度看，该系统对微扰的响应是病态的，并且数值计算结果随着网格的变化而变化(Ben-Zion and Huang, 2002).因此，要得到稳定的解，需要进行正则化处理.根据前人数值模拟的结果(Brietzke et al., 2009; Cochard and Rice, 2000; Ranjith and Rice, 2001; Rubin and Ampuero, 2007; 朱守彪和袁杰, 2016)，计算时通过延缓断层面上的剪切应力随正应力变化的过程来实现正则化，从而提高计算的稳定性.实际数值模拟中，断层面上的剪切应力随时间的变化率为(DeDontney et al., 2011)

(2)

其中σ、τ、μ_f分别为正应力、剪切应力以及摩擦系数；t^*为正则化时间尺度，取值为：4Δx/C_s1(其中，Δx为滑移面上的单元尺度，C_s1为断层两侧中较柔软一侧介质的S波速度).理论上t^*应远大于计算的时间步长，通常在实际计算时大一个数量级即可(朱守彪和袁杰, 2016)，具体数值见表 1.

表 1 (Table 1)

表 1 模拟中所采用的共同模型参数 Table 1 Common parameters for simulations in the models 物理量描述 数学符号(单位) 数值 M1的X方向杨氏模量 E1x(GPa) 70 M1的Y方向杨氏模量 E1y(GPa) 70 M2的X方向杨氏模量 E2x(GPa) 90 M2的Y方向杨氏模量 E2y(GPa) 70 泊松比 v 0.25 剪切模量 G12(GPa) 30 初始正应力 σ0(MPa) 200 初始剪应力 τ0(MPa) 69 介质密度 ρ(kg·m-3) 2700 静态摩擦系数 μs 0.6 直接演化系数 αf 0.01 演化系数 βf 0.40 特征速度 Vc(m·s-1) 0.7 特征时间 tc(s) 0.3 速率弱化特征长度 Dc(m) 0.21 成核区大小 Lnucl(km) 1.0 正则化时间尺度 (s) 1.0×10-3

表 1 模拟中所采用的共同模型参数 Table 1 Common parameters for simulations in the models

本研究利用ABAQUS/Explicit商业有限元软件(有限元显式计算程序)来模拟断层自发破裂的动力学过程，该计算方法与朱守彪和袁杰(袁杰和朱守彪, 2014a, 2014b; 朱守彪和袁杰, 2016; Zhu and Zhang, 2010)中相同；计算中时间步长取为0.0001 s，从而保证模拟结果的正确性.

为深入考察断层两侧介质的各向异性分布对破裂动力学过程的影响，文中我们将构建多种情况的模型来一一分析在不同的情况下断层的破裂特征.模型Ⅰ是断层两侧介质都为正交各向异性，并且材料主轴方向分别与坐标轴X轴、Y轴方向平行(见图 2).然后，将模型Ⅰ中的断层下侧材料(M₂)主轴方向分别旋转45°和90°，便得到模型Ⅱ及模型Ⅲ.下面将通过有限元模拟，依次给出这三个典型模型中，断层面上的位错、断层面上的应力变化以及地震波产生的地震动峰值加速度空间分布.

图 2

Fig. 2

图 2 模型Ⅰ中断层两侧材料主轴方位分布 Fig. 2 Spatial distribution of principal orientation of orthotropic materials across the fault in Model Ⅰ

需要指出的是，三个模型中，除了特别说明之外，所有模型参数、初始应力场、摩擦本构关系等都是相同的.此外，三个模型中，破裂都是在断层中心处成核，并且成核方式都与朱守彪等(朱守彪和袁杰, 2016; 朱守彪等, 2017)的方法相同，此处不再赘述.

在模型Ⅰ中，断层上侧及下侧介质皆为正交各向异性材料，并且材料主轴的方向都是一样的，即：上侧M₁中的E_1x与下侧(M₂)中的E_2x是平行的(实际上E_1y也平行于E_2y)(见图 2).但是，模型的上、下侧中，不同方向的杨氏模量大小不等，其中：E_1x =70 GPa, E_1y =70 GPa, E_2x=90 GPa, E_2y=70 GPa, ν=0.25, G₁₂=30 GPa(除特别说明外，下文中的泊松比及剪切模量不变).所以该模型中材料的各向异性表现为下侧介质中水平向(断层走向)的杨氏模量明显高于垂直向的杨氏模量.

有限元模拟中，一旦破裂在断层中央成核，在不需要任何外界因素的影响下，仅仅由于摩擦关系的控制，破裂可以自发地沿着断层向两侧传播.

图 3是有限元计算给出的破裂沿着断层向两侧传播过程中产生的位错空间分布图.图中可见，断层两侧位错分布很不对称.X轴正向(成核中心右侧)的位错量明显小于X轴负方向的位错，X轴负向的最大位错量可达~4 m, 而X轴正向的最大位错量小于2 m.此外，破裂从成核中心向两侧传播的距离也存在明显的不对称性，呈现单侧破裂的情形.此外，图 3还显示，断层面上任一点的位错持续时间很短，远远小于整个断层破裂的持续时间，因此该破裂为脉冲型(pulse-like)破裂.

图 3

Fig. 3

图 3 模型Ⅰ中断层两侧位错空间分布 Fig. 3 Slip distribution along the fault in Model Ⅰ

为探求不对称破裂的力学机制，文中将导出破裂过程中断层面上的正应力、剪应力的数值，其结果显示在图 4中.该图是断层面上不同时刻的正应力、剪应力分布的快照.图 4左侧是断层面上的正应力、右侧为剪应力分布.由图结果可见，正应力、剪应力在不同时刻的空间分布是不同的.通过仔细观察不难发现，正应力、剪应力在成核中心两侧空间分布也是不对称的，特别是正应力在成核中心左侧(X轴负方向)，正应力的绝对大小在缩小，即该方向上的断层面之间的压力在变小，有拉张的趋势，从而有助于断层破裂的发生；但在X轴的正方向，正应力的绝对值在增大，即断层处在不断压紧的状态，这样不利于断层破裂.同时，图 4还显示，断层面上的剪应力成核中心几乎是对称分布的.因此，从不同时刻断层面上的应力状态，容易理解为什么断层破裂从成核区发生后，出现空间上不对称的单侧破裂.

图 4

Fig. 4

图 4 断层面上的正应力及剪切应力在不同时刻的分布(左侧为正应力，右侧为剪应力) Fig. 4 Distributions of normal and shear stresses at different times along the fault. (The normal stress is denoted on the left side and the shear stress on the right side)

断层在发生破裂的同时，除了产生地震位错外，由于断层及其周围介质质点的振动而释放地震波(Yuan et al., 2015; Zhu, 2018).通常离开断层不同距离的质点其地震波的振动方式不同(震相不同).一般来说，离开地震震源越远，地震波衰减越严重，地震动就越小，相应地地震灾害也就越轻.为研究模型Ⅰ中断层破裂造成的地震灾害情况，本研究给出了模型中所有质点振动的峰值加速度(PGA)空间分布云图(图 5).由图可见，峰值加速度在空间的分布具有很强的空间不对称性.在成核中心(震源)附近，峰值加速度的数值可达~10 m·s^-2, 但随着离开断层距离的增大而很快衰减，即：地震灾害在震源附近最强，随着离开震源越远，震害越小.图 5还显示，峰值加速度关于震源呈现明显的不对称性分布，震源右侧，峰值加速度衰减快，而在震源左侧方向，峰值加速度的高值区分布范围大，可见这些地区的地震灾害可能会很严重.由此可见，研究断层的破裂特征，对于分析、估计地震灾害及防震减灾等有重要的科学意义及实际应用价值.

图 5

Fig. 5

图 5 模型Ⅰ中所有质点振动峰值加速度(PGA)的空间分布云图 Fig. 5 Spatial distribution of peak acceleration of strong ground motion (PGA) in Model Ⅰ

将模型Ⅰ中断层下侧的材料主轴旋转45°，断层上侧物质材料主轴方向保持不变，这样便可得到模型Ⅱ，见图 6.

图 6

Fig. 6

图 6 模型Ⅱ中断层两侧材料主轴方位分布 Fig. 6 Spatial distribution of principal orientation of orthotropic materials across the fault in Model Ⅱ

模型Ⅱ中所有模型参数、边界条件、初始应力场及摩擦关系等都与模型Ⅰ中的一致.下面是利用有限元方法对模型Ⅱ进行模拟的计算结果.

图 7是模型Ⅱ中断层两侧位错空间分布.图中可以明显看出，向左侧传播的破裂其位错量比向右侧方向的位错量大，破裂传播的距离相等.对比图 2可见，模型Ⅱ的结果与模型Ⅰ明显不同.由计算可知，这种结果的差异都是由于断层下侧材料主轴的旋转导致的.由此可见，各向异性材料的材料主轴方位对震源破裂过程影响很大.

图 7

Fig. 7

图 7 模型Ⅱ中断层两侧位错空间分布 Fig. 7 Slip distribution along the fault in Model Ⅱ

图 8表示模型Ⅱ中断层面上的正应力及剪切应力在不同时刻的分布图.图中显示，成核中心左侧方向上的正应力在不断变小，而剪应力变化不大，从而导致了左侧方向的位错大于右侧方向的位错.图 9是模型中各点峰值加速度分布云图.由图可见，峰值加速几乎是空间对称分布的，这与图 5中模型Ⅰ的结果明显不同，所以材料各向异性的方位对于地震破裂及其震害的影响是非常重要的.

图 8

Fig. 8

图 8 模型Ⅱ中断层面上的正应力及剪切应力在不同时刻的分布(左侧为正应力，右侧为剪切应力) Fig. 8 Distributions of Normal and shear stresses at different times along the fault (The normal stress on the left side and the shear stress on the right side)

图 9

Fig. 9

图 9 模型Ⅱ中所有质点振动峰值加速度(PGA)的空间分布图 Fig. 9 Spatial distribution of peak acceleration of strong ground motion (PGA) in Model Ⅱ

让模型Ⅱ中断层上侧物质其材料主轴方向保持不变，将断层下侧的材料主轴再旋转45°，便得到模型Ⅲ，见图 10所示.同样，模型Ⅲ中所有模型参数、边界条件、初始应力场及摩擦关系等参数都与模型Ⅰ中的一致.下面给出利用模型Ⅲ进行的计算结果.

图 10

Fig. 10

图 10 模型Ⅲ中断层两侧材料主轴方位分布 Fig. 10 Spatial distribution of principal orientation of orthotropic materials across the fault in Model Ⅲ

图 11是模型Ⅲ中断层两侧位错空间分布.图中结果可见，向左侧传播的破裂其位错量比向右侧方向的位错量大，破裂传播的距离相等(与模型Ⅱ相同).进一步对比图 2可见，模型Ⅲ的断层位错结果与模型Ⅰ的明显不同.图 12是模型中各点峰值加速度分布云图.该图中峰值加速的空间分布几乎是对称的，这与图 5中模型Ⅰ的结果也是区别很大.

图 11

Fig. 11

图 11 模型Ⅲ中断层两侧位错空间分布 Fig. 11 Slip distribution along the fault in Model Ⅲ

图 12

Fig. 12

图 12 模型Ⅲ中所有质点振动峰值加速度(PGA)的空间分布图 Fig. 12 Spatial distribution of peak acceleration of strong ground motion (PGA) in Model Ⅲ

实际上，模型Ⅲ中断层下侧的材质主轴的两个正交方向正好与模型Ⅰ的进行了交换，本来沿着断层方向的材料主轴变换成为垂直于断层方向.所以不难看出，这种结果的差异说明各向异性材质沿着断层方向的材料属性比垂直于断层方向的材料属性对于地震破裂过程的影响更为显著.

上文的模型里，断层两侧物质的各向异性性质存在差异，但在实际地球介质中断层两侧出现各向异性，有可能断层两侧的异性材料是相同物理性质的(相同的各向异性).所以，本文还将对断层两侧相同的各向异性材料情况进行研究，为此构建另一模型：模型Ⅳ.同样，在模型Ⅳ中，所有模型参数、断层几何、初始应力场、摩擦关系等均与模型Ⅰ中的相同，只是在模型Ⅳ中，E_1x=80 GPa, E_1y =60 GPa, E_2x=80 GPa, E_2y=60 GPa；即：断层两侧材料均为各向异性，并且它们各向异性的性质是相同的.

图 13a是计算给出的模型Ⅳ中断层两侧位错空间分布.图中显示，断层两侧的位错关于成核中心(震源)呈现很好的对称分布.

图 13

Fig. 13

图 13 有限元计算给出的不同模型断层面上的最终位错分布图 (a)表示模型Ⅳ中断层两侧位错空间分布；(b)是模型Ⅴ中断层两侧位错空间分布；(c)为模型Ⅵ中断层两侧位错空间分布 Fig. 13 Final slips distributed along the fault in different models (a) Slip distribution along the fault in Model Ⅳ; (b) Slip distribution along the fault in Model Ⅴ; (c) Slip distribution along the fault in Model Ⅵ

按照上文同样的方法，在模型Ⅳ之上，让断层下方的介质其材料主轴旋转45°，得到模型Ⅴ.在模型Ⅴ中，所有模型参数皆与模型Ⅳ相同.通过有限元计算，可以获得模型Ⅴ中断层破裂过程的详细信息.图 13b是模型Ⅴ中断层破裂位错的空间分布图.由图可见，位错在空间的分布也是近似对称的.但与图 13a相比，可以看出，由于模型Ⅴ中断层下侧材料主轴方向的旋转，位错分布的空间对称性程度不如图 13a的结果，同时其位错幅度比图 13a的要大.由此说明，各向异性材料主轴的旋转对断层破裂还有一定的影响，但其影响程度不是很显著.

此外，在模型Ⅳ之上，让断层两侧介质各向异性材料的剪切模量增大，即从原模型Ⅳ的G₁₂=30 GPa变化为38 GPa，作为模型Ⅵ.在其他所有模型参数都保持与模型Ⅳ中的一样时，通过有限元方法计算的断层最终位错结果如图 13c所示.显然，图 13c中破裂向两侧对称传播，但破裂很快就被终止，这种破裂所形成的地震震级会远远低于模型Ⅳ中的震级，这样造成的地震灾害也会大大减小.通过考察多个模型发现，随着正交各向异性材料参数的变化，断层破裂的总体特征也会发生很多复杂的变化形状，模型Ⅵ的结果就是其中典型的代表.

最后，研究中还对摩擦关系对破裂结果的控制，网格密度大小对计算结果稳定性的影响，计算时间步长等影响进行了试验计算，发现文中的结果是稳定、可靠的.

另外需要说明的是，本研究中模型Ⅰ里断层出现单侧破裂，而实际的天然大地震中又有很多是单侧破裂的实例(如2001年发生的中国昆仑山口西M_S8.1级特大地震，2008年四川汶川M_S8.0大地震，等).那么，这些大地震的单侧破裂是否是由于断层两侧介质各向异性差异很大造成的？现在还无法给予回答，需要今后进行深入系统的研究.虽然，目前地震学方法对震源区介质的各向异性的分辨率还较低，还远没有达到可以利用其结果进行精细模拟断层破裂过程的阶段，但随着地震科学的发展，更为精细的材料模型结果会给震源过程的模拟带来新的科学发现.

另外，文中是二维模型，还有很多完全各向异性材料的弹性常数与破裂之间的复杂关系仍需继续深入研究.尽管如此，本研究作为尝试，首次对断层两侧材料各向异性对地震破裂特征的影响进行了定量分析，研究结果对深入认识震源动力学过程有重要的科学意义.

根据文中的计算与分析，得出如下初步结论：

(1) 利用有限单元计算方法，可以在断层两侧材料存在各向异性性质时模拟断层的自发破裂动力学过程.尽管断层两侧介质为各向异性材料时的破裂问题是个高度的非线性力学问题.通过加密网格、缩小计算时间步长、摩擦关系的正则化等处理，发现模拟结果是收敛、稳定、可靠的.

(2) 断层两侧介质物性的各向异性对断层破裂动力学过程有重要影响，并且材料参数与破裂行为之间的关系非常复杂.对于正交各向异性材料，沿着断层走向的材料主轴上的杨氏模量对于破裂过程的影响大于垂直于断层走向的杨氏模量，特别是当平行于断层走向的材料其一侧杨氏模量显著大于另一侧时，断层出现不对称的双侧破裂(或称为单侧破裂)；但随着材料主轴方位的变化，这种破裂的空间不对称性分布会变化，说明各向异性材料的主轴方位对破裂特征也产生影响.

(3) 若断层两侧材料仅仅为相同的各向异性时(即断层两侧为同样的各向异性材料)，则不会影响断层破裂的空间对称性分布，这种情况与断层两侧介质均匀各向同性材料时的效果有类似之处.当其中一侧的各向异性材质主轴方位发生变化时，断层破裂的空间对称性会受到一定程度的影响，但其影响很小；然而，随着正交各向异性材料剪切模量的增加，断层破裂容易被终止，这样情况下就不易产生特大地震.

可以期待，随着模拟技术的改进、计算设备的升级以及断层附近介质物性分辨率的提高，符合实际地震震源过程的模拟对于地震灾害的评估及地震物理预测等会发挥越来越重要的作用.