地震研究的终极目标是对地震事件的发生进行精准预报, 但地震的发生并非“无源之水, 无本之木”, 而是断层系统长期演化的结果.研究地震成因, 眼界不应仅仅局限在动态破裂过程(同震阶段), 而应该考虑整个地震周期(包括同震、震后、震间及震前)的断层行为是如何发展演变的.自20世纪以来, 随着观测技术、理论研究及物理学实验的不断发展, 人们对地震演化物理机制研究已经取得长足进步.地震的演化过程就是断层在不同状态下克服摩擦的过程, 准确描述断层摩擦行为的摩擦准则就成为研究者关注的重点.事实上找到这种摩擦准则非常困难, 因为需要该准则既能描述断层长达数百至数千年的长期准静态缓慢滑动, 又能描述地震发生时几十秒至数分钟的瞬间动态破裂过程.到20世纪七八十年代, 研究者(如, Dieterich, 1978, 1979; Ruina, 1980, 1983)经过大量岩石学物理实验总结出一个能够描述断层摩擦中滑移弱化、速率弱化等现象的摩擦准则, 即速率-状态摩擦定律(Rate-State dependent Friction law, RSF).RSF为我们提供了一个基本框架以对断层内部复杂属性进行定量描述.具体地讲, 该摩擦定律表达了断层滑移过程中, 摩擦与观测到的瞬时速率、状态变量之间的相关性.

基于速率-状态摩擦定律对断层演化进行数值建模时, 一般将断层划分为若干子断层, 每个子断层应力既可以表示为初始应力与其他子断层对其影响之和, 同时又要遵从RSF准则, 将两者联立再采用微分方程数值解法进行求解.早期的数值模型一般仅适用于滑动速率非常小(即准静态)的条件下, 但在该模型同震阶段滑动速率会趋向无穷大.Rice(1993)在RSF中引入一个黏滞项近似表示地震波能量的损失, 从而避免了出现滑动速率无穷大的现象, 能够成功模拟出地震循环周期性黏滑现象, 这一模型称为准动力学模型.RSF近些年来已经被广泛用于地震学研究的各个方面(张龙等, 2013), 如Parkfield地震循环的数值模拟(Barbot et al., 2012), 应力扰动对发震时间的影响(Gallovič, 2008; Perfettini et al., 2003; 解孟雨和史保平, 2016a), 俯冲断层上非震滑移行为研究(Liu and Rice, 2005, 2007), 慢地震研究(Shibazaki and Iio, 2003; Shibazaki and Shimamoto, 2007; 解孟雨和史保平, 2016b), 余震触发机制(Dieterich, 1994), 震后形变机制(付真等, 2013)等等方面.而在其他一些涉及到摩擦行为研究的非地震领域, 如山体滑坡的研究中, 该摩擦定律亦被广泛应用(Helmstetter et al., 2004).

速率强化(velocity-strengthening)是指断层摩擦应力随速率增加而增加的行为, 对应着断层摩擦参数a>b的区域, 在该区域会促进产生稳定的滑动或蠕动, 不会产生自发的成核过程, 动态破裂传播至该区域也会停止(Scholz, 1998).震间期在地层浅部的蠕动现象(Lyons et al., 2002)、大地震后浅地表震后余滑(Hsu et al., 2006)以及浅层地表的地震活动性缺失(Shearer et al., 2005)等现象都间接表明在浅层地表存在速率强化区.实验室岩石学实验中, 在低正应力下的岩石摩擦通常由于松散断层泥的影响而表现出速度强化行为(Marone et al., 1991; Marone, 1998).另一方面, 判断断层是否会发生失稳的关键参数—临界刚度系数k_c—也受制于正应力的影响.Ariyoshi等(2007)的研究表明正应力是控制大地震震后滑移传播最重要的因素, 而震后滑移又与余震有紧密联系.由此可见, 正应力是影响断层演化以及地震活动性的关键参数.因此深入研究正应力影响下的断层演化行为有助于了解其对断层演化参数、地震孕育过程、震后滑移传播等方面的影响, 而这些方面直接关系到断层附近强地面运动水平和地震危险性评估.

