2008年5月12日在中国四川汶川发生8.0级强烈地震，造成重大人员伤亡和破坏.从动力学背景来看，在印度板块与欧亚大陆板块碰撞挤压作用下，喜马拉雅造山带东构造结向NNE 方向顶挤、楔入青藏高原东北缘，迫使高原深部物质向东流，在龙门山地区受到四川盆地阻挡而发生汶川地震^[1].在大的构造背景基础上，人们往往也在寻找直接的诱发因素，如在《科学》(Science)杂志上就发表过一篇报道，介绍了“紫坪铺水库可能诱发了四川汶川地震" 的说法^[2-3].水库蓄水放水能够诱发地震已经被人们所证实.而这次汶川大地震的微观震中离库区最近距离只有约6km,库区跨越汶川大地震的主要发震断层之一的映秀-北川断层的破碎带^[4].因此，是否是紫坪铺水库诱发了汶川大地震在科学界也引起了不少争论.

雷兴林根据紫坪铺水库的蓄放水过程，考虑了水库载荷的弹性效应和水的渗透效应，计算得到在震源深度处库仑应力能够达到1bar(10⁵Pa)的量级，因此他认为紫坪铺水库对龙门山断层有明显的作用^[4-5].陈颙通过分析水库地震的特点和力学上的解释，认为汶川地震不是水库蓄水造成的^[3].Ge等通过二维模型计算得到在震源深度处库仑应力最大可以达到0.05MPa,结合龙门山地区荷载增长率小于0.005 MPa/yr,因此认为紫坪铺水库使得汶川地震提前10年到100 年到来^[6].Deng等同样通过计算弹性效应和扩散效应讨论了此问题，他们认为在10km 和15km 深度，库仑应力分别为10^-2和10^-3bar (1k Pa和0.1k Pa)，这个量级的库仑应力小于固体潮效应，可以忽略掉，因此紫坪铺水库对汶川地震的影响非常小^[7].Kalpna Gahalaut和V.K. Gahalaut 根据水位变化，计算水库载荷的Boussinesq解和水的扩散效应，得到在震源位置处的库仑应力为-0.1k Pa,紫坪铺水库对汶川地震没有任何触发效应^[8].周斌等通过二维孔隙弹性模型分析了紫坪铺水库诱发地震的时空演化与库水加卸载及渗透过程的关系，认为水库诱发地震与库底岩体有效应力的变化密切相关^[9].

以上学者的研究结果和结论存在比较大的差异，对这些研究成果进行分析，我们认为主要存在以下三种情况可能对计算结果有比较大的影响，而对此问题的详细探讨也许能更深入地得到紫坪铺水库对汶川地震的影响程度.(1) 二维模型和三维模型的差别.二维模型将紫坪铺水库的效应夸大了，因此势必会使计算结果夸大，如文献^[6].(2) 模型中采用的断层分布和计算库仑应力时采用的接收断层模型.水库和断层的位置与地震触发效应有着密切关系^[10-11].从波形反演得到的单倾角断层面结果看(USGS或中国地震局)，紫坪铺水库处于汶川逆断层地震的上盘.因此，单独从弹性效应来看，紫坪铺水库蓄水使得断层更加稳定，放水使得断层更加危险^[12].但是，从地表断裂分布来看，紫坪铺水库坐落在龙门山地区的三条断裂(汶川-茂县逆断裂、映秀-北川逆断裂和灌县-安县逆断裂)之一的映秀-北川逆断裂^[13]，如图 1.而且三条断裂带呈“铲" 状分布^{[1, 14-15]}.如果“铲"状断层贯通震源和紫坪铺水库库区，那么由“铲"状断层(文献^[5-6])和单倾角断层(文献^[7-8])这两种情况计算得到的库仑应力必然存在差别.因此进一步探测断层的产状对研究结果有重要意义.(3) 断层存在的效应.以上学者的研究都是建立在水库的弹性载荷效应和水的扩散效应基础上，弹性载荷效应和水的扩散效应分别通过Boussinesq解和扩散方程单独进行计算.但是实际情况中，采用联合两种效应的孔隙弹性模型进行计算无疑更适合.同时，计算中往往忽略了断层破裂带的效应，因为断层破裂带，无论是弹性常数还是扩散常数都和正常地层有着很大的差别.而以往学者单独进行计算，只能整体变化材料参数，很难考虑进去断层存在的效应.

