紫坪铺库区小震产生的库仑应力变化及其与汶川地震的关系

尹迪,董培育,石耀霖

1 中国科学院计算地球动力学重点实验室，中国科学院大学，北京 100049 2 中国地震局地震研究所，地震大地测量重点研究室，武汉 430071

0 引言

2008年5月12日，四川省汶川县发生了MS8.0大地震，引起了国内外地震学家的广泛关注.随着对龙门山断裂带成因、活动性及构造演化的深入研究，区域构造运动和地球深部动力学过程被认为是造成汶川地震发生的主要原因(邓启东等，1994;张培震等，2003,2008;郑勇等，2009;李平恩等，2019).

然而，距离汶川震源不足10 km处为紫坪铺水库，其库尾位于汶川地震发震构造带上，因此紫萍铺水库蓄水是否对汶川地震具有触发作用成为了国内外学术界广泛争议的热点.争议主要存在两种观点，第一种观点是紫坪铺水库具备诱发地震的基本条件(Kerr and Stone,2009;Ge et al.,2009；Lei,2011)，其通过建立有限元二维均匀介质数值模型，模拟紫坪铺水库水体静态荷载和孔隙压力扩散情况下的库仑应力变化，获得蓄水后的库仑应力增量达到了地震触发阈值，因此得出紫坪铺水库促使汶川地震提前 10～100 年发生.相反的观点则认为紫坪铺蓄水造成的库仑应力变化低于地震诱发阈值，尚未达到触发汶川地震的条件(Deng et al.,2010;Gahalaut and Gahalaut,2010；孙玉军等，2012).Shi等(2013)和程慧红等(2015)研究发现，由于不同研究者采用的参数不同，如汶川地震震源深度、震源机制、断层和岩体渗透率等不同，使用的方法不同，如二维模型还是三维模型、解析解还是数值解等，均会导致模拟结果存在较大差异，甚至得到两种截然不同的结论.因此，回答紫坪铺水库是否对汶川地震具有触发作用需要更加深入的野外调查和物理数值实验研究，也需要多方位多角度的讨论.

然而紫坪铺水库建成蓄水之后，其周边区域内发生不少小震，而且在时间和空间上存在成簇的特征，有研究认为它们是水库蓄水诱发的地震(程万正等，2010)(图1).关于这一系列诱发小震及其动态响应是否对汶川地震具有促进作用也引起了人们的思考(Ruan et al.,2017；刘远征等，2014)，但尚未有相关的定量研究.因此，为了定量分析库区小震是否影响了汶川地震的发生，我们搜集了紫坪铺库区ML震级在2.1～3.9的65次地震事件的震源机制解(胡先明等，2009；李春政等，2018)，将小震滑动破裂作为扰动源，采用Wang 等(2006)开发的PSGRN/PSCMP计算程序，计算这些小震事件产生的同震静态应力扰动，根据库仑破裂应力变化理论，获得紫坪铺水库诱发地震在汶川地震震源处的应力扰动量，同时考虑了小震造成的动态应力改变及振动对岩层渗透率、强度的影响，综合分析了其对汶川地震发生的作用.

1 地震观测资料

1.1 汶川地震震前紫坪铺库区小地震事件

已有的研究资料表明，大中型水库蓄水会对库区周围造成一定的应力影响，使得库区周围小型地震频次增加(程惠红等，2015)，甚至会触发一些较大的地震(Stein et al.,1982;Magadza,2006;Durá-Gómez and Talwani,2010).据张致伟等(2009)，程万正等(2010)对紫坪铺水库地震台网和地震台记录的地震资料的分析，认为水库库区蓄水后的小震活动可以视为诱发小震.胡先明等(2009)利用PS垂直分量振幅比资料测定了紫坪铺库区2004年8月至2007年12月262次小震震源机制参数；张永久和张致伟(2010)运用PS振幅比法计算了紫坪铺库区及周边2004年8月至2008年5月486次震级大于1.6级的震源机制参数.马文涛等(2011)研究表明小震位于紫坪铺库区十几公里范围内，从紫坪铺水库蓄水(2005年9月30日)到汶川地震发生前，地震震级虽未出现明显增大，但地震次数明显增多，且小震在空间上呈现出丛集分布的特征，时间集中于几个不同的时间段，分别为2006年1月的水磨震群、2006年3月深溪沟震群及2008年2月至4月的都江堰震群等(胡先明等，2009；卢显等，2010).