本文以二维半空间垂直走滑断层为研究对象, 首先阐述建立基于RSF的准动力学数值模型的方法和原理, 给出模型参数, 通过比较两种正应力随深度变化模型的模拟结果, 研究浅层正应力变化对断层演化参数、地震孕育过程、震后滑移传播等方面的影响.

考虑一个半空间垂直纯走滑断层, 所研究问题就变为断层性质随深度变化的二维反平面断层问题.将断层面划分为N_L×N_W个子断层, 断层周围则以背景场板块运动速率V_pl沿走向方向运动.在任意时间节点, 第i个子断层上的摩擦应力τ_i可以表示为:

(1)

其中, τ₀为初始应力, d_j为第j个子断层上的滑动位移, K_ij则是弹性静力学核, 表示第j个子断层上的单位位移在第i个子断层上引起的应力变化.G和β分别为剪切模量和S波速度.(1)式等号右边第三项即为黏滞项, 近似表示地震波能量的损失.第i个子断层的摩擦应力还要遵从速率-状态摩擦准则:

(2)

(3a)

(3b)

其中, σ_i为有效正应力, θ_i是与时间和滑动速率相关的状态变量, 与时间t有相同的量纲, 具体地描述了断层内部和接触时间相关的老化过程.μ_i^*和V^*分别是参考摩擦系数和滑动速率.a_i和b_i为实验所测得断层摩擦参数.D_c为特征滑动距离, 表示滑块受扰动后恢复到稳定状态所需滑动距离.(3a)和(3b)式均被称为状态演化函数, 其中(3a)式由Dieterich(1978, 1979)给出, 被称为慢度准则(Slowness Law)或者老化准则(Ageing Law).(3b)式由Ruina(1980, 1983)给出, 被称为滑动准则(Slip Law).两种状态演化函数各有千秋(Marone, 1998), 我们曾分别利用两种状态演化函数进行过数值模拟, 结果显示在相同子断层划分条件下, 滑动准则能够产生稳定的周期性地震事件, 而在慢度准则下模型会产生不稳定的慢地震事件, 这可能是因为慢度准则对断层离散的网格尺寸更加敏感产生的(Mitsui and Hirahara, 2011).如果想在慢度准则下产生稳定周期地震事件, 必须减小子断层尺寸, 增加子断层数目, 这样一来计算时间会大幅增加.因此, 综合考虑模型的有效性和计算效率, 本文选择滑动准则进行模拟.在稳定状态下, θ_ss=D_c/V_i, 此时(2)式可简化为:

(4)

从(4)式中可以看到, a-b控制着摩擦滑动的稳定性(张慧茜等, 2014).当a-b>0, 摩擦应力随速率增加而增加, 被称为速率强化; 当a-b < 0, 被称为速率弱化, 摩擦应力随速率增加而减小, 因此速率弱化区产生不稳定黏滑或地震, 而速率强化区则会产生稳定的滑动或蠕动.

为计算方便, 我们对参数进行简化, 设Θ=ln(θV_pl/D_c), η=G/(2β), 让(1)式和(2)式相等, 并对滑动速率V_i求导, 然后与状态演化函数(3b)式联立, 得到:

(5)

(5) 式代表的微分方程系统即为断层准动力学模型.给定合适的初始条件, 可采用自适应时间步长的四阶Runge-Kutta积分算法求解(Perfettini et al., 2003).求解过程中最耗费时间的计算是子断层之间相互作用部分, 而由于在水平方向上的速率和应力场是周期性的, 该部分可以写为卷积形式, 因此可以采用二维快速傅里叶变换算法进行求解(Dublanchet et al., 2013).许多研究者曾给出过K_ij的解析表达式(Perfettini et al., 2003; Kato, 2002), 本文使用基于Andrews(1974)表示定理得到的波数域内弹性静力学核表达式.Andrews的方法无法给出自由表面的影响, 为解决这一缺陷, 我们采用Gallovič(2008)给出的镜像法策略, 即将计算范围设定为2L×4W的区域, 实际断层面及其镜像面位于计算区域中间对称位置, 其他区域均以背景场滑动速率V_pl运动, 具体见图 1.图 1中x为断层走向方向, z为垂向方向.