图 1

Fig. 1

图 1 研究区域的构造背景.WM汶川一茂县逆断裂；YB映秀一北川逆断裂；GA灌县一安县逆断裂[13].小震震源机制解为2004年8月到汶川地震发震前2.2级以上地震[16-17].虚线为本次计算模型的区域.蓝色五角星为汶川地震的震中位置，蓝色实线圈定的范围为紫坪铺水库正常蓄水位的水域范围 Fig. 1 Tectonic features in the Longmenshan area. WM - Wenchuan-Maoxian fault. YB-Yingxiu-Beichuan fault. GA- Guanxian-Anxian fault[13]. The magnitude of the small earthquakes shown with focal mechanism is larger than 2.2 from the August,2004 to the occurrence of Wenchuan earthquake[16-17]. The dotted line shows the research region in the numerical model.The blue real line shows the region of Zipingpu reservoir.The blue pentagram indicates the epicenter of Wenchuan earthquake.

本文中我们根据三维孔隙弹性模型，采用有限元方法考虑断层存在的影响，详细分析随着水库水位变化，紫坪铺水库对汶川地震的影响.对于断层模型，我们假设存在一条贯通水库库区和震源位置的断层带，这样充分考虑了断层存在的效应.而计算库仑应力时，我们仍然采用波形反演得到的单倾角断层面模型作为接收断层.

孔隙弹性体的一般本构关系为^[18-19]

(1)

其中G是剪切模量(G= ，单位:Pa,E和υ 分别是弹性模量和排水泊松比)，υ_u 是非排水泊松比，B是Skempton 系数，p为孔隙压力，单位:Pa,δ_ij是Kronecker delta函数.相比弹性体，孔隙弹性体中考虑了孔隙压力和弹性的耦合关系.同样，对于水的渗流仍需要考虑弹性压力的效应，这里我们将水的渗流问题考虑为扩散方程问题，即暂不考虑弹性压力对水的渗流的影响^{[11, 20]}:

(2)

(3)

其中k是介质的渗透率，单位m²，η 是流体动力学黏滞性系数，单位Pa·s,c为水力学扩散系数，以下简称扩散系数，单位m²/s,本文中将c作为计算参数.

在小形变的几何方程和平衡方程约束下，施加合理的边界条件就可以进行求解.

几何方程:

(4)

平衡方程(忽略体力项):

(5)

我们研究的时间段为2004年到2008年汶川地震发震时刻，根据这一时间段的水库水位变化给定相应的边界条件.对于库区水位以上，给定力的边界和孔隙压力边界均为0，对于库区水位以下地区给定力的边界和孔隙压边界条件为:

(6)

ρ_w 为水的密度，g=9.8m/s² 为重力加速度，h(x，y，z，t)为t时刻(相对计算初始时刻，即2004 年1 月1 日)水库水位减去不同地点的地形高程值，p(x，y，z，t)为对应的水库蓄水载荷或孔隙压力.对于力的边界:模型侧面和底面均为滑动边界条件.对于孔隙压力:模型侧面给定为等渗流梯度，底面固定孔隙压力为0.