紫坪铺水库台网记录到的2004年8月16日至2008年5月10日的小震有1569个(卢显等，2010)，到汶川地震发生前，库区内发生的最大地震为2008年2月的都江堰3.7级地震.由于区域小震数目较多，且根据地震波能量与震级关系可知，震级相差1级，释放的能量相差约30倍，小震释放的能量有限(Gutenberg,1956).因此本文仅选取紫坪铺库区2004—2008年ML>2.1级的65次小震用以计算分析，见图1.

计算紫坪铺小震对汶川震源处产生的应力扰动，需将这些小震的破裂滑动分布作为扰动源,每个小震的破裂参数可根据龙锋等(2006)给出的面波震级MS与破裂长度(L)及面积(A)的经验公式(如式(1)所示)计算得到.而本文中的地震震级为近震震级ML，根据我国常规震级标度与矩震级之间的对应关系(李莹甄等，2014)，小震震级ML<4.0级，震级转换关系可视为ML≈MS≈MW.因此根据公式(1)，可得到每个小震的破裂长度(L)、破裂宽度(W)以及滑动量(slip)，其中G为剪切模量，参照前人对龙门山地区的研究成果(易桂喜等，2013)，取G=3.0×1010Pa，各个小震的破裂模型见附表1.

MS=3.821+1.860 lg(L),

MS=4.134+0.954 lg(A),

MO=G×slip×A.

(1)

1.2 汶川地震震源参数

汶川地震发生后，不同的研究机构和学者根据不同的数据和方法，获取了汶川地震的震源参数信息.中国地震台网中心根据93个地震台在汶川震后一个月的记录资料，将汶川地震的震源参数修订为：震中30.95°N,130.57°E,震源深度14 km；陈九辉等(2009)综合利用川西流动台站观测数据和震后应急地震观测台站的震相数据，采用双差定位方法确定的震中位置与台网中心结果基本一致，但震源深度测定为18.8 km；王卫民等(2008)结合远场体波波形记录及近场同震位移数据，利用反演技术重建地震破裂过程，得到汶川地震震源深度为15.4 km，并给出震源机制解结果：走向229°、倾角32°、滑动角118°；杨智娴等(2012)综合运用四川省地震台网的观测资料以及紫坪铺水库地震台网的记录资料反演了汶川地震更精准的震源位置，指出震中位置为31.018°N,103.365°E，震源深度为15.5 km.大部分学者普遍认为汶川震源深度位于14 km以下，但Gong等(2019)通过对紫坪铺库区近源地震记录到的P波相位进行反演，指出汶川地震初始破裂位置为31.013°N，103.391°E，震源深度较浅，仅为8.4 km，详见表1.本文综合考虑前人的研究成果，分别计算在这些不同参数条件下，紫坪铺小震在汶川震源处造成的库仑应力变化量.

2 方法与模型

2.1 计算方法

地震在释放自身能量的同时会对周边区域应力场产生扰动，将该应力变化张量投影到接收断层面上，可定量计算对该断层面是否会产生促进或抑制滑动的作用.应力扰动张量通过Wang等(2006)开发的PSGRN/PSCMP计算程序，将地震破裂模型作为输入源计算得到，由于库区内地震震级较小，可以忽略震后的黏弹性松弛效应，仅考虑同震应力变化，随后采用库仑破裂应力变化理论计算库仑应力变化，如公式(2)所示：

ΔCFS=Δτ+μ(Δσn+ΔP).

(2)

采用弹性力学定义，张应力为正值，其中Δτ为断层面上剪应力变化量(沿着断层面滑动方向为正)，Δσn断层面上正应力变化量，ΔP为孔隙压力变化，μ为断层摩擦系数(King et al.,1994;Cocco and Rice,2002).若断层面上剪应力增加或正应力增加，那么ΔCFS>0，断层更易滑动，反之则断层面更加稳定.由于无法知道孔隙压力变化量，常利用有效摩擦系数μ′来反映孔隙压力变化，即(2)式变为(3)式，研究中通常取μ′=0.2-0.6,本文取值0.4.

ΔCFS=Δτ+μ′Δσn.

(3)

2.2 地壳分层结构

龙门山断裂带西侧为活跃的巴颜喀拉块体，东侧为古老而稳定的四川盆地，两侧地质构造和地壳速度结构均存在明显差异(邓起东等，1994;Burchfiel et al.,1995)，本文重点讨论小震对汶川震源处的影响，参考前人对汶川地区地震波速的研究成果(雷建设等，2009；杨智娴等，2012)确定模型的速度结构(见表2).