图 1

Fig. 1

图 1 二维准动力学数值模拟计算区域示意图 L和W分别为断层长度和宽度, x为断层走向方向, z为垂向方向. Fig. 1 Computational domain for the 2-D quasi-dynamic numerical simulation L and W are the length and width of the fault, respectively, x represents the direction of strike, and z represents the down-dip direction.

表 1给出了本文模拟中所采用的相关断层参数, 其中断层长宽为一个中强地震断层尺度, 背景场加载速率V_pl设定为平均板块收敛速率32 mm·a^-1≈10^-9 m·s^-1(Lapusta and Liu, 2009).剪切模量、剪切波速、特征滑动距离、参考摩擦系数、参考滑移速率的参数值均参考其他学者模型参数的取值(Gallovič, 2008; Perfettini et al., 2003; Lapusta and Liu, 2009).摩擦参数a和b的数值依赖于环境因素而随深度变化, 温度或断层泥从固结到松散的转换等因素都会影响其变化(Marone, 1998).本文采用Gallovič(2008)使用的a和b随深度变化模型, 在该模型中a值除了随深度变化, 在断层两侧边缘亦有横向变化, a和a-b的值空间分布如图 2所示.断层的浅部、深部及两侧边缘部分显示速率强化性质(a-b>0), 而在中间部分断层具有速率弱化特性(a-b < 0), 该区域即为孕震区, 在该区域a-b=-0.004.为了研究不同正应力影响下的断层行为, 我们给出两种正应力随深度变化模式:第一种为常数模式, 有效正应力设定为75 MPa不随深度变化.第二种模式中, 考虑到在地表断层的剪切加热对孔隙压力影响可以忽略不计, 地表有效正应力设定为3 MPa, 自地表至6 km深度线性增加至75 MPa, 假定在6 km深度以下孔隙压力与地壳压应力以相同比率随深度增加, 则在6 km深度以下有效正应力恒定为75 MPa.两种有效正应力变化模型的对比如图 2所示.后文中, 我们称第一种正应力常数模型为模型Ⅰ, 称第二种模型为模型Ⅱ.

表 1 Table 1

表 1 模型参数表 Table 1 Parameters used in this study 参数名 物理意义 数值 L 断层长度 40 km W 断层宽度 20 km NL 走向方向子断层数目 160 NW 垂向方向子断层数目 80 G 剪切模量 33 GPa β 剪切波速度 3 km·s-1 Vpl 远场加载速率 10-9 m·s-1 Dc 特征滑动距离 4 cm μi* 参考摩擦系数 0.6 V* 参考滑移速率 1 μm·s-1

表 1 模型参数表 Table 1 Parameters used in this study

图 2

Fig. 2

图 2 摩擦参数a、a-b及两种有效正应力在断层面上的分布示意图 Fig. 2 Distribution of friction parameters a, a-b and two types of effective normal stress model on fault plane