图 2为建立的有限元模型及网格剖分图，模型考虑的深度为25km,范围为103.3°E-103.7°E，30.8°N-31.2°N，如图 1 虚线所示.地表地形采用SRTM90m 数字高程数据(http://srtm.csi.cgiar.org/2008/version4[2008-12-10]).采用三棱柱单元进行网格剖分，并对紫坪铺库区进行局部加密，整个模型节点数为520548，总单元数为995220.模型中不同颜色代表不同的介质，虽然我们的模型可以进行分层计算，但是为了考虑断层存在的影响，本文中我们考虑整个模型有两种材料组成，即断层和非断层地区.断层的存在在孔隙弹性体中表现为弹性模量比较低，在渗流过程中表现为扩散系数比较高.模型中三条断层大致以60°倾角延伸至15km,因此对于映秀-北川逆断裂在深部大致经过汶川地震震源.我们建立若干模型分别讨论断层存在对计算结果的影响(表 1).

图 2

Fig. 2

图 2 本研究所建立的有限元模型. 灰色五角星为汶川地震震源 Fig. 2 Finite element model for this study.The hoar pentagram indicates the hypocenter of Wenchuan earthquake

表 1(Table 1)

表 1 各个模型中计算参数的选取 Table 1 The parameters used in different models

表 1 各个模型中计算参数的选取 Table 1 The parameters used in different models

库仑应力的变化量定义为^[21-22]

(7)

Δτ 为所考虑的断层面上剪应力变化量，μ 是摩擦系数，Δσ_n 为断层面上的正应力改变量(拉伸为正)，Δp为孔隙压力改变量.本文中我们重点考虑紫坪铺水库从建成蓄水到汶川地震发震时刻这一时间段内对龙门山地区的影响，特别是对震源位置处的影响. 对于水库蓄水产生的载荷效应和渗流孔隙压力，相对原来的构造背景均为应力或者孔隙压力的改变量.因此，利用式(1) 和式(2) 计算出来的为应力张量和孔隙压力的改变量，而根据式(7) 计算得到的便是由于水库存在相对于原来构造背景的库仑应力改变量.

汶川地震后，国内外学者都在第一时间进行波形反演给出了地震的震源破裂模型，各个模型在断层参数上差别不大^{[15, 23]}(http://earthquake.usgs.gov/eqcenter/eqinthenews/2008/us2008ryan/finite_fault.php[2008-12-10]).但是对于震源深度，各个研究结果差别较大.之后，四川省地震局根据已有的所有资料又重新进行了定位，给出了更可靠的震源位置.本文中我们选取USGS 的断层参数，同时采用新的地震定位结果.计算库仑应力采用的断层参数如表 2.需要说明的是对于模型1-3，摩擦系数均为0.6，而模型4中摩擦系数采用了不同的值.

表 2(Table 2)

表 2 计算库仑应力所采用的断层参数 Table 2 The fault parameters for calculating the Coulomb stress

表 2 计算库仑应力所采用的断层参数 Table 2 The fault parameters for calculating the Coulomb stress

图 3-6为采用不同模型在震源位置处得到的孔隙压力和库仑应力变化图.首先根据模型1 的计算结果来分析孔隙弹性模型计算结果.模型1 中我们通过改变断层部分的弹性模量来分析其对库仑应力计算的影响，因为扩散系数同时也与弹性模量有关，为了使问题简单化，突出弹性模量改变的影响，这里我们保持扩散系数为常数.可以看出在紫坪铺水库2005年9月蓄水之前，水库渗流作用引起的孔隙压力非常小，几乎可以忽略.而随着蓄水时间的增加，渗流引起的孔隙压力也随之升高.同时，可以看出，由于水库蓄水引起的载荷效应使得震源位置处的库仑应力降低，断层变得更安全.但是对于库仑应力的计算，孔隙压力总是使得库仑应力增加，断层变得更危险.因此这两种效应的结合就会使得结果变得较为复杂，从计算结果看，随着孔隙压力的增加，负值的库仑应力逐渐增大，在汶川地震发震前，库仑应力变为正值.

图 3

Fig. 3

图 3 模型1计算得到的震源处孔隙压力和库仑应力随水库水位变化图 Fig. 3 The pore pressure and Coulomb stress vary with the water level at the hypocenter of Wenchuan earthquake for the model 1

图 4

Fig. 4

图 4 模型2计算得到的震源处孔隙压力和库仑应力随水库水位变化图 Fig. 4 The pore pressure and Coulomb stress vary with the water level at the hypocenter of Wenchuan earthquake for the model 2

图 5

Fig. 5

图 5 模型3计算得到的震源处孔隙压力和库仑应力随水库水位变化图 Fig. 5 The pore pressure and Coulomb stress vary with the water levelat the hypocenter of Wenchuan earthquake for the model 3

图 6

Fig. 6

图 6 模型4计算得到的震源处孔隙压力和库仑应力随水库水位变化图 Fig. 6 The pore pressure and Coulomb stress vary with the water level at the hypocenter of Wenchuan earthquake for the model 4

从模型1中可以看出，从弹性角度来看，断层弹性模量的改变对震源位置的库仑应力影响很小，因此在后续模型中断层带的弹性模量均取为12.5GPa.

模型2是保持断层带的扩散系数为5.0m²·s^-1，改变非断层地区的扩散系数计算的结果.当非断层区的扩散系数增大时，所对应的孔隙压力增大速率明显增高(图 4a)，当扩散系数小于0.01 m² ·s^-1 时，震源位置处的孔隙压力几乎不变化.库仑应力变化情况与之前分析的一样，不过随着扩散系数的增大，震源位置处不同时间点的库仑应力也逐渐增大. 同时，也可以看出由于扩散作用，当水库水位在最高时，震源位置处的孔隙压力不在最高点(图 4a 红线)，而是有一个滞后时间，这与震源深度和扩散系数有关.而当蓄水量经过最高水位后又重新经历一个低水位时，载荷效应最小，而同时孔隙压力有可能达到最大，因此，此时的库仑应力也可能最大，也是最危险的时候.由此也可以看出，对于水库诱发地震，在水库蓄水位最高处不一定对应着诱发地震最活跃的时刻，而是在蓄水量经过最高水位后的一段时间里.

模型3是保持非断层区的扩散系数为0.1m²·s^-1，改变断层区的扩散系数.这与模型2 中对震源位置处产生的孔隙压力和库仑应力效果一样，随着断层扩散系数的增大，震源位置处变得更危险.

目前，除了对龙门山地区有一些跨断层的渗透率测量外^[24]，还没有公开发表的关于龙门山断裂带扩散系数的参考数值，为了一般起见，我们参考了前人关于水库地震计算过程中选取的部分扩散系数^[25].最终选取的模型4 具有一般意义，其中弹性模量和泊松比是通过地震波的结果进行换算得到^[26].在此模型中为了分析摩擦系数对库仑应力计算的影响，分别给定不同的摩擦系数.我们计算得到的震源位置处的孔隙压力和库仑应力随水库水位变化结果如图 6所示.可以看出在2008年之前，震源处的库仑应力基本上为负值，而随着孔隙压力的增加，在汶川地震发震前，库仑应力均达到了正值.随着摩擦系数的改变，库仑应力有所改变，可以看出，摩擦系数在弹性载荷效应大(对应蓄水量达到最大) 或孔隙压力大的时候影响比较明显.由公式(7) 也可以看出，对于库仑应力的计算结果，改变摩擦系数主要影响到了正应力和孔隙压力两项.若摩擦系数在0.6以上，在汶川地震发震前，库仑应力达到约+1.0k Pa.根据前人的一些研究结果，0.01 MPa的库仑应力改变就有可能触发地震^[21]，但是是否更小量级的库仑应力改变也能触发地震呢?这个问题仍需要不断深入的研究.

图 7和图 8为根据模型4计算得到的汶川地震发震时刻不同深度和不同剖面的孔隙压力和库仑应力分布图，其中摩擦系数取为0.6.结合从2004年8 月到汶川地震发震前2.2级以上小震的分布情况来看，紫坪铺库区周边的小震震源分布几乎都在库仑应力增加的区域，因此可以说明这些小震发生的原因与紫坪铺水库有直接关系.同样，汶川地震的震源位置也处于库仑应力增加的区域，但是库仑应力的增加量为1k Pa左右.孔隙压力的分布状态大致以紫坪铺库区为中心随着距离的增加而逐渐减小(图 7a,7c和图 8a,8c).库仑应力由于受到孔隙压力的影响，在地层浅部，只有在库区中心位置库仑应力为负，其它地区均为库仑应力增加的区域；而在地层深部，库区中心位置下部库仑应力为负的区域逐渐增大，相应的库区周边库仑应力为正的区域逐渐缩小(图 7b,7d和图 8b,8d).同时可以看出，在断层带地区计算出来的库仑应力也较高，这是由于断层区高扩散系数的存在使得断层区的孔隙压力相对较高导致的结果(图 8b,8d).因此，如果断层破碎带连通了震源和库区，那么在震源位置处可能会更危险.

图 7

Fig. 7

图 7 模型4计算得到的汶川地震发震时刻不同深度上的孔隙压力和库仑应力分布图，其中摩擦系数取为0.6.蓝色五角星为汶川地震震中，(b)图中地震震源机制解为研究区域内震源深度小于10 km、震级大于2.2级的所有地震，(d)图中地震震源机制解为研究区域内震源深度大于10 km、震级大于2.2级的所有地震 Fig. 7 Distribution of the pore pressure and Coulomb stress at different depth when the occurrence of the Wenchuan earthquake for the model 4 (the friction coefficient is 0.6).The blue pentagram indicates the epicenter of Wenchuan earthquake.The depth of the earthquakes shown in Fig. 7b and Fig.7d is less than 10 km and more than 10 km,respectively.The magnitude of the small earthquakes shown with focal mechanism is larger than 2.2 from the August,2004 to the occurrence of Wenchuan earthquake

图 8

Fig. 8

图 8 模型4计算得到的在汶川地震发震时刻不同剖面上的孔隙压力和库仑应力分布图，其中摩擦系数取为0.6.白色五角星为汶川地震震源位置，白色圆圈为距离剖面0.06°范围内的大于2.2级地震的震源位置投影在剖面上的位置 Fig. 8 Distribution of the pore pressure and Coulomb stress along different profiles when the occurrence of the Wenchuan earthquake for the model 4 (the friction coefficient is 0.6).The white pentagram indicates the hypocenter of Wenchuan earthquake. The white circles shown in the profiles indicate the projection position of earthquake (≥2.2) in the range of 0.06 degree

本文中采用三维孔隙弹性模型，根据紫坪铺水库开始蓄水到汶川地震发震前这一段时间内的水位变化，来研究其对库区周边的影响.如前所述，断层的产状对最后的计算结果会产生比较大的影响，但是目前为止，断层深部具体情况还需要不断的探索，没有肯定的结果.因此在本文中我们没有讨论断层的形态对库仑应力计算结果的影响，而在模型中加入了断层的大致分布，通过探讨不同的参数变化来分析断层对震源处的影响.初步得到如下结论: 从弹性角度来看，断层地区的弹性模量对震源位置的库仑应力计算影响不大.断层和非断层地区的扩散系数对震源位置的库仑应力计算影响较大，即使非断层区的扩散系数较低而断层区的扩散系数很大也有可能造成震源位置处的库仑应力增大.分别给定断层和非断层区一般性的参数计算得到，在汶川地震发震时刻，震源位置处的库仑应力变化量为+1k Pa左右，这表明紫坪铺水库使得汶川地震的发震断层更加危险.但是否这个量级的库仑应力能够触发汶川地震还需要进一步探讨.通过分析库区周边的小震分布，发现小震几乎都分布在库仑应力增加的区域.因此，库区周边小震的发生应该与紫坪铺水库有直接关系.

感谢两位匿名审稿人提出的建设性意见.