表2 本文采用的地壳速度结构模型Table 2 Velocity structure model used in this study revised

3 计算结果

3.1 典型小震的计算结果

由于小震数量较多，不一一呈现每个地震的计算结果.仅以震级较大且距离汶川震中较近的两次地震为例，分别为2005年7月20日庙子坪3.6级地震及2008年2月14日都江堰3.7级地震(见图1)，其震源深度分别为：10 km和7 km，具体震源破裂参数见表3，分析其在汶川地震震源处(31.018°N,103.365°E,深度15.5 km)产生的应力张量变化，如图2所示.

表3 两次小震的震源破裂参数Table 3 Source rupture parameters of two small earthquakes

图2 两次典型小震在汶川地震震源(15.5 km)水平截面上造成的应力扰动其中红色五角星为汶川地震位置，黄色圆圈为小震发震位置.(a)—(c)为2005年7月20日庙子坪3.6级地震产生的应力扰动量ΔSxx，ΔSyy，ΔSxy；(d)—(f)为2008年2月14日都江堰3.7级地震产生的应力扰动量ΔSxx，ΔSyy，ΔSxy.Fig.2 The stress change on the horizontal section of the Wenchuan earthquake hypocenter (15.5 km)caused by the two typical small earthquakesThe red star and yellow circle represent the location of Wenchuan earthquake and small earthquake respectively.(a)—(c)ΔSxx,ΔSyy,ΔSxy of Miaoziping M 3.6 earthquake on July 20,2005 ；(d)—(f)ΔSxx,ΔSyy,ΔSxy of Dujiangyan M 3.7 earthquake on February 14,2008.

由图2可知，小震的影响范围仅在震源数公里内且很快衰减，两次小震在汶川震源处造成的应力扰动ΔSxx分量分别为1.22 Pa和1.80 Pa；ΔSyy分量分别为0.87 Pa和-1.64 Pa;ΔSxy分量分别为1.11 Pa和1.09 Pa.将两次小震产生的应力张量变化投影到汶川地震破裂面上，走向/倾向/滑动角分别为：229°/32° /118°(王卫民等，2008)，得到库仑应力变化(表4),分别仅为0.75 Pa和-1.78 Pa，量值很小，影响微弱.

表4 两次典型小震在汶川震源处产生的应力扰动Table 4 Stress changes caused by two typical small earthquakes to the Wenchuan earthquake source

另外最近有学者指出，汶川地震的破裂可能始于8.4 km深度处(Gong et al.,2019),为了考察接收断层不同震源深度对结果的影响，本文过汶川震中做了一条长为50 km的剖面(图1蓝线所示)，计算了两次小震沿着剖面0～16 km深度的应力张量变化，并投影到王卫民等(2008)的断层接收面参数下，得到了库仑应力随深度的变化图(图3)，两次小震在8.4 km深度处造成的库仑应力变化量的值分别为2.53 Pa和2.22 Pa，与15.5 km深度处的结果(0.75 Pa和-1.78 Pa)虽然存在明显差异甚至发生正负的改变，但库仑应力变化绝对值的影响仍非常微弱.

图3 两次典型小震在汶川震源剖面上产生的ΔCFS(a)2005年7月20日庙子坪3.6级地震；(b)2008年2月14日都江堰3.7级地震.Fig.3 ΔCFS on the profile of Wenchuan hypocenter caused by the two typical small earthquakes(a)Miaoziping M 3.6 earthquake occurred on July 20,2005;(b)Dujiangyan M 3.7 earthquake on February 14,2008.

3.2 全部小震的计算结果

依次计算了紫坪铺库区65次小震在汶川震源处造成的应力张量变化，经统计，有34次小震在汶川震源处产生的同震库仑应力变化为正值，其中影响最大的是第12次小震，ΔCFS约为3.2 Pa；有31次小震产生的ΔCFS为负值，影响最大的是第60次小震，ΔCFS约为-1.78 Pa.将全部小震产生的应力张量变化累加求和，ΔSxx、ΔSyy和ΔSxy变化总量分别为：-2.56 Pa、5.80 Pa、11.37 Pa，计算库仑应力变化量为4.97 Pa，与单次小震产生的影响量级相当.表明这些小震对汶川震源处产生的影响十分有限.