Rice和Ruina(1983)基于弹簧-滑块模型对断层的稳定性分析认为, 断层若发生失稳则要求等效弹性刚度k满足k < k_c, 其中k_c为临界刚度系数, k_c=σ(b-a)/D_c.在RSF理论框架下, 对于一个给定长度L的断层, k < k_c等价于L>L_c, L_c被称为成核尺度, 实际表示了能够发生断层失稳的最小尺度.多位学者都曾给出过L_c的理论表达式(Rice and Ruina, 1983; Dieterich, 1992; Ampuero and Rubin, 2008), 本文由于采用的是滑动准则下的状态演化函数, 因此采用Rice(1993)给出的表达式L_c= (2/π)[G/k_c].本研究中L_c理论估算值约为2.8 km.数值模拟中, 离散化后的子断层尺寸ΔL必须远小于成核尺度L_c, 否则会导致个别子断层单独失稳的情况发生, 进而影响断层成核过程(Rice, 1993), 但是在实际模拟中“远小于”这一条件并无定量化约束, 甚至有学者给出的ΔL/L_c比值达到0.6(Shibazaki and Iio, 2003).采用Perfettini等(2003)给出的约束条件, L_c≥10ΔL, 对子断层离散化进行约束.走向和垂向子断层数目分别设置为160和80, 这样ΔL=250 m, 能够满足L_c≥10ΔL的条件.

本文中所有数值计算均基于中国地震局地壳应力研究所地壳运动与遥感应用研究室所属的IBM Bladecenters HS23刀片服务器集群(14台), 每一个刀片机配置有12个2 GHz至强处理器, 单个服务器总内存达256 GB, 能够满足计算需要.实际计算过程显示, 在现有硬件条件下, 每模拟一次大地震断层失稳事件计算时间需要3~4天.我们选择以断层稳定状态下的速率和状态变量作为初始条件, 在速率弱化区给定V_i=V_pl/1000, 而在速率强化区给定V_i=V_pl, 取θ_i= D_c/V_i.模拟结果显示, 在该初始条件下断层能够产生稳定性的周期性地震重复过程.

图 3、图 4分别展示了模型Ⅰ和模型Ⅱ两种模型下, 断层面上一个地震周期内不同演化阶段的速率V和摩擦应力τ分布.模拟结果显示, 我们的数值模型在给定上述模型参数和约束条件下, 能够完整模拟出地震周期的全过程, 包括震间、震前、同震以及震后多个特征阶段.

图 3

Fig. 3

图 3 基于模型Ⅰ模拟得到的一个地震周期内断层面速率和摩擦应力演化 Fig. 3 Snapshots of spatial slip velocity and stress distribution during a typical earthquake cycle for Model Ⅰ

图 4

Fig. 4

图 4 基于模型Ⅱ模拟得到的一个地震周期内断层面速率和摩擦应力演化 Fig. 4 Snapshots of spatial slip velocity and stress distribution during a typical earthquake cycle for Model Ⅱ

在震间期, 速率弱化区中速率水平是背景场加载速率V_pl的百分之一乃至千分之一, 可以被视为闭锁静止状态.在速率弱化区周围的速率强化区, 断层基本以V_pl速率水平进行运动, 这种缓慢滑动会在滑动区和闭锁区接触边缘产生应力集中, 导致闭锁区周边的稳定滑动慢慢向速率弱化区渗透扩展, 导致速率弱化区内越来越多的区域加载到V_pl水平.当断层闭锁区域逐渐减小到几乎消失时, 断层中心区域的滑动速率继续增加, 开始进入成核阶段.成核完成后, 动态破裂开始从成核区向断层两侧对称性传播, 并很快拓展至几乎整个断层面.在动态破裂传播过程中, 破裂前缘会形成较高的应力前锋, 其经过的区域应力水平迅速降低.动态破裂后, 应力前锋渗透至深部速率强化区, 引起余滑的产生.由于应力前锋在速率强化区传播速度远远低于动态破裂速度, 因此在震后很长时间内速率强化区应力水平要高于速率弱化区应力水平.仔细观察会发现, 地震之后速率弱化区内应力水平与滑移速率同步开始下降, 但一定时间后应力水平开始上升, 而滑动速率则继续下降.这一现象与Miyazaki等(2004)实际观测是相符的, Helmstetter和Shaw(2009)认为速率弱化区内应力水平的后期上升是由于周边速率强化区内大量余滑的加载效应.

