2008年发生的汶川大地震,断层破裂过程非常发杂,地表所能观测到的破裂长度达300余千米.地震发生后,国内外很多学者陆续利用多种方法和不同资料(如：地震波形、地表破裂、强地面运动、GPS和InSAR资料等)对汶川地震的破裂过程的运动学模型进行了研究(Ji and Hayes, 2008；王卫民等,2008；张勇等,2008；杜海林等,2009；Shen et al., 2009；Xu et al., 2009；赵翠萍等,2009； Nakamura et al., 2010；Zhang and Ge, 2010；Wang et al., 2011；Fielding et al., 2013; Pang et al., 2014; Zhang et al., 2014; Zhang and Wang, 2015),取得了多种反演结果.尽管不同作者给出的结果不尽相同,但描述的汶川地震破裂的运动学特征大体如下：

破裂首先在位于映秀—北川断裂带上的震源处起始,在最先的开始阶段,表现为双侧的破裂形式,在东北及西南的两个方向同时产生破裂,但向西南方向的破裂很快就终止.而在断层的东北方向,破裂继续前行,跨越了多个断层阶区,破裂经过的断层面形状复杂(如弯曲、拐折、交叉等),破裂最后在青川断裂带上终止.汶川大地震破裂的总长度大约为300 km,破裂持续时间大约100 s; 破裂速度在不同的地方差异很大(但基本上都小于剪切波速度),静态位错分布也很不均匀,断层的倾角由南向北越来越大,位错方向由近纯逆冲逐渐过渡到走滑类型.

尽管运动学模型告诉我们断层破裂的长度、破裂速度、位错大小等很多物理量,也为认识汶川地震及抗震救灾发挥了重要作用.但是,是什么因素造成了这样的破裂格局？或者说,汶川地震破裂过程中,什么物理量起着决定性作用？这些问题至今依然没有给出回答.实际上,这些问题在物理上是动力学问题,只有通过对汶川地震的破裂过程进行动力学分析,才有可能加以解决.汶川地震发生后,其震源破裂的动力学过程倍受人们的关注,国内外不少专家对此进行了探索,获得了一些认识.

Duan(2010)利用三维有限元方法对汶川地震的破裂动力学过程进行了研究.模型中,汶川地震的破裂面假定为一条平直断层面,其倾角为33°；介质材料均匀,初始应力场不均匀并随深度变化.通过研究发现,随着初始应力场方向的改变,断层的破裂速度会发生很大的变化,初始应力场的方向还可以导致截然不同的位错格局及地表运动形态.特别是当初始场的指向与断层走向之间的夹角处在最不利于断层破裂的方位时,破裂会自然地停止.汶川地震在断层东北端的终止以及两处破坏最强的区域,也许就是初始场的方向在起着作用.Wen等(2012)利用显示的三维有限元方法,所用模型与Duan(2010)的模型基本相同,设定了不均匀的初始应力条件,模拟结果表明初始应力场决定了破裂速度、位错的空间分布等主要破裂特征.

Zhu和(2010,Zhang2013)利用改进的有限单元模型,分两步(即先反演,后正演)分别模拟了汶川地震震间变形及同震破裂过程.第一步利用大地测量资料(水准测量及GPS观测)作约束,反演震间的最佳数值模型,获得汶川地震前的应力、应变状态；然后据此最佳模型,进行正演计算模拟断层的自发破裂过程.计算结果显示：下部缓倾角断层上的破裂对上部陡倾角上的错动有触发作用,破裂过程中,断层面上的剪切应力和法向应力都有利于断层的错动,这可能也是高倾角逆冲断层上发生大地震的原因之一.

虽然上述学者对汶川地震破裂动力学机制研究方面作了不少探索,取得一定的进展,但汶川地震破裂的最基本特征为单侧破裂(或者认为是不对称的双侧破裂),对这个问题至今未见有文献报道.此外,我们注意到龙门山断裂带两侧介质物理性质差异很大(Lei and Zhang, 2009; Zhu and Zhang, 2013; Liu et al., 2014),而断层两侧介质物理性质差异(bimaterial contrast)会严重影响断层的破裂过程和破裂特征(Weertman,1980; Adam,1995; Andrews and Ben-Zion, 1997; Ben-Zion,2001; Ranjith and Rice, 2001; Ben-Zion and Shi, 2005; Ben-Zion et al., 2007; Rubin and Rice, 2007; Ampuero and Ben-Zion, 2008; Duan,2008; Ma and Beroza, 2008; Olsen-Kettle et al., 2008; Brietzke et al., 2009; DeDontney et al., 2011; Langer et al., 2012; Scholz,2014; Yuan and Zhu, 2016).所以本文将利用有限元方法,考虑龙门山断裂带两侧介质物性的差异,模拟汶川地震的破裂过程,试图定量给出汶川地震为什么会出现单侧破裂的动力学机制.