4 讨论

4.1 不同参数的影响

库仑应力变化的计算与接收断层面的几何参数及有效摩擦系数密切相关.本文选取不同机构给出的汶川地震震源机制解(表1)作为接收断层面，计算震源处深度15.5 km结果，如表5所示，除了USGS参数下计算得到的库仑应力增量为负，其余情况下得到的库仑应力增量均为正值.尽管略有差异，但由于小震本身产生的应力扰动量(绝对值)较小，因此在不同震源机制解下，小震对汶川震源处造成的库仑应力变化值的影响均很微弱，不足以触发断层滑动.值得注意的是，以王卫民等(2008)的震源机制解下的结果为参照，不同震源机制解下的应力变化误差最高可达90%；另外若计算深度为8.4 km，全部小震造成的库仑应力变化为-4.0 Pa，与同样计算条件下15.5 km处的结果正负相反(表5).由此可见，虽然小震产生的应力扰动绝对值很小但库仑应力受震源机制解及计算深度影响较大.

表5 全部小震汶川震源处造成的库仑应力变化Table 5 Coulomb stress changes at the Wenchuan source of all small earthquakes

另外，关于有效摩擦系数对计算结果的影响，已有研究结果表明虽然会影响到计算结果，但量级和正负方向上不会有变化(董培育等，2020；黄禄渊等，2019).本文计算结果仅为几个Pa左右，远未达到触发阈值10 kPa的量级，如公式(3)所示，有效摩擦系数的改变不会造成结果的改变，因此忽略有效摩擦系数对计算结果的影响.

4.2 小震的近距离影响

虽然小震产生的应力扰动随距离衰减很快，但是在数公里内仍有较显著的影响.本文选取了2005年7月20日的3.6级地震和2005年9月10日的2.6级地震(附表1中第15和20次小震)，震源深度均为10 km.计算了距各自震源3 km范围内产生的同震应力变化，由图4所示，可以看出3.6级地震在3 km范围内应力变化可高达50 kPa，2.6级地震在1 km范围内应力也可高达50 kPa，大于一般认为的10 kPa 的触发地震所需要的阈值(King et al.,1994).图2中的小震震源深度分别是10 km和7 km，而计算深度在15.5 km处，在小震震中周围的同震应力变化仅为几个Pa的量级，说明随着深度的衰减较快于随水平面的衰减.考虑小震间近距离(3～5 km尺度内)产生的应力扰动可能会触发地震，例如2019年四川长宁地震前，许多小震发生在震中10 km的范围内，这些小震的相互影响也是今后值得注意的问题(易桂喜等，2019).

附表1 小震震源破裂模型Appendix Table 1 Small earthquake source rupture model

续附表1

图4 两次典型小震在小震震源10 km水平截面上造成的应力扰动(3 km范围)(a)—(c)2005年7月20日3.6级地震产生的应力扰动量ΔSxx，ΔSyy，ΔSxy；(d)—(f)2005年9月10日2.6级地震产生的应力扰动量ΔSxx，ΔSyy，ΔSxy.Fig.4 The stress change on the horizontal section of the earthquake hypocenter (10.0 km)caused by the two typical small earthquakes (3 km range)(a)—(c)ΔSxx，ΔSyy，ΔSxy of M 3.6 earthquake on July 20,2005；(d)—(f)ΔSxx,ΔSyy,ΔSxy of M 2.6 earthquake on September 10,2005.

5 讨论与结论

本文计算了汶川地震发生前，紫坪铺库区的65个小震在汶川震源处产生的静态应力扰动量.计算结果主要受到投影面即汶川地震震源机制解参数的影响.另外在相同条件下，不同深度面上的计算结果也可能出现反向变化，如8.4 km和15 km深度处的库仑应力正负不同.但整体上小震在汶川震源处产生的同震应力扰动总量仅为几个Pa的量级.前人在龙门山地区开展的数值模拟结果表明，该区域构造应力年增长速率约为1～2 kPa(柳畅等，2012；万永魁等，2017)，与构造应力积累速率相比，小震产生的静态应力扰动基本可忽略，认为其对汶川地震没有直接触发作用.