除上述共同演化特征外, 两种正应力变化模型的模拟结果也存在不同之处, 其中最明显区别是在模型Ⅱ下, 浅层摩擦应力明显小于模型Ⅰ, 且在同震阶段的动态破裂能够传播至自由表面.而在模型Ⅰ下, 由于浅层相对较大的正应力导致速率强化行为效应增强, 动态破裂在浅层速率强化区停止, 并没有到达自由表面(Lapusta and Liu, 2009).这一点从我们给出的断层上x=20 km(x为断层走向方向)处各点滑动量随时间的变化图中也可以验证(图 5).图 5给出了两种模型下5次地震演化过程.

图 5

Fig. 5

图 5 断层上x=20 km处各点滑动量随时间的变化图 上图为模型Ⅰ结果, 下图为模型Ⅱ结果.蓝色实线代表每隔5年震间期对应的滑动量积累, 绿色实线为主震后一年内震后阶段每7天对应的滑动量积累, 红色实线是同震阶段当最大滑动速率达到1 mm·s-1后每2 s代表的滑动量积累 Fig. 5 Slip accumulation along the line x=20 km The above figure shows the result of model Ⅰ, and the figure below shows the result of model Ⅱ.Blue lines are plotted every 5 years, representing slip accumulation during interseismic period.Green lines illustrate rapid postseismic afterslip which is plotted every 7 days.Red lines illustrate coseismic slip which are plotted every 2 s when maximum slip velocity over the fault exceeds 1 mm·s-1.

图 5中显示在速率强化区主要通过震间和震后阶段累积滑移, 而滑移弱化区主要通过同震过程累积滑移.蓝色实线的线间距表明震间期断层在速度强化区域稳定累积滑移, 绿色实线显示在震后阶段快速累积了震后滑移.密集的红色实线代表成核阶段的过程, 稀疏的红色实线代表动态破裂传播过程(Lapusta and Liu, 2009).可以看到, 两种模型动态破裂在深度方向上均渗透至16 km深度左右, 成核区的深度也都在7~10 km左右.仔细观察成核过程可以发现, 在模型Ⅰ下5次地震事件成核阶段滑移累积是自下而上, 而在模型Ⅱ下成核阶段滑移累积除第一次事件是自上而下外, 随后的四次事件又恢复为自下而上, 但与模型Ⅰ有区别, 说明断层正应力在浅层的非均匀性能够显著影响其成核过程.Kaneko等(2008)发现平均正应力大的情况下反而其成核尺度更大, 而这与Ampuero和Rubin(2008)理论公式中正应力与成核尺度成反比的结论相悖.Lapusta和Liu(2009)的研究显示, 一个占断层面不足1%大小的高正应力区域就能非常显著地影响断层长期演化并改变断层成核位置.我们发现, 尽管浅层正应力存在差别, 两种模型下只是成核过程有明显差异, 成核区尺度并无显著差异, 因此本文研究结论倾向于支持Kaneko等(2008)的数值模拟结论, 也就是说正应力的非均匀性对地震成核有显著影响, 且该影响有时与直觉相悖.因此在研究涉及到非均匀正应力条件下的地震成核问题时, 不能凭经验简单套用理论公式, 应该由数值模拟进行验证.

为更好对比两种模型下断层演化的差异, 我们给出了断层面三个位置上子断层速率和位移随时间变化图, 这三处位置分别为P1(x=10 km, z=0.5 km)、P2(x=10 km, z=9 km)和P3(x=10 km, z=16.5 km), 其中P2位于速率弱化区而其余两点位于速率强化区.

