地球物理学报  2010, Vol. 53 Issue (1): 102-110   PDF    
库仑应力计算及应用过程中若干问题的讨论-以汶川地震为例
石耀霖1 , 曹建玲1,2     
1. 中国科学院研究生院计算地球动力学实验室, 北京 100049;
2. 中国地震局地震预测研究所, 北京 100036
摘要: 大地震发生后,估计后续地震发展趋势是人们关心的问题.目前,学术界经常利用大地震造成的库仑应力变化研讨大震对后续地震的影响,但库仑应力的计算和应用中尚有一些问题被忽视而有待探讨.本文修正了传统的库仑应力计算中沿地震破裂面滑动方向计算剪应力变化的近似方法,考虑震后主应力方向可能改变对剪应力变化量计算的影响,对改进方法和传统方法获得的结果进行了对比分析.本文基于龙门山地区背景应力场测量资料,根据不同研究者反演的地震破裂模型,计算了汶川地震造成的静态库仑应力变化,考察了不同地震破裂模型下库仑应力分布差异.在断层面附近,使用改进方法与传统方法计算的库仑应力分布差异相当大;利用多个地震破裂模型计算大地震库仑应力判断大震后余震发展趋势时,要注意破裂模型不确定性对危险性估计的影响.
关键词: 库仑应力      修正方法      应力触发      地震危险性     
Some aspects in static stress change calculation-case study on Wenchuan earthquake
SHI Yao-Lin1, CAO Jian-Ling1,2     
1. Laboratory of Computational Geodynamics, Graduate University of Chinese Academy of Sciences, Beijing 100049, China;
2. Institute of Earthquake Science, China Earthquake Administration, Beijing 100036, China
Abstract: When a major earthquake occurs, how aftershocks distribute catch the attention of the public. Recently, the concept of Coulomb failure stress is used frequently to study the earthquake triggering of main shock. However, there are some pendent problems in the calculation and application of Coulomb stress. In conventional calculation of Coulomb stress, the shear stress change is approximated to be the changed value along the slip on the rupture of main shock. In fact, the stress orientations may change after the event. It affects the shear stress on ruptures. We improved the formula of Coulomb stress and made comparison. We applied the improved method to Wenchuan earthquake. According to the stress measurement in Longmenshan and the three rupture models, we calculated the static Coulomb stress of Wenchuan earthquake. The results show big differences between conventional and improved calculation around the fault. Our investigation reminds us to concern the uncertainty of the rupture models when we use Coulomb stress to predict the aftershocks..
Key words: Coulomb failure stress      Improved method      Earthquake trigger      Potential seismicity     
1 引言

近年来, 人们广泛采用静态库仑破裂应力变化(ΔCFS)来讨论余震与主震关系以及强震间的应力触发作用[1~10].Stein等考察了1979年加州地震造成沿破裂面的库仑应力变化和后续余震分布关系[1].King等计算了1992年Landers7.4级地震造成破裂面附近沿最优方向的库仑应力变化, 并考察了该地震造成周边断层面上的应力变化[2].Stein等计算了在1939~1992年间发生在土耳其North Anatolian断裂带10个6.7级以上地震的库仑破裂应力演化过程, 发现其中的90%的地震是被先前地震触发[11].沈正康等利用黏弹性模型, 模拟了1937年M7.5花石峡地震、1963年Ms7.1都兰地震、1973年Ms7.3玛尼地震、1997年Mw7.5玛尼地震和2001年Mw7.8可可西里地震之间的应力转移过程, 认为前面的4个地震可能触发了2001年可可西里大地震[12].郝平等对2001年昆仑山口西大地震及其后发生的5个强余震进行了库仑破裂应力的研究, 发现其中4个强余震发生在主震造成的库仑破裂应力增加区域[13].万永革等对发生在青藏高原东部20个7级以上大地震应力演化和地震触发作用进行了研究, 发现85%的地震是由于库仑破裂应力增加而被触发[14].万永革等利用黏弹性模型研究了唐山地震的应力触发效应, 发现绝大多数余震发生在库仑破裂应力增加区[15].上述研究表明强震发生造成的静态库仑应力增加有利于附近区域余震以及后续地震的发生, 从而可以通过计算静态库仑应力变化对后续地震趋势进行评估.

许多研究者计算了汶川地震造成的静态库仑破裂应力变化, 对周边断层地震危险性进行了评价[16~18].Parsons等利用USGS (美国地调局)的汶川地震破裂模型计算了周边断层的静态库仑应力, 其中雅安断裂的库仑应力增加可达0.1 MPa[16]. Toda等的计算结果认为鲜水河断裂南端、东昆仑断裂和岷江断裂今后地震危险性增加[17].万永革等利用由GPS和InSAR资料反演的汶川地震破裂模型, 计算了周围断层的静态库仑应力, 结果表明龙门山断裂两端、鲜水河断裂南端、东昆仑等多条断裂库仑应力增加[18].