然而地震除了产生静态应力扰动之外，其地震波的振动还会引起动态应力的变化，尽管不会像静态应力变化在近场触发较大震级的地震，但也可能会触发较远距离的地震(石耀霖等，2021).现有的研究观测并证实了动态应力触发现象的存在(Gomberg et al.,1998;Kilb et al.,2002)，例如1995年亚喀巴海湾MS7.3地震发生2.8 h后，距离地震震中500 km的叙利亚黎巴嫩地区发生20多次的地震群，其中最大震级3.5级，有学者认为是受到动态应力的触发作用(Mohamad et al.,2000)；2015年尼泊尔MS8.1大地震远距离触发了中国西藏、重庆地区的一系列地震活动，造成25天后区域地震活动再次明显(万永革等,2015;张贝等,2015).我国西北和东北境内的两个四分量钻孔应变仪，记录到2018年发生在相距数千公里外的四次大震(斐济群岛M8.1、印尼M7.4、巴布亚新几内亚M7.1和日本千岛群岛M6.6)产生的动态应力变化，峰值可达数百帕(李富珍等，2021).2016年山西原平地区发生了ML4.7地震，距震源19 km处的原平台站钻孔应力仪记录到的动态应力峰值在百帕量级(石耀霖等，2021).王琼等(2016)计算2014年于田MS7.3地震对30 km外区域产生的动态应力变化峰值约3.0 MPa,而静态应力仅为0.8 MPa.动态应力峰值随着距离r-1衰减(Steacy et al.,2005)，而线弹性理论计算结果显示地震产生的静态应力变化随着距离r-3衰减(Hill et al.,1993)，衰减速度远快于动态应力.当震中距大于1～2个断层长度时，可认为动态应力峰值较静态应力峰值大，起主导作用，因此远距离的地震触发通常认为是动态应力的作用.紫坪铺水库所诱发的一系列震级在2.2～4之间的小震，其断层长度多数低于0.1 km，最长仅约为0.5 km,而距离汶川震中为10～20 km范围，参考山西原平台的记录数据，可能同样会对汶川地震处造成几十到上百帕范围内的动态应力扰动，尽管动态应力峰值明显大于静态应力变化，但量值仍然远低于构造应力积累年速率.尽管如此，本文小震造成应力变化的量级也不足以触发汶川地震.

实验室开展的岩石力学实验表明，随时间变化的波动可以通过改变孔隙岩石颗粒及破坏胶结物的聚集，从而引发渗透率及岩石内部结构的改变(郑茂盛等，2008；朱立等，2012).有观测现象表明地震波引起的低至10-4Jm-3的能量密度仍然可以引发渗透率的变化，如2002年的M7.9的Denali地震造成的距离其5000 km远的Iowa的地下水渗流的变化就发生在10-4Jm-3的能量密度量级下(Wang and Manga,2010,2014).地质观测也发现了地震波振动可能会改变岩体或断裂带的力学性质从而触发地震活动(Moyer et al.,2018)，如Taira等(2009)利用美国圣安德烈亚斯断层20年(1987—2008年)的观测资料，分析断层强度的变化，发现1992年兰德斯地震可能降低了1993年帕克菲尔德地震震源区的岩石强度，进而触发了此次地震，甚至2004年苏门答腊M9.1地震都对圣安德烈斯断层强度有重要影响.Zhang等(2017)通过观测我国东北地区某一口井水位变化，发现三次大震(2008汶川MW7.9地震，2011日本东北MW9.1地震以及2012苏门答腊MW8.6地震)后均出现水位阶变现象，且与体应变变化呈正相关关系，分析认为大地震可能导致远场介质的渗透率发生变化，进而引起孔隙压变化，导致同震体应变远大于理论计算值.那么小震是否也会引起局部区域断层强度或岩体渗透率变化呢？石耀霖等(2021)利用山西原平钻孔应力仪观测到的相距19 km处发生的ML4.7地震数据，认为小震振动也可能引起地下水位的变化，进而造成孔隙水压变化，推测这可能是观测数据与理论值存在偏差的原因之一.尽管小震的震动持续时间较短，仅在几秒到十几秒钟，且引起的振动较小，也可能会造成周边区域介质的渗透性发生改变，进一步加速了岩层中水的渗透，从而增加孔隙压力.

本文定量计算结果表明小震群在汶川震源处产生的静态应力扰动与汶川地震没有直接关系，但动态应力扰动量级比静态应力大1～2个数量级，且小震在近距离(3～5 km)内产生的静态应力扰动量级较大(数十至百千帕)，是否小震之间相互作用，连续振动影响了临近区域的介质性质，进而与汶川地震的发生具有间接关系？这些问题需要进一步深入定量计算研究.