为研究汶川地震单侧破裂的物理机制,首先必须对实际复杂问题进行抽象与简化.

汶川地震中,断层破裂规模巨大,断层几何形状复杂多变,既有平面也有曲面,断层面倾角自南向北不断发生变化,断层在空间上还存在不连续分布等多种特征.此外,地震中断层的位错主要表现为由逆冲向走滑性质过渡：从震源处的几乎纯逆冲,然后不断地由南向北过渡到几乎纯走滑类型.研究中,为抓住主要矛盾,重点突出汶川地震单侧破裂过程的物理本质,研究中将实际的三维地质体简化为二维模型(通过震源的水平面),将复杂的断层几何简化为平面,复杂的断层破裂过程简化为简单的走滑类型.模拟中,同时还假定汶川大地震之前龙门山断裂带及其附近区域处在均匀的应力环境里,这样文中研究的破裂问题就简化为标准的滑开型(Ⅱ-型)破裂模式.计算时,首先让破裂在汶川地震的震源处成核,然后让破裂自由地向东南—西北两个方向自发地传播、演化.

对于二维模型,当模型四周利用无限单元来吸收地震波时,破裂动力学过程的能量表达形式为(Hibbitt et al., 2006)

(1)

其中t为时间,ρ为介质密度,v为速度场矢量,U为内能,f是体力,T^l为表面分布载荷,T^qb为无限单元的辐射拉力,T^f为摩擦力,E·_W为外力做功速率,E·_QB为无限单元阻尼衰减的能量速率,_F为摩擦接触的能量衰减速率.

此外假定断层面上的法向位移为连续,即

(2)

断层面上的质点所处的运动状态表达式为

(3)

其中σ_n法向接触应力,τ为剪切应力,μ为摩擦系数.

利用有限元方法对于上述(1)—(3)的方程进行求解,即可得到模型中位移、速度、加速度,应力,应变以及能量等物理量.

为保证计算精度及计算效率,本研究计算采用显示的直接积分方法.尽管基于条件稳定的中心差分方法不存在稳定性问题,但时间步长必须小于临界值,其大小为：(L为模型中最小网格边长).根据本文模型中的模型参数(表 1),该临界值为Δt_cr≈0.019 s；实际计算时,取时间步长为1.0×10^-4s. 显然,该值远远小于临界值,从而保证计算过程的稳定性与计算精度的要求.

表 1(Table 1)

表 1 模拟中所用的模型参数 Table 1 Model parameters used in simulation 参量物理性质 参数(单位) 数值 四川盆地S波速度 Vs1(m·s-1) 4120 四川盆地P波速度 Vp1(m·s-1) 7150 川西高原S波速度 Vs2(m·s-1) 3380 川西高原P波速度 Vp2(m·s-1) 5850 初始剪应力 τ0(MPa) 69 初始正应力 σ0(MPa) 200 介质密度 ρ(kg·m-3) 2900 静摩擦系数 μ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 正则化时间尺度 t*(s) 1.2×10-3

表 1 模拟中所用的模型参数 Table 1 Model parameters used in simulation

根据汶川地震的震中位置以及龙门山断裂带的实际构造,研究中构造了二维有限元模型,如图 1所示.模型空间几何尺度为500 km×500 km.为减少地震波辐射对破裂行为的影响,模型四周施加了吸收边界(利用ABAQUS中的无限单元来实施).地震破裂的成核位置位于模型的中心(图中粗黑线段).有限元网格全部为三角形单元,单元最小边长为100 m,单元最大边长为500 m.破裂带长度为500 km,利用ABAQUS中接触单元来表示.由于初始应力场影响破裂行为,所以模拟中根据对龙门山断裂带及附近地区地应力的观测结果(Wu et al., 2009)来选择初始应力场.

图 1

Fig. 1

图 1 有限元模型图中坐标原点处的粗实线为汶川地震的震中(成核区). 龙门山断裂带的走向沿X轴正方向(东北方向). Fig. 1 Sketch of the finite element model The thicksolid line at the origin of coordinates stands for the epicenter of the 2008 Wenchuan earthquake. The strike of the Longmen Shan fault is in the direction of X-axis,namely northeast.