位于速率弱化区内的P2点显示出明显的黏滑(stick slip)行为, 即周期性静止-失稳, 失稳过程中最大滑移速率超过1 m·s^-1, 其滑移的累积近似于阶梯型.P1和P3点大部分时间以接近V_pl速率水平运动, 模型Ⅰ下P1点最大滑移速率达到0.007 m·s^-1, 模型Ⅱ下P1点最大滑移速率达到0.06 m·s^-1, 两者相差近一个数量级.在速率强化区深部, 两种模型下P3点最大滑移速率均为6.5×10^-5 m·s^-1左右, 并且随着深度增加最大滑移速率减小.经过测算, 模型Ⅰ下平均地震周期T_r1=82.9 a左右, 在模型Ⅱ下平均地震周期T_r2=82.1 a左右.Kato和Tullis (2003)的研究显示, 无论是在慢度准则、滑动准则或者其自己提出的复合状态演化准则(Kato and Tullis, 2001)下, 地震复发周期均正比于断层的相对刚度(k_c/k).模型Ⅰ和模型Ⅱ的断层系统刚度系数是相同的, 而前者的临界刚度系数k_c要大于后者(前者的平均正应力更大), 因此理论上推测T_r1>T_r2, 这与模拟结果相符.接下来我们研究单次地震事件地震矩和平均静态应力降的差异.

(6)

(7)

我们规定断层失稳的速率阀值为0.01 m·s^-1, 若速率超过该阀值即认为动态破裂发生.分别按照(6)式和(7)式估算地震矩M₀和断层平均静态应力降Δτ, 其中Ω为发生动态破裂的断层面积, t_ini和t_end分别为动态破裂起始时间和终止时间(Lapusta and Liu, 2009).经过计算, 模型Ⅰ和模型Ⅱ下地震事件释放的地震矩分别为2.60×10¹⁹ Nm、2.58×10¹⁹ Nm, 换算为矩震级均为M_w=6.88左右; 两种模型下最大同震滑移分别为2.34 m和2.32 m, 平均同震滑移量分别为1.27 m和1.21 m.用远场滑移速率与地震平均复发周期相乘就可以估算两种模型一次地震周期内滑移量累积分别为2.61 m和2.58 m, 因此平均同震滑移量占总滑动量百分比分别为49%和47%.图 7给出了两种模型下同震滑移分布以及x=30 km处滑移和静态应力降随深度的变化情况.

图 7

Fig. 7

图 7 模型Ⅰ和模型Ⅱ下同震滑移分布以及x=30 km处同震滑移和静态应力降随深度的变化情况 Fig. 7 The distribution of coseismic slip of model Ⅰ and model Ⅱ and the variation of the coseismic slip and static stress drop along the depth at x=30 km

从图 7中可以看到, 两种模型下的同震滑移分布总体特征非常相似, 在成核区外等值线均呈椭圆形向外拓展.模型Ⅱ由于破裂至自由表面处, 因此在该处存在同震位移, 而模型Ⅰ在地表则无同震滑移.模型Ⅰ在速率弱化区(深度范围4~14 km)静态应力降为常数, 而模型Ⅱ显然受到浅层正应力变化影响, 其常数静态应力降深度范围为6~14 km.利用(6)式求得模型Ⅰ和模型Ⅱ下地震事件平均静态应力降分别为18.46 MPa和15.54 MPa.总体上讲, 虽然模型Ⅱ在浅层同震滑动量要高于模型Ⅰ, 但是其总体地震矩释放量和平均静态应力降相差并不明显, 浅层正应力对常数应力降深度范围存在影响.从空间域比较同震滑移模型很难看出两者差别, 我们将在讨论部分讨论在波数域下两种同震滑移模型的差异及其可能对地表强地面运动产生的影响.

大量研究表明, 无论是持续时间还是衰减特征, 余震与震后形变具有相似特征(Perfettini and Avouac, 2007), 由于震后地表形变主要来源于断层震后滑移(又称为余滑, afterslip), 因此许多研究者开始关注震后库仑应力变化与余震之间的关系, 并且越来越多的证据表明震后滑移对余震的触发起着关键作用, 尤其在震后初期数百天内, 震后滑移对余震的影响最重要(Perfettini and Avouac, 2007; Hsu et al., 2006; Bourouis and Bernard, 2007).因此, 研究震后滑移传播的空间变化及其影响因素能使我们更好评估余震活动性和余震迁移的物理机制.

图 6

Fig. 6