然而, 在计算静态库仑应力变化过程中, 需要考虑影响计算结果的因素.首先, 区域构造应力场会对库仑应力计算结果产生影响[12, 13, 19~22].King等注意到区域构造应力不仅使得地震发生后断层面上的应力状态更为复杂, 而且还影响到库仑应力分布[2]. Hardebeck等研究了Landers地震和Northridge地震对余震触发作用, 发现考察静态库仑应力变化的方法适用于探讨Landers地震和余震关系, 但对Northridge地震则不适用; 并指出应用这种方法需要考虑到主震发生前的区域应力场和岩石强度等因素[23].其次, 由于反演问题复杂性, 不同研究者反演的地震破裂过程往往会在"大同"的情况下存在"小异".不同地震破裂模型会造成计算的库仑应力分布差异.另外, 断层位置、走向和倾角, 孔隙压力大小等, 都会对库仑应力变化计算结果产生影响.最后, 传统计算地震破裂面上库仑应力方法, 通常假定断层面还会沿主震错动方向活动, 而只考虑沿主震滑动方向剪应力变化, 但余震震源机制解表明, 主震会改变断层面附近主应力方位[24], 而使断层在主震后不会严格按照主震滑动方向活动.本文修改了传统计算库仑应力方法, 不限定剪应力只能作用在原主震错动方向上, 观察与传统算法的差别; 另外, 讨论同一地震不同研究者反演的错动模型对库仑应力计算, 以及对地震活动性探讨的影响.

2 改进的库仑应力计算方法及结果

大地震造成断层面及邻区库仑应力变化被广泛应用于研究其对周围断层的应力触发作用, 和探讨大地震与后续地震关系.通常定义库仑破裂应力Δσc

(1)

其中, Δτ是沿断层面剪应力变化, μ为断层面摩擦系数, Δσn为断层面上正应力变化(鉴于弹性力学与岩石力学中对应力张量的正负定义不同, 本文全部遵循弹性力学定义, 即张应力为正), Δp为孔隙水压变化.如果断层面上有效正压应力减小, 剪切应力增大, 则断层面上库仑应力增加, 断层更容易发生错动.按照库仑应力定义(1), 断层面上库仑应力变化为正值, 表明断层向破裂状态靠拢, 存在滑动危险; 相反, 如果库仑应力为负值, 则表明断层更加安全.

如前所述, 在计算库仑应力之前, 需要指出所求库仑应力是对应于哪一个断层面的.通常, 在地表资料揭示特定断层的走向倾角时, 可以计算该断层面上的库仑应力; 对于通过深部余震震源机制了解其破裂面的情形下, 计算沿余震破裂面的库仑应力; 在研究主震断层及其延伸段的危险性时, 计算主震断层破裂面库仑应力分布; 在断层不发育或不清楚断层分布情况下, 计算沿最优破裂方向库仑应力[1].

在传统计算库仑应力变化方法中, 通常由于不了解构造应力的大小、仅仅通过震源机制了解错动方向, 因此对断层面上剪应力的计算往往要做一些简化, 假定地震发生后震源附近主应力方向不发生改变, 计算库仑应力时, 只考虑断层面上剪应力在主震滑动方向上大小的变化[2].然而研究表明, 地震虽然对远场区域构造应力方位改变甚为微小, 但可以改变震源附近主应力的方位, 余震的T, P, N轴方位和主震会有显著不同[24].本文对以往计算库仑应力方法进行了修正.考虑到地震以后主应力方位可能会发生变化, 不限定断层面按主震错动方向活动, 也就是计算断层面剪应力绝对值变化.

设震前的应力为σij0.在大地震发生后不久, 地震学家就可以求出其震源机制解, 结合野外地质考察资料, 可以判断发震断层产状, 包括断层走向、倾向和倾角.经过简单推导, 可以根据断层面走向S和倾角D求出断层面法向量(n1, n2, n3), 其中, n1=cosa1, n2=cosa2, n3=cosa3, a1, a2, a3分别为断层面法向量与X轴, Y轴, Z轴夹角, 则断层面法向量为

(2)

任意给定一个断层面走向和倾角, 就可以计算沿断层面上的初始应力矢量T0和初始正应力σn0和剪切应力τ0:

(3)

在确定发震断层几何形态和地震破裂面尺寸及错动量后, 可以计算地震造成的扰动应力场Δσij, 最简单的模型是利用弹性半无限空间理论位错模型(Okada位错模型)计算[25, 26], 固然也可以用分层黏弹性模型, 甚至更复杂的模型直接进行数值计算.将Δσij按分量叠加到构造应力场后得到新的应力张量(σij=σij0σij), 根据(3)式计算出任意产状断层面上新的正应力σn和剪切应力大小τ.然后可以计算正应力和剪应力的同震变化量:

(4)

(5)

需要注意的是, 虽然Δσn只与同震变化量Δσij有关, 但Δτ还与初始应力状态σij0有关, 而在初始应力未知的条件下, 往往根据震源机制解的滑动方向S(s1, s2, s3), 近似认为ΔτσijnjSi.这种近似会造成多大的影响许多研究者并没有真正予以考虑和讨论.

主震造成沿破裂面或平行破裂面断层的库仑应力变化可以按(6)式进行计算:

(6)

孔隙水压力的作用可以简化为降低了的断层面等效摩擦系数μ′来考虑.μ′取值通常在0~0.75之间, 平均值约0.4[6].地震位错产生的区域应力变化采用USGS基于Okada表达式开发的边界元程序3D-DEF[27, 28]计算, 杨氏模量取值70 GPa, 泊松比取值0.25.

为了便于比较改进方法与传统方法得到的库仑应力结果, 选取一个简单逆冲断层为例, 分别按这两种方法计算沿断层破裂面库仑应力变化.该逆冲断层长300 km, 宽30 km, 走向180°, 倾角45°, 滑动量1m, 滑动角90°(图 1a).首先按传统不考虑背景构造应力场方法计算沿断层面滑动方向库仑应力变化.然后参照实测结果的量级[29], 假设该逆冲断层处于一个东西向挤压构造环境下, 最大主压应力为10MPa, 其余主应力为零, 用改进方法计算沿断层面库仑应力变化.断层面摩擦系数取0.85, 计算地表情况.由于地震发生不会改变背景应力场, 但主要是影响到断层面附近应力场方向, 所以我们着重观察两种方法获得的断层面附近库仑应力差异.在该断层面一端选取一个10 km长、10 km宽区域, 来比较两种方法的结果差异(图 1b, 图 1c).

图 1 (a)断层在地面迹线及放大的区域所在位置;(b)按传统方法计算的沿地震破裂面滑动方向的库仑应力;c)按改进方法计算的地震破裂面的库仑应力 Fig. 1 (a) Fault plane projected on the surface and the location of the zoomed zone; (b) Coulomb stress change along the slip direction (conventional method); (c) Coulomb stress change along the fault plane (improved method)

在断层两端几公里范围内, 库仑应力分布差异十分显著(图 1).传统方法计算的断层已破裂部分附近库仑应力均为降低, 仅在断层端点以外未破裂部分库仑应力升高; 而使用改进方法计算结果显示库仑应力增加区的范围明显增大, 断层已错断部分靠近端点附近的库仑应力符号与传统方法的结果相反, 原来认为库仑应力降低的区域, 实际应为库仑应力增高.也就是说, 靠近断层端点, 已破裂部分倾向于继续滑动产生更大位错, 未破裂部分更容易发生沿原断裂延伸的破裂.这对靠近端点部位断层的传播和终止会有很大的影响, 值得在今后研究中给予足够重视.

3 汶川地震造成的库仑应力变化

2008年5月12日发生在龙门山断裂的汶川大地震是中国大陆自1976年唐山地震以来发生在人口密集区又一次特大地震, 此次地震造成了当地人民巨大生命和财产损失.这次地震主震面波震级达到了8.0级, 后续余震频发.据中国地震台网测定, 仅震后一个月(截止2008年6月11日)发生在汶川及周边地区的面波震级在4.0以上余震就有208个之多(图 2), 其中包括33个5级以上地震和5个6级以上地震.

图 2 研究区域构造简图及地震分布情况 蓝色线条表示断裂:1.汶川-茂汶断裂;2.映秀-北川断裂;3.彭县-灌县断裂;4.虎牙断裂;5.平武-青川断裂;f6.岷江断裂;7.鲜水河断裂;8.安宁河断裂.红圈是中国地震台网给出的4级以上余震分布.USGS公布的部分余震震源机制解中,红色的对应主震震源机制解. Fig. 2 Tectonic sketch and Wenchuan earthquake sequence Blue lines denote fault: f1. Wenchuan-Maowen fault; f2. Yingxiu-Beichuan fault; f3. Pengxian-Guanxian fault; f4. Huya fault; f5. Pingwu-Qingchuan fault; f6. Minjiang fault; f7. Xianshuihe fault; f8. Anninghe fault. Red circles are earthquakes (Ms≥4) recorded in one month after the main shock. Mechanism solutions given by USGS, the red one is the main shock.