摩擦本构关系决定着破裂的形式及其演化特征,根据前人对断层两侧物性不同发生破裂情况的模拟经验(Rubin and Ampuero, 2007; Ampuero and Ben-Zion, 2008; Ma and Beroza, 2008; Finzi and Langer, 2012; Gabriel et al., 2012；Langer et al., 2012),研究中选择速度-状态相依的摩擦本构关系,具体数学表达式为

(4)

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

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

(5)

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

在数值模拟中,要使破裂过程启动,通常利用减小摩擦系数来实现.具体计算中是将成核区内的静态摩擦系数突然减小,表达式为

(6)

破裂在成核区一旦形成,在其外的区域里,破裂行为由摩擦本构关系及应力环境来决定.破裂可以自发地向断层两侧传播、演化,直至最后终止.

为模拟2008年汶川地震的破裂过程,我们根据前人对龙门山地区的地震层析成像、以及其他地球物理反演结果等资料,建立具体的有限元模型.模型中龙门山断裂带及附近地区介质的地震波速度根据地震层析成像结果(以汶川地震震源深度处)来选取(Lei and Zhao, 2009; Liu et al., 2014),数值见表 1.这样龙门山断裂带两侧(川西高原与四川盆地)介质的地震波速度相差大约为20%.模拟中关于摩擦行为的基本参数参考前人的模拟结果来选定(Rubin and Ampuero, 2007; Ampuero and Ben-Zion, 2008; Ma and Beroza, 2008; Finzi and Langer, 2012).

图 2为模拟给出的汶川地震震中(成核中心)两侧断层面上的滑移量分布(为清楚起见,该图只展示了局部区域,以下同).由图可知,破裂为明显的不对称分布,向北东方向的破裂距离明显大于向西南方向的距离,位错量也是在东北方向的大于向西南向的.因此许多学者认为汶川地震为单侧破裂类型(实际为不对称的双侧破裂).破裂向东北方向传播距离远,而向西南方向传播的距离很小,仅10余千米.

图 2

Fig. 2

图 2 震中两侧断层面上的位错在不同时刻的空间分布 Fig. 2 Spatial distribution of fault displacements on both sides of the epicenter at different times

此外,图 2还显示,断层面上任一点的位错持续时间很短,远远小于整个断层破裂的持续时间,因此该破裂为脉冲型(pulse-like)破裂(袁杰和朱守彪,2014a).另外,我们还可以从图 2中大体知道破裂的传播速度为2.7 km·s^-1,因此属于亚剪切破裂(破裂速度小于剪切波速度).该破裂速度与张勇等(2008)以及Xu等(2009)通过反演给出的运动学震源破裂模型的破裂速度大体相当.

根据图 1所示的应力边界条件以及川西高原所在一侧是物质相对较柔软的格局,我们知道本研究中破裂的正方向应该是朝龙门山断裂带的东北方向,即图中的x轴正向；那么破裂的负方向就是x轴的反向(Andrews and Ben-Zion, 1997; Ben-Zion and Andrews, 1998; Ben-Zion,2001; Ranjith and Rice, 2001; Ben-Zion and Shi, 2005; Shi and Ben-Zion, 2006; Ben-Zion et al., 2007; Rubin and Ampuero, 2007; Ampuero and Ben-Zion, 2008).对于双层不同性质的断层,破裂的正向应为破裂的优势取向,即破裂在东北方向应有较好的破裂效果,图 2的结果正是说明了这一点.

图 3为滑移速率在断层上的不同位置的分布.由图可见,在破裂的正方向上(东北向),速度脉冲的幅值随着破裂传播距离的增大而越来越大,说明破裂具有自己放大的能力,这样地震破裂就能持续很长距离而不会衰减；此外,位错速率脉冲的宽度(最大振幅一半处的脉冲宽度)在正方向上是越来越小、越来越尖锐,说明破裂能量向低波长方向转移.但在破裂的负方向,脉冲变小,最后终止.图 3还显示,每个脉冲起始后,在很短的时间里就达到最大值,然后迅速降低为0(意味着断层很快就愈合).可见,这种脉冲型破裂具有自己愈合的能力,这正是脉冲型破裂(pulse-like)的典型特征.图 3所示的破裂脉冲就是通常所称谓的Weertman脉冲,它能自己愈合、自行维持、自行放大,对地震破裂过程的影响非常深远(Weertman,1980; Andrews and Ben-Zion, 1997; Ben-Zion,2001; Shi and Ben-Zion, 2006; Ma and Beroza, 2008; Ampuero and Ben-Zion, 2008; Scholz,2014; Brietzke et al., 2009; DeDontney et al., 2011; Finzi and Langer, 2012; Langer et al., 2012; Scholz,2014; Yuan and Zhu, 2016),一直受到国际地震学界的广泛重视.由于破裂过程中产生的Weertman 脉冲具有自己放大、自行维持的属性,所以这种破裂不易于停止,这就是为什么汶川地震的震级是如此巨大的重要因素.