图 6 断层上P1(上)、P2(中)和P3(下)点速率和位移随时间变化图 黑色实线为速率-时间变化, 红色实线为位移-时间曲线. Fig. 6 Long-term histories of velocity and slip at three representative fault locations, P1 (upper), P2 (middle) and P3 (lower) Black lines represent velocity-time histories, and the red lines represent slip-time curve.

震后滑移在速率强化区内的传播速度远低于剪切波速度, 我们以主震后子断层剪应力达到峰值时刻, 即应力集中前缘(Stress Concentration Front, SCF)到达时刻, 为震后滑移开始时刻, 记为t_max(Ariyoshi et al., 2007).对于靠近自由表面速率强化区内的点, 由于动态破裂传播至该区域, 其摩擦应力值在同震过程中即达到最大, 因此我们将只讨论在断层深部(深度范围:16.5~20 km)震后滑移传播.以各点到震源的距离除以t_max作为震后余滑传播的平均速率, 记为V_post.

图 8给出了两种模型下主震后t_max和V_post在空间分布情况.从总体特征上两种模型下的t_max和V_post在空间分布非常类似, 两者均呈辐射状向断层深部传播, 断层中间位置传播最快, 在主震后4个月到达断层最底部, 而两翼传播则相对较慢.震后滑移传播平均速率随深度迅速衰减, 且深度越深衰减越快.Jiang和Lapusta(2016)的研究显示, 震后滑移传播与子断层震后滑移速率、震间期滑移速率、aσ以及静态应力降有关.两种模型下, 这几项参数在断层深部类似, 因此t_max和V_post在空间分布上呈现出相似特征, 这种相似性说明浅层正应力变化对断层深部余滑传播影响不大.

图 8

Fig. 8

图 8 模型Ⅰ和模型Ⅱ下震后滑移开始时间tmax及平均传播速度Vpost空间分布 Fig. 8 The spatial distribution of tmax and Vpost after the earthquake in model Ⅰ and model Ⅱ

由于强地面运动尤其是高频地震波的强地面运动对于同震滑移分布的不均匀性非常敏感, 因此断层破裂面上滑移分布的不均匀性作为定量化分析断层面滑移分布重要特征就成为预测发震断层强地面运动水平的重要参数.拐角波数K_c控制着滑移分布不均匀性的强弱, 对一个给定大小的断层尺度, 拐角波数越小就意味着滑移分布越均匀, 滑移面越光滑, 而拐角波数越大就意味着滑移分布越粗糙, 不均匀性越强(Somerville et al., 1999; Mai and Beroza, 2002).

我们分别求取了两种模型下同震滑移模型(图 9)沿走向和下倾方向滑移分布的拐角波数K_cx和K_cz.在计算过程中, 先将同震滑移分布通过二维傅里叶变换转化到波数域得到波数谱(D_obs), 并将其与理论的波数谱(D_mod)进行比较, 最终通过最小化参数E来求取拐角波数.E的表达式如下:

(8)

图 9

Fig. 9

图 9 模型Ⅰ、模型Ⅱ空间拐角波数计算示意图 其中黑色实线为理论的波数谱, 圆圈为滑移模型通过傅里叶变换得到的波数谱, 虚线对应的波数是表达式E取最小值时的波数值, 即拐角波数. Fig. 9 Calculation of the corner wavenumber for model Ⅰ and model Ⅱ The black solid line indicates the theoretical 1-D K-2 model; the circles indicate the slip spectrum and the dashed lines denote the values of corner wavenumber

理论的波数谱的表达式为K^-2模型(Herrero and Bernard, 1994):

(9)

式中K为波数, K_N为尼奎斯特波数.图 9给出了两种模型同震滑移分布在走向和下倾方向上空间拐角波数的计算值.从图 9可以看到无论是走向方向还是下倾方向, 两种滑移模型的傅氏振幅谱与理论K^-2模型傅氏振幅谱均符合较好, 说明两种模型均符合地震同震滑移模型的运动学特征.Causse等(2010)的研究发现高频地震波的强地面运动水平同拐角波数存在正相关关系, 即峰值加速度PGA∝K_c^3/2, 而就拐角波数值的大小而言, 模型Ⅱ的拐角波数值明显高于模型Ⅰ, 这也就意味着在其他参数一致的条件下, 模型Ⅱ最终应该产生高于模型Ⅰ的高频强地面运动水平.