在龙门山地区实测背景应力场的基础上, 按照修正的库仑破裂应力计算方法, 根据国内外相关研究人员反演得到的汶川地震震源破裂模型, 计算了不同地震破裂模型在不同孔隙压力条件下造成的区域静态库仑应力变化, 并讨论不同因素对计算结果影响.根据计算的库仑破裂应力分布结果, 不但按常规对余震分布以及未来邻近区域潜在地震危险性问题作一些初步探讨, 更着重探讨不同作者的不同破裂模型对地震危险性讨论会有什么影响, 这对了解库仑应力讨论地震危险性时包含的不确定性有更明确的认识, 也才能更恰当地运用库仑应力.

3.1 汶川地震破裂模型

在计算地震造成的应力变化时, 需要确定地震破裂面几何形态和破裂面滑动量.一般情况下, 在地震后几分钟到数十分钟时间内, 地震学家可以根据观测到地震波很快确定出震中位置, 在地震发生后几个小时以后, 研究人员就可以根据地震仪记录的波形反演出地震破裂过程和破裂面几何形态.中国科学院地质与地球物理研究所姚振兴等、中国地震局地球物理研究所陈运泰等以及美国地调局(USGS) Ji Chen等都在震后不久公布了各自反演的汶川地震破裂模型(表 1).尽管有些作者根据更丰富的资料随后又对自己的破裂模型进行了改进, 但本文主要是讨论不同模型对库仑应力和后续地震活动性估计可能产生的影响, 因此并没有采用更新的模型.汶川地震发生在龙门山断裂带中部, 龙门山断裂带主要由三条北东走向断裂组成, 分别是汶川-茂汶断裂、映秀-北川断裂和彭县-灌县断裂(图 2).虽然这些断层几何形态存在一些差异, 倾向和倾角也存在分段差异[30], 但它们总体走向与发震断层破裂面走向比较接近.

表 1 汶川地震破裂模型参数 Table 1 Rupture models of Wenchuan earthquake