图 3

Fig. 3

图 3 断层面上滑移速率沿着断层上不同位置时的分布(时间间隔为1.5 s) Fig. 3 Distribution of slip rates at different locations along the fault at regular time interval(1.5 s)

图 4给出了断层面上的正应力及剪切应力在不同时刻的分布.图中清楚地显示,断层上的正应力是随着时间变化的.但在模拟时模型施加的初始应力场为均匀场,即断层面上的正应力应该是处处相等,不会发生变化.图 4结果表明,随着破裂传播距离的增大,张性正应力变化会越大越大,特别是当断层很长,破裂传播的距离很远时,可以让断层处在拉张的应力状态,这样破裂就会变得更加容易,摩擦产生的热几乎为0,断层的双侧介质效应会愈加明显.所以震级俞大,Weertman脉冲的影响就会愈加显著(Ranjith and Rice, 2001; Ben-Zion and Huang, 2002).

图 4

Fig. 4

图 4 断层面上的正应力及剪切应力在不同时刻的分布 (a)正应力；(b)剪切应力. Fig. 4 Distribution of normal and shear stresses on the surface of the fault at different times (a)Normal stress;(b)Shear stress.

为了搞清楚断层面上的正应力变化的机理,我们又构建了另外一个模型.该模型中断层两侧介质的物理性质完全一样,即整个模型介质均匀(其他模型参量均保持不变).这时,模拟给出的破裂传播过程中断层面上的应力分布如图 5所示. 通过与图 4的比较可以发现,当模型中介质均匀时,破裂传播过程中断层面上的正应力不随时间及位置的变化而变化.

图 5

Fig. 5

图 5 当模型中介质性质均匀时,破裂传播过程中断层面上的正应力及剪切应力在不同时刻的分布 (a)正应力；(b)剪切应力. Fig. 5 Distribution of normal and shear stresses on the surface of the fault at different times in rupture propagation when the material in the model is not uniform (a)Normal stress;(b)Shear stress.

实际上,早在1980年,对于断层两侧介质性质不同的破裂传播问题,Weertman(1980)就通过数学解析推导,给出了断层面上的剪切应力变化及正应力变化的数学表达式为

(7)

(8)

其中δ为断层面上的位错,和μ^*分别是弹性模量及破裂速度的函数,对于均匀介质,μ^*= 0,₀=μ/(1－ν),这里μ是剪切模量,ν为泊松比,下标0表示准静态变形过程.若模型中断层两侧介质物性有差异,在广义瑞利波速度有定义时,即：在断层两侧介质的S波速度相差小于35%的范围(Ben-Zion and Andrews, 1998)的情况下,则破裂在传播过程中断层面上的正应力会发生变化,即会产生Weertman脉冲.此外,根据图 4及图 5可知,我们模拟所展示的断层面上的剪应力及正应力变化的结果与Weertman的数学解析解是完全吻合的,这也从一个方面说明文中模拟结果的可靠性.

Andrews和Ben-Zion(1997)将介质在断层物性较软弱一侧的运动形象地比作为地毯在地板上的移动.当地毯在地面上移动时则会出现褶皱(wrinkle),介质瞬时从接触面离开,因此在褶皱传播方向上就产生了一个净的剪切位移,导致正应力变小(Andrews and Ben-Zion, 1997; Ben-Zion and Andrews, 1998; Cochard and Rice, 2000; Ranjith and Rice, 2001; Ben-Zion and Huang, 2002; Ben-Zion and Shi, 2005).

事实上,正是由于龙门山断裂带两侧介质物性的差异,导致了破裂在东北方向的传播过程中断层面上的正应力减小,特别是当传播距离越来越远时,甚至会出现完全的拉张正应力状态(Rubin and Ampuero, 2007),这非常有利于破裂的传播(破裂的正方向),使得破裂能够自己维持、自己愈合、自行放大,破裂过程中摩擦产生的热能极小,导致破裂传播的距离很远,最终导致地震的震级很大(Brietzke et al., 2009).