在本研究中我们使用的准动力学数值模型与动力学模型的最大差别即在于前者将地震波辐射近似为一个黏滞项, 而忽略了瞬态地震波传播对动态破裂的影响, 这一差别导致两者模拟结果虽然在地震演化形态上存在近似, 但是在演化数量上存在明显差别.研究表明准动力学数值模型下断层发震的平均破裂速度、最大滑动速率、最终滑动量、应力降等都要低于动力学模型的结果(Lapusta et al., 2000; Lapusta and Liu, 2009).Lapusta和Liu(2009)研究显示, 可以通过增加剪切波速的手段减小准动力学模型与真实动力学模型结果的差别, 但这样做的负面效果是计算量的增加.我们认为应当根据具体问题来选择模型以及合适的剪切波速度, 这样才能在保证结果准确前提下有效提高计算效率.

RSF由于其本身来源于室内的岩石物理学实验, 在该模型中一些重要本构参数的取值是否能直接应用于天然地震仍存在较大争议, 另外这些本构参数是否为常数又或者与构造环境有关系也存在不确定性.本文中, 我们仅考虑了两种正应力变化模型, 且仅在浅地表存在差异, 也有一些研究者在数值模型中考虑过有效正应力随深度线性增加的情况(Kaneko et al., 2013), 除了有效正应力外, 其他模型参数对断层演化带来的影响还需要研究.

本文基于结合RSF的准动力学数值模型, 以二维半空间垂直走滑断层为研究对象, 通过比较两种正应力随深度变化模型的模拟结果, 研究了浅层正应力变化对断层演化参数、地震孕育过程、震后滑移传播等方面的影响.使用基于Andrews(1974)表示定理得到的弹性静力学核计算子断层间相互作用, 并通过镜像法得到存在自由表面影响的结果.结果显示, 我们的数值模型在给定模型参数和约束条件下, 能够完整模拟出地震周期中震间、震前、同震以及震后多个特征阶段.常数正应力模型下, 由于浅层相对较大的正应力导致速率强化行为效应增强, 动态破裂在浅层速率强化区停止, 而在浅层变化正应力模型下动态破裂可以传播至自由表面, 从而导致更高的最大滑移速率和同震滑移量.另外我们发现两种模型下成核尺度近似, 但成核过程有明显差异.经过计算, 在一次地震事件中两种模型下的地震矩、地震周期、平均应力降和同震滑移分布特征的值存在差异但差别不大, 浅层正应力变化对常数应力降深度范围存在影响.两种模型下, 深部震后滑移开始时刻及其空间传播速率在空间分布上呈现出相似特征, 说明浅层正应力变化对断层深部余滑传播影响有限.我们分别求取了两种模型下同震滑移模型沿走向和下倾方向滑移分布的拐角波数, 结果显示, 两种滑移模型的傅氏振幅谱与理论K^-2模型傅氏振幅谱均符合较好, 说明两种模型均符合地震同震滑移模型的运动学特征.模型Ⅱ的拐角波数值明显高于模型Ⅰ, 意味着在其他参数一致的条件下, 模型Ⅱ最终应该产生高于模型Ⅰ的高频强地面运动水平.综上所述, 我们的研究结果表明浅层正应力变化对地震成核过程、浅层滑动速率及同震滑移均有影响, 但对深部震后余滑传播影响有限.基于RSF的准动力学模型在地震学领域有着广泛的应用前景, 选用不同的模型参数对最终结果存在显著影响, 如何将观测数据与物理模型结合对断层演化参数进行精确约束, 合理判断断层演化所处状态, 进而进行准确的地震危险性和危害性分析是未来需要解决的科学问题.