模型一来自姚振兴、王卫民等[31]给出的汶川地震破裂面, 最大滑动量可达11 m.模型二根据陈运泰等人(个人通信)给出的汶川地震破裂面, 最大滑动量为6.7 m.模型三是美国地调局公布的汶川地震破裂资料, 最大滑动量约为9 m (http://earthquake.usgs.gov/eqcenter/eqinthenews/2008/us2008ryan/finite_fault.php).上述3个地震破裂模型中, 模型一与模型三在破裂面产状和静态滑动量上比较相近, 断层长度差别约100 km, 模型二倾角最大, 破裂面最长, 滑动量最小.3个模型都将破裂面按走向和倾向划分成若干个子单元, 每个子单元有不同滑动量, 模型一和模型三的子单元不仅有不同滑动量, 而且滑动角也各不相同, 模型二的每个子单元滑动角相同但滑动量不同.为考察3个地震破裂模型上述区别会造成库仑应力分布的差异, 在基于弹性半无限空间假设条件下, 本文计算了3个地震破裂模型造成的汶川大地震及邻近区域沿主震破裂面库仑应力变化.

3.2 考虑区域背景应力场的库仑应力计算

在计算断层面上剪应力变化时, 应该考虑原来的应力张量, 即背景构造应力场.岩体压力随着深度增加, 各个主应力量值均会随深度变化, 应力球张量数值随深度增加.在计算库仑应力时, 剪应力在主震前后变化主要与背景应力场偏应力相关, 应力球张量不影响库仑应力的计算.沿龙门山断裂带两侧进行的水压致裂地应力测量结果表明[32], 该地区构造应力场与龙门山断裂带逆冲性质一致, 最大主压应力和中间主压应力水平, 垂直应力最小, 最大主压应力的方向为北西60°.综合该区6个百米深度钻孔水压破裂应力测量结果, 在地表最大水平偏应力(在弹性力学定义中为最小主应力σ3)约为-7.5 MPa, 最小水平偏应力σ2约为-3 MPa, 垂直偏应力σ1为0.这些地应力测量数据仅能代表近地表地应力状态, 所以只能用于计算地表附近库仑应力.

实测主应力需要转换到计算同震应力所用坐标系.如果选取东向和北向分别为坐标系X轴和Y轴, Z轴指向上方, 由于实测最大水平主压应力方向为北西60°, 与X轴正方向夹角α为30°, 将主轴坐标系绕Z轴旋转了30°得到构造应力场6个分量:

(7)

由于在龙门山地区地应力测量深度浅, 测量结果仅代表近地表处地应力状态, 我们只考虑地表库仑应力变化(图 3a, 图 3b).根据模型一计算的沿破裂面库仑破裂应力分布显示, 库仑应力增加大多出现在地震破裂面上和破裂面两端向外延展的区域, 而在破裂面两侧大范围区域库仑应力下降明显.

图 3 根据不同模型计算的地表库仑应力增量 (a) (b)是根据模型已计算的结果,(c) (d)是根据模型二计算的结果,(e) (f)是根据模型三计算的结果;(a) (c) (e)为等效摩擦系数取0.85的情况,(b)(d)(f)为等效摩擦系数取0.4的情况;黑色圆圈为4级以上余震分布. Fig. 3 Coulomb stress change on the surface calculated from three fault models (a)(b) Results from the first model; (c)(d) Results from the second model; (e)(f) Results from the third model. (a)(c)(e) are cases that the friction coefficient is 0.85; (b)(d)(f) are cases that the friction coefficient is 0.4. Black circles are aftershocks (Ms≥4).

为了比较孔隙水压力对库仑应力计算值的影响, 考虑等效摩擦系数μ′分别取0.85和0.4两种情况, 根据上述3个地震破裂模型和公式(4)、(5)、(6)计算了汶川大地震造成的库仑破裂应力变化.

图 3的计算结果表明, 同一模型等效摩擦系数不同时, 虽然计算的库仑应力细节有差异, 但大体图案仍然相似; 然而对于不同破裂模型计算的库仑应力分布, 虽然也具有相似的模式, 即在破裂面上和破裂面两端外延的方向库仑应力增加, 而在断层面两侧区域库仑应力下降, 但在具体图形上、特别在接近断层的近场的分布, 存在相当显著的差别.

对于逆掩断层, 地震的发生释放了垂直于断层走向的水平压应力, 因此断层两侧再发生类似逆掩断层的危险性大大减小, 也就是计算中各个模型都显示出的断层两侧(特别在断层逆掩成分为主的中、南段)库仑应力减小区域.但在断层上盘由于垂直运动和水平运动综合作用造成的复杂变形, 库仑应力计算表明地震危险性反而增加, 这就是3个模型中沿地表断层西北一侧库仑应力增加的窄长条.吴建平等[33]给出的余震分布剖面表明, 在逆掩为主的中南段, 余震几乎都发生在断层的上盘地壳中, 与库仑应力计算结果吻合.断层两端一般是应力集中的地方, 容易进一步发生地震, 3个模型的端部均为应力增加区.

模型二断层破裂长度接近模型一的两倍, 是3个模型中最大的, 而滑动量是3个模型中最小的.模型二将整个破裂面沿走向划分成31段, 沿倾向划分成5段, 形成155块子断层, 各子断层有不同的滑动量和相同的滑动角, 根据模型二计算的库仑破裂应力分布显示因活动断层较长, 因此两侧库仑应力下降区域也较大(图 3c, 图 3d).

模型三与模型一相比, 断层破裂向东北方向延伸较长, 东北部的库仑应力增加区也较大(图 3e, 图 3f), 与该段余震活动分布较为吻合.

4 讨论

本文改进了库仑应力计算方法, 根据汶川大地震3个地震破裂模型计算了静态库仑应力变化.通常研究者们计算地震库仑应力变化, 目的之一是考察主震对余震的触发作用, 往往统计发生在库仑应力增加区和库仑应力减小区内余震数目和比率, 如果大多数余震发生在库仑应力增加区, 常常被认为是支持库仑应力增加对触发余震的证据.我们这里也对发生在库仑应力增加区域内的余震(由于后续强余震陆续发生会改变应力分布情况, 因此仅限于讨论主震后一个月内面波震级大于4级的余震)作了一个简单统计, 发现3个模型结果存在一些差别.例如, 有将近80%的余震发生在用模型一计算的库仑应力增加区域内, 而只有约60%的余震发生在模型二和模型三计算的库仑应力增加区域.考虑到实际应用中人们可能希望对余震范围作出估计, 从保险一些、尽量不遗漏余震来讲, 我们将各点在3个模型得到的库仑应力最大值抽取出来作为综合模型, 称之为模型四, 也参与了统计, 96%以上的余震发生在合成的模型四的库仑应力增加区内.图 4是各个模型比较的结果.由于地球物理反演问题复杂性, 不同研究者反演得到的地震破裂模型之间的差异是难以避免的.库仑应力计算过程中摩擦系数选取会影响到计算的库仑应力分布, 使同一个模型对余震分布拟合程度产生差异, 通过考察库仑应力来预测余震分布还存在一些问题.然而, 引人注意的是, 的确均能有60%~80%的余震发生在这3个震后数天不同研究者利用有限资料做出的破裂模型的库仑应力增加区内; 而且把它们的库仑应力增加区合成为一体的模型四, 余震发生在库仑应力增加区的比率竟高达96%以上, 这也提示了, 在大地震后综合各家利用有限资料做出的破裂模型和相应的库仑应力分布图, 有可能对余震发生区域做出较好的判断.

图 4 发生在不同地震破裂模型库仑应力增加区域的余震占余震总数的比例 Fig. 4 The ratio of the number of aftershocks occurred where Coulomb stress increases to the number of aftershocks in one month after the main shock

虽然库仑破裂应力的计算容易受到一些因素影响, 但在大地震发生后, 人们倾向于准备应付最危险的情况, 会对潜在地震危险性做出最保守估计, 将不同地震破裂模型综合起来考察库仑应力变化, 可以避免个别模型可能忽视掉的问题, 有利于人们在大震之后判断后续地震的趋势.我们将3个模型得到的库仑应力综合后提取最大值, 观察其库仑应力分布(图 5).虽然同在库仑应力增加区内, 但不同地点的库仑应力增加值相差甚多, 可达几个数量级, 真正能达到MPa水平的应力变化区域并不大.为了真实反映库仑应力增加区内不同的变化量级, 在图 5中绘的是库仑应力增加值对数的等值线图.

图 5 在0 km深度的库仑应力正值的常用对数值 (a)对应等效摩擦系数取0.85的情况;(b)对应等效摩擦系数取0.4的情况. Fig. 5 Common logarithm of Coulomb stress increase at zero depth (a) Given friction coefficient as 0.85; (b) Given friction coefficient as 0.4.

图 5可以看出大多数余震都发生在综合模型库仑应力增加且增加值较大的区域, 尽管如此, 还应该注意, 现在计算的库仑应力是假设存在与主震破裂面平行的潜在断层面上的库仑应力分布.倘若余震发生的断裂面与主震破裂面具有不同的走向和倾角, 还要对这些细节作深入研究, 才可以对主震触发余震问题得到更深层次的认识.

另外, 我们计算库仑应力时采用的构造应力取值源自近地应力测量结果, 只适合用来讨论近地表的库仑应力变化, 地下几公里或者十几公里的应力大小如何, 虽然国际上个别钻孔(例如德国大陆深钻到9 km)取得了一些认识, 但在龙门山断裂带我们还缺乏类似资料, 所以应用库仑应力讨论地震的触发作用也有待深入研究.

5 结论

由于人们往往不知道区域应力场的大小和方向, 因此试图计算库仑应力变化量, 回避初始应力场未知的问题, 讨论地震同震变形和应力变化对后续地震的影响.在这样的计算中, 假定震前和震后地应力方向没有变化, 发震断层错动方向上应力分量的变化就是剪应力的变化.然而, 本文考虑震前后主应力方向的变化, 修正了传统的库仑应力计算中仅沿地震破裂面滑动方向计算剪应力变化的近似方法, 考虑震后主应力方向改变对剪应力变化量计算的影响, 通过对比改进方法和传统方法获得的结果, 认识到虽然在远场研究中原有方法能给出很好近似, 但在邻近断层端点数公里的小区域内, 二者库仑应力变化有显著差异, 这种关键部位的应力差异对研究断层的错动、发展和终止会有不可忽视的影响, 值得在今后进一步研究.

本文还基于龙门山地区背景应力场测量资料, 根据不同研究者利用震后资料数天内反演的地震破裂模型, 计算了汶川地震造成的静态库仑应力变化, 考察了不同地震破裂模型、不同孔隙水压或等效摩擦系数下, 计算出的库仑应力分布的差异, 以及利用各自的库仑应力判断大震后余震发展趋势是否能得到一致结果的问题.认识到孔隙水压虽然有一定影响, 不会对库仑应力整体分布图像有太大影响; 但不同破裂模型会使库仑应力分布产生显著的差异.不过, 各个模型在不同程度上都显示了库仑应力增加较高值的区域, 与余震发生有较好的对应关系.建议在利用多个地震破裂模型计算库仑应力和探讨未来余震区域时, 从增大安全系数出发, 可以建立多个模型的综合模型, 对余震危险区域进行估计.

致谢

感谢中国科学院地质与地球物理研究所姚振兴院士和中国地震局地球物理研究所陈运泰院士为本研究提供数据.

参考文献
[1] Stein R S, Lisowski M. The 1979 Homestead Valley earthquake sequence, California:control of aftershocks and postseismic deformation. J. Geophys. Res. , 1983, 88(B8): 6477-6490. DOI:10.1029/JB088iB08p06477
[2] King G C P, Stein R S, Lin J. Static stress changes and the triggering of earthquakes. Bull. Seismol. Soc. Amer. , 1994, 84(3): 935-953.
[3] Stein R S, King G C, Lin J. Change in failure stress on the southern San Andreas fault system caused by the 1992 magnitude=7.4 Landers earthquake. Science , 1992, 258(5086): 1328-1332. DOI:10.1126/science.258.5086.1328
[4] Stein R S, King G C P, Lin J. Stress triggering of the 1994 M=6.7 Northridge, California, earthquake by its predecessors. Science , 1994, 265(5177): 1432-1435. DOI:10.1126/science.265.5177.1432
[5] Cianetti S, Giunchi C, Cocco M. Three-dimensional finite element modeling of stress interaction:an application to Landers and Hector Mine fault systems. J. Geophys. Res. , 2005, 110(B05S17). DOI:10.1029/2004JB003384
[6] Freed A M. Earthquake triggering by static, dynamic, and postseismic stress transfer. Annu. Rev. Earth. Pl. Sc. , 2005, 33: 335-367. DOI:10.1146/annurev.earth.33.092203.122505
[7] Harris R A. Introduction to special section:stress triggers, stress shadows, and implications for seismic hazard. J. Geophys. Res. , 1998, 103(B10): 24347-24358. DOI:10.1029/98JB01576
[8] Freed A M, Lin J. Delayed triggering of the 1999 Hector Mine earthquake by viscoelastic stress transfer. Nature , 2001, 411(6834): 180-183. DOI:10.1038/35075548
[9] 石耀霖. 关于应力触发和应力影概念在地震预报中应用的一些思考. 地震 , 2001, 21(3): 1–7. Shi Y L. Stress triggers and stress shadows:how to apply these concepts to earthquake prediction. Earthquake (in Chinese) , 2001, 21(3): 1-7.
[10] 万永革, 吴忠良, 周公威, 等. 地震应力触发研究. 地震学报 , 2002, 24(5): 533–551. Wan Y G, Wu Z L, Zhou G W, et al. Research on seismic stress triggering. Acta Seismologica Sinica (in Chinese) , 2002, 24(5): 533-551.
[11] Stein R S, Barka A A, Dieterich J H. Progressive failure on the North Anatolian fault since 1939 by earthquake stress triggering. Geophys. J. Int. , 1997, 128(3): 594-604. DOI:10.1111/gji.1997.128.issue-3
[12] 沈正康, 万永革, 甘卫军, 等. 东昆仑活动断裂带大地震之间的黏弹性应力触发研究. 地球物理学报 , 2003, 46(6): 786–795. Shen Z K, Wan Y G, Gan W J, et al. Viscoelastic triggering among large earthquakes along the east Kunlun fault system. Chinese J. Geophys. (in Chinese) , 2003, 46(6): 786-795.
[13] 郝平, 傅征祥, 田勤俭, 等. 昆仑山口西8.1级地震强余震库仑破裂应力触发研究. 地震学报 , 2004, 26(1): 30–37. Hao P, Fu Z X, Tian Q J, et al. Large aftershocks triggering by Coulomb failure stress following the 2001 Ms=8.1 great Kunlun earthquake. Acta Seismologica Sinica (in Chinese) , 2004, 26(1): 30-37.
[14] 万永革, 沈正康, 曾跃华, 等. 青藏高原东北部的库仑应力积累演化对大地震发生的影响. 地震学报 , 2007, 29(2): 115–129. Wan Y G, Shen Z K, Zeng Y H, et al. Evolution of cumulative Coulomb failure stress in northeastern Qinghai-Xizang (Tibetan) Plateau and its effect on large earthquake occurrence. Acta Seismologica Sinica (in Chinese) , 2007, 29(2): 115-129.
[15] 万永革, 沈正康, 曾跃华, 等. 唐山地震序列应力触发的黏弹性力学模型研究. 地震学报 , 2008, 30(6): 581–593. Wan Y G, Shen Z K, Zeng Y H, et al. Study on visco-elastic stress triggering model of 1976 Tangshan earthquake. Acta Seismologica Sinica (in Chinese) , 2008, 30(6): 581-593.
[16] Parsons T, Ji C, Kirby E. Stress changes from the 2008 Wenchuan earthquake and increased hazard in the Sichuan basin. Nature , 2008. DOI:10.1038/nature07177
[17] Toda S, Lin J, Meghraoui M, et al. 12 May 2008M=7.9 Wenchuan, China, earthquake calculated to increase failure stress and seismicity rate on three major fault systems. Geophys. Res. Lett. , 2008, 35(L17305). DOI:10.1029/2008GL034903
[18] 万永革, 沈正康, 盛书中, 等. 2008年汶川大地震对周围断层的影响. 地震学报 , 2009, 31(2): 128–139. Wan Y G, Shen Z K, Sheng S Z, et al. The influence of 2008 Wenchuan earthquake on surrounding faults. Acta Seismologica Sinica (in Chinese) , 2009, 31(2): 128-139.
[19] 陈兵, 江在森, 车时, 等. 玛尼7.9级地震对昆仑山口西8.1级地震的触发作用及动力背景初探. 中国地震 , 2003, 19(1): 1–7. Chen B, Jiang Z S, Che S, et al. Study of triggering action between Mani (Ms7.9) and Kunlun (Ms8.1) great earthquake and their dynamic background. Earthquake Res.China (in Chinese) , 2003, 19(1): 1-7.
[20] 刘桂萍, 傅征祥. 1976年7月28日唐山7.8级地震触发的区域地震活动和静应力场变化. 地震学报 , 2000, 22(1): 17–26. Liu G P, Fu Z X. Regional seismicity triggered by the Ms7.8 Tangshan event of July 28, 1976 and the static stress field change. Acta Seismologica Sinica (in Chinese) , 2000, 22(1): 17-26.
[21] 张秋文, 张培震, 王乘, 等. 鲜水河断裂带断层间相互作用的触震与缓震效应. 地震学报 , 2003, 25(2): 143–153. Zhang Q W, Zhang P Z, Wang C, et al. Earthquake triggering and delaying caused by fault interaction on Xianshuihe fault belt, southeastern China. Acta Seismologica Sinica (in Chinese) , 2003, 25(2): 143-153.
[22] 周仕勇. 川西及邻近地区地震活动性模拟和断层间相互作用研究. 地球物理学报 , 2008, 51(1): 165–174. Zhou S Y. Seismicity simulation in western Sichuan of China based on the fault interactions and its implication on the estimation of the regional earthquake risk. Chinese J. Geophys. (in Chinese) , 2008, 51(1): 165-174.
[23] Hardebeck J L, Nazareth J J, Hauksson E. The static stress change triggering model:constraints from two southern California aftershocks equences. J. Geophys. Res. , 1998, 103(B10): 24427-24437. DOI:10.1029/98JB00573
[24] Zhao D, Kanamori H, Wiens D. State of stress before and after the 1994 Northridge earthquake. Geophys. Res. Lett. , 1997, 24(5): 519-522. DOI:10.1029/97GL00258
[25] Okada Y. Surface deformation due to shear and tensile faults in a half-space. Bull. Seismol. Soc. Amer. , 1985, 75: 1135-1154.
[26] Okada Y. Internal deformation due to shear and tensile faults in a half-space. Bull. Seismol. Soc. Amer. , 1992, 82(2): 1018-1040.
[27] Gomberg J, Ellis M. 3D-DEF:a User′s Manual, in Open-File Report. Denver , 1993: 93-547.
[28] Gomberg J, Ellis M. Topography and tectonics of the central New Madrid seismic zone:results of numerical experiments using a three-dimensional boundary element program. J. Geophys. Res. , 1994, 99(B10): 20299-20310. DOI:10.1029/94JB00039
[29] 郭啟良, 王成虎, 马洪生, 等. 汶川Ms8.0级大震前后的水压致裂原地应力测量. 地球物理学报 , 2009, 52(5): 1359–1401. Guo Q L, Wang C H, Ma H S, et al. In-situ hydro-fracture stress measurement before and after the Wenchuan Ms8.0 earthquake of China. Chinese J. Geophys. (in Chinese) , 2009, 52(5): 1359-1401.
[30] 周荣军, 李勇, DensmoreA L, 等. 青藏高原东缘活动构造. 矿物岩石 , 2006, 26(2): 40–51. Zhou R J, Li Y, Densmore A L, et al. Active tectonics of the eastern margin of the Tibet Plateau. J. Mineralogy Petrology (in Chinese) , 2006, 26(2): 40-51.
[31] 王卫民, 赵连锋, 李娟, 等. 四川汶川8.0级地震震源过程. 地球物理学报 , 2008, 51(5): 1403–1410. Wang W M, Zhao L F, Li J, et al. Rupture process of the Ms8.0 Wenchuan earthquake of Sichuan. China.Chinese J. Geophys. (in Chinese) , 2008, 51(5): 1403-1410.
[32] 安美其, 丁立丰, 王海忠, 等. 龙门山断裂带的性质和活动性研究. 大地测量与地球动力学 , 2004, 24(2): 115–119. An M Q, Ding L F, Wang H Z, et al. Research of property and activity of Longmen mountain fault zone. J. Geodesy Geodynamics (in Chinese) , 2004, 24(2): 115-119.
[33] 吴建平, 黄媛, 张天中, 等. 汶川Ms8.0级地震余震分布及周边区域P波三维速度结构研究. 地球物理学报 , 2009, 52(2): 320–328. Wu J P, Huang Y, Zhang T Z, et al. Aftershock distribution of the Ms8.0 Wenchuan earthquake and three dimensional P wave velocity structure in and around source region. Chinese J.Geophys. (in Chinese) , 2009, 52(2): 320-328.