本文运用完全动力学方法模拟了断层的破裂过程,显示了龙门山断裂带两侧介质物理性质差异是导致汶川地震单侧破裂的决定性因素.为研究介质物性差异对地震破裂过程的影响,下面再展示一下上文提到的将模型中断层两侧材料性质视为均匀的模型中位错在空间的分布结果.图 6是当模型中介质均匀时,位错在空间的分布图.图中显示,当断层两侧介质性质完全相同时(其他模型参量不变),破裂确实是向两侧传播,并且在空间上是对称分布.但传播的距离非常之短,总共只有16 km,因此在这种情况下,产生不了像汶川地震(M_s8.0)那样的特大地震.由此可见,断层两侧物性的差异对于破裂过程的影响是多么巨大.通过图 6与图 2的比较,我们清楚地认识到龙门山地震带两侧介质物理性质的差异导致了汶川大地震的单侧传播过程(或不对称的双侧破裂).

图 6

Fig. 6

图 6 当断层两侧介质均匀时,位错随在不同时刻的空间分布图 Fig. 6 Spatial distribution of displacements on both sides of the fault from at different times when the material on both sides of the fault is uniform

在此之上,文中还专门针对断层两侧介质物性差异程度对破裂传播过程的影响进行了模拟分析.当断层两侧介质物性差异分别为5%、10%以及15%时,断层破裂都会出现不对称的双侧破裂行为,但不对称性随着物性差异程度的增大而不断增大(Yuan and Zhu, 2016).

龙门山断裂带及周边地区的应力状态应该是不均匀的.但Duan(2010)以及Wen等(2012)的研究结果表明,初始应力场的不均匀性会严重影响破裂的静态位错、应力降以及破裂速度的分布,所以文中所用的均匀初始场不会影响对汶川地震单侧破裂传播过程这一物理本质的认识.此外,文中的简单模型导致断层破裂很难终止,之所以汶川地震最后终止于青川断裂带,Duan(2010)认为是初始应力场的方位改变造成了破裂的终止.所以初始应力场的非均匀性,不影响我们对汶川地震破裂过程不对称的数值模拟.但是,需要特别指出的是,汶川地震的向南终止绝对不是初始应力场的方位发生变化造成的,因为在汶川地震震中附近一直向其西南方向很大的范围里,应力场方位几乎是相同的(Wu et al., 2009).

此外,本文在模拟过程中,将龙门山断裂带的逆冲兼走滑性质的断层简化为纯走滑断层,将实际的三维断层结构简化为二维.这种简化的合理性主要是：首先从地图上看,汶川地震的破裂是由南向北传播,所以从整体上看该地震可以等价为一个走滑型事件.其次,在地震的起始阶段,破裂表现为逆冲型,但随着时间的推移,走滑特征越来越明显,特别是到破裂的后期,几乎就是纯走滑性质.所以,文中抓住汶川地震单侧破裂这个根本性的科学问题,将其简化处理,这不会从本质上影响我们对这个问题的认识.当然,今后可以利用三维模型,对破裂过程进行更为符合实际的模拟,对文中的结果进行进一步的验证. 顺便指出,震源破裂过程不对称性的运动学研究成果非常丰富(高原等,1997； Gao et al., 2000),这些结果对于深入研究震源破裂的动力学过程提供非常重要的参考.

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

利用有限单元计算方法,可以模拟断层的自发破裂动力学过程.尽管对于断层两侧介质性质不同的破裂问题,模拟中会出现Adams不稳定现象,但通过正则化可以得到稳定的数值解.

龙门山断裂带的西北侧(青藏高原)介质柔软,而东南侧(四川盆地)介质刚硬.由于断层两侧介质物性的不同,导致了汶川地震破裂的不对称性空间分布.

模拟结果表明,汶川地震在破裂过程中,断层面上的正应力在东北方向(破裂的正方向)随着传播距离的增大而不断减小,位错速率随着破裂的传播距离而越来越大,其脉冲变得越来越尖锐,即产生了能够自己放大、自己愈合、自行维持的Weertman脉冲,所以汶川地震在东北方向破裂的距离很大.但在西南方向,破裂过程中断层面上的正应力在不断增大,不利于破裂继续扩展.所以,在龙门山断裂带及附近区域的应力环境不变时,断裂带两侧的物性差异造成了汶川大地震中发生了单侧传播现象.研究断层两侧物性差异影响断层的破裂行为,对于深入认识地震动力学过程、分析未来地震的震级大小、预测地震灾害等有重要的科学意义.

两位审稿专家提出了十分宝贵的建议,在此表示衷心的感谢！