位于川西地区的紫坪铺水库于2005年9月首 次蓄水，在2006年12月蓄水达到最高值海拔875 m. 此后的两年半，水位一直保持在817 m以上(Lei，2011；雷兴林等，2008).蓄水2.7年后，2008年5月12日M_w7.9级汶川地震在龙门山断裂带上发生，震中距离库区最近距离仅约6 km(见图 1).对其震源位置和深度，不同机构或研究者给出的结果略有不同：美国地质调查局(USGS)为19 km(http://earthquake.usgs.gov/earthquakes/eqinthenews/2008/us2008ryan/)；陈九辉等(2009)为18.8 km，中国地震局(CEA)公布数据为 14 km(http://www.csndmc.ac.cn/newweb/data/).地震发生后引发了紫坪铺水库蓄水是否触发了汶川地震的激烈争论(Klose，2008；Ge et al., 2009；Kerr and Stone， 2009； Deng et al., 2010；Gahalaut and Gahalaut， 2010；Lei，2011；雷兴林等，2008；周斌等，2010；孙玉军等，2012).雷兴林等(2008)，Lei(2011)率先计算了由于紫坪铺水库蓄水的重力加载和孔隙压力扩散引起的库仑应力变化，结果显示，在水库下方地壳10 km以上库仑应力增长达100 kPa以上，在震源区深度(按CEA结果：14 km)库仑应力增长达几十千帕量级.Ge等(2009)的研究认为水库蓄水导致震源处高达50 kPa 的库仑应力变化，如此量级的库仑应力增长相对于龙门山断裂带1~2 kPa/a的长期应力累积水平意味着将潜在的汶川地震提前了数十 年到数百年.孙玉军等(2012)的计算结果表明在汶川地震发震时刻，震源处的库仑应力变化量达+1 kPa，增加了汶川地震的发震危险性.周斌等(2010)通过二维孔隙弹性模型分析了水库诱发地震与库底岩体有效应力的关系，计算了与紫坪铺水库有关的库水加卸载及渗透过程，并分析了其与库区小震(深度在4~10 km范围内)的时空演化关系，但并未涉及汶川地震震源深度区域的讨论.然而，Gahalaut K和Gahalaut V K(2010)和Deng 等(2010)的研究却认为紫坪铺水库蓄水导致的汶川地震震源处的库仑应力加载为负值或仅10 Pa左右，对促进汶川地震发生的影响微乎其微甚至相反.

图 1

Fig. 1

图 1 紫坪铺水库和汶川地震区构造背景 红点表示从2004年9月到2008年4月紫坪铺水库周围M>1.0级地震.蓝色、亮蓝和紫色三角分别表示用于重新定位紫坪铺水库周围小震的紫坪铺、成都和四川地震台网台站位置.深红色线表示汶川同震破裂地表行迹.蓝色、绿色和黄色五角星分别表示USGS，CEA和陈九辉等(2009)发布的汶川地震主震位置.亮蓝色不规则区域表示出紫坪铺水库区域.Fig. 1 Tectonic setting of the Zipingpu Reservoir(ZR) and Wenchuan Earthquake(WE) The red dots show M>1.0 earthquakes near the ZR between September 2004 and April 2008. The Zipingpu，Chengdu， and Sichuan seismic networks used to relocate the earthquakes are shown as blue，light blue， and purple triangles respectively. Dark red lines indicate surface traces of WE coseismic rupture. The mainshock location from United States Geological Survey(USGS)，China Earthquake Administration(CEA)， and Chen et al.(2009)are indicated by blue，green， and yellow stars，respectively. The light blue irregular shaped area indicates the ZR.

有关水库诱发地震的研究早已有之(Hagiwara and Ohtake， 1972；Gupta and Chadha， 1995；Talwani，1997；Guha，2000)，最早认知的诱发机制为：水 库的水体重力加载触发了地震(Bell and Nur， 1978；Simpson，1986；Simpson et al., 1988).然而，观测结果表明库区小震密集发生时间总是在水库蓄水后的几年，与水库重力加载的弹性响应在时间上不相符(Talwani and Acree， 1984；Talwani，1997).Roeloff(1988)评估了在均匀半空间孔隙弹性介质中水体重力加载效应和孔隙压力扩散效应，Kusalara和Talwani(1992)计算了水库蓄水的弹性响应和孔隙压力扩散引起的应力场变化，结果均显示，弹性响应对断层的影响依赖于先存断裂展布，既可以使得库区附近断裂趋于稳定，也可以使得断裂趋于不稳定，但是孔隙压力的扩散肯定会导致断裂趋于不稳定.因此，在多数情况下，水库蓄水导致的区域孔隙压力扩散才是水库诱发地震的主要原因(Talwani，1997).

水库诱发地震是孔隙压力场与弹性骨架应力场相互耦合作用的物理过程.孔隙压力的改变会作用于介质弹性骨架造成应力场改变，同时应力场的改变也会反作用于孔隙压力场.本文基于完全耦合孔隙弹性理论，考虑孔隙压力场和骨架应力场的相互影响，利用二维有限元模型模拟紫坪铺水库蓄水导 致的龙门山断裂带上库仑应力随时间的演化，探讨紫坪铺水库蓄水与汶川地震发生的关系.

库仑破裂准则被广泛应用于判断先存断裂面上应力变化是否有利于发生摩擦滑动(King et al., 1994；Cocco and Rice， 2002).库仑应力变化(ΔCFS)定义为:

许多研究在计算库仑应力时，由于无法知道孔隙压力增量，不得不将(1)式简化为：ΔCFS=σ_s+μ′σ′ _n，其中σ′ _n为有效正应力，而μ′为视摩擦系数.在此情况下，视摩擦系数μ′已经不是岩石的固有介质常数，它融合了岩石摩擦性质和孔隙压力改变两部分的影响，虽然许多研究中经常令μ′=0.4或0.6；但理论上来说视摩擦系数的范围为：-∞≤μ′≤+∞(Beeler et al., 2000).而本文利用公式(1)，独立考虑孔隙压力的影响，需要利用岩石固有摩擦性质，因此选择岩石固有摩擦系数(0.6~0.85)的中间值μ=0.7.

本文研究紫坪铺水库蓄水对龙门山地区的影响，即研究由水库蓄水产生的载荷效应和渗流孔隙压力增量，均为相对于视为平衡态的背景构造应力和孔隙压力场的改变量.公式(1)中的σ_n+P为有效正应力(σ_effn).与断层错动方向一致的剪切应力和有效正应力的增加使得ΔCFS增加，促进断层破裂.本文通过计算得到剪切应力、正应力和孔隙压力；并通过讨论剪切应力和有效正应力的演化，来分析断裂带上ΔCFS的演化过程.

King等(1994)在对L and ers地震及其余震以及附近地震的触发研究中发现100~600 kPa的ΔCFS升高可触发地震；Stein(1999)研究了圣安德烈斯断层上大震引起的ΔCFS对区域地震危险性的影响，进一步得出10 kPa的应力升高，能使得地震活动性提高20%~25%，统计上十分显著.Cochran 等(2004)根据对全球较大地震事件(M>5.5)与强固体潮力之间相关性的统计研究，发现浅层逆冲性质地震与固体潮幅值有强相关性，而此幅值仅为5 kPa.这些研究成果为我们研究紫坪铺水库蓄水可能造成的地震触发提供了重要参考.

在饱和孔隙弹性介质中的弹性骨架应力场和孔隙流体压力场是相互耦合作用的，完全耦合孔隙弹性控制方程为(二维平面应变假设)：

通过计算水库蓄水导致的应力场和孔隙压力场的改变进行区域断层稳定性和地震活动性评估是经典孔隙弹性问题.在水库蓄水初期，水体重力加载的弹性响应影响区域岩石骨架变形和孔隙流体运移；同时水库作为孔隙压力源，形成区域孔隙压力梯度.随着时间进程，孔隙压力梯度的存在导致孔隙压力由水库向周围地区逐渐扩散，使得区域弹性应力场和孔隙压力场随时间演化，则区域内断层稳定性也随时间演化.

紫坪铺水库紧邻汶川地震主破裂线，汶川地震地表破裂走向约为N(42° ± 5°)E(Xu et al., 2009)(见图 1).从剖面上看，龙门山断裂带倾向西北，断层上部倾角相对陡，在底部以近似水平的角度并入低速层(Burchfiel et al., 2008；Jia et al., 2010).本文设计的二维平面应变有限元几何模型为横穿紫坪铺水库，与汶川地震地表破裂近垂直、横跨青藏高原和四川盆地的一个剖面.为消除边界效应影响，设计模型宽280 km、厚 100 km.因专注于水库蓄水到汶川地震前的时间范围，可以不考虑岩石圈内介质黏弹性的影响.模型区域分成地壳、地幔两部分，地幔部分为纯弹性区域，地壳部分为孔隙弹性区域.根据库区资料，紫坪铺水库形态约为不规则椭圆型，长轴大致沿主震地表破裂走向约5 km长，短轴大致与地表破裂垂直约2 km宽.根据紫坪铺水库的蓄水历史(雷兴林等，2008)，最高水位为877 m，最低水位约为830 m，库底高度为728 m，蓄水以来保守估计的最低水深为100 m.因此，模型中将水库形态简化为横跨2 km，深100 m的半椭圆型，嵌入地表相应位置，模型共有24964个节点，24570个平面应变单元(图 2).孔隙弹性介质的渗透率随深度的增加而减小.水库附近单元尺度约为10 m，模型其他区域单元尺度约为800 m.

图 2

Fig. 2

图 2 有限元模型 (a)二维模型分成地壳和地幔两部分；地壳部分进一步分成5层，分别对应不同孔隙弹性介质，参数如表 1所示.(b)模型网格划分成24964个节点，24570个单元.(c)对应断层位置和水库附近的加密网格.红色箭头表示紫坪铺水库加载的位置.Fig. 2 The FEM configuration (a)The 2D domain is partitioned into mantle and crust regions. The crust is further partitioned into 5 layers. Corresponding material properties are summarized in Table 1.(b)FEM mesh contains 24964 nodes and 24570 elements.(c)Blow-up view of the FEM mesh near the fault and reservoir. Red arrows indicate the location of the ZR load.

渗透率为岩石的固有介质属性.Ingebritsen和Mannning(1999，2010)根据热液模型和由流体驱动的变质反应过程估算了大陆地壳随深度变化的渗透 率.研究显示，在构造活跃的大陆地壳中，平均渗透 率随深度递减，遵从公式logk≈-14-3.2 log(-z)，其中k为渗透率(m²)，z为深度(km)，向上为正.据 此，在近地表处渗透率范围为10^-11~10^-13 m²，在地壳约40 km深处渗透率达到10^-19 m².Shmonov 等(2002)通过野外和实验室测量也得到类似的关系式：logk≈-12.56-3.225(-z)^0.223，在近地表到 约40 km地壳深处渗透率范围约为10^-14~10^-20 m².

本研究利用有限元软件Abaqus software(2009)进行模拟计算，将渗透率转换为液导率作为模型的输入参数.液导率与渗透率的关系为：

表 1(Table 1)

表 1 模型各部分的介质参数Table 1 EParameter of mechanic character model 泊松比 杨氏模量(Pa) 液导率(m·s-1) 花岗岩不排水体积模量(Pa) 参考文献 Crust_1 0.25 3.75×1010 1.0×10-8 4.2×1010 Wang(2000)， Ingebritsen et al.(2010) Crust_2 0.25 3.75×1010 2.0×10-10 4.2×1010 Crust_3 0.25 3.75×1010 3.6×10-11 4.2×1010 Crust_4 0.25 3.75×1010 1.0×10-11 4.2×1010 Crust_5 0.25 3.75×1010 2.5×10-12 4.2×1010 Mantle 0.36 1.4×1011 Turcotte and Schubert(2002)

表 1 模型各部分的介质参数Table 1 Parameter of mechanic character model

地壳岩石的孔隙弹性参数来自Westerly 花岗岩的孔隙弹性参数.因与地幔特征松弛时间相比，紫坪铺水库注水到汶川地震发生的时间可近似为瞬态，因此可忽略地幔的黏弹性反应.地幔近似为弹性介质，利用了橄榄岩的物性参数(Turcotte and Schubert， 2002).模型中地幔和地壳的介质参数总结在表 1中.

根据紫坪铺水库的蓄水历史(雷兴林等，2008)，蓄水以来保守估计最低水深为100 m，因此库区中心位置施加h=100 m深的重力加载和p=ρgh(p=1 MPa)的水头压力.非库区处地表自由、孔隙压力p=0.由于库区为半椭圆型嵌入边界，则沿库区边界按深度施加水压边界条件.由于模型横向尺度足够大，可近似设定模型两侧边界横向位移为0、孔隙压力p=0.地壳底部流量为0，即假设地幔与地壳之间没有流体交换；地幔底部，垂向位移为0.

图 3给出了孔隙压力随时间演化过程.随着紫坪铺水库开始蓄水，初始弹性响应之后，孔隙压力向周围扩散，周边地区、特别是库底的孔隙压力逐渐升高.附近的龙门山断裂带由浅至深逐渐受到以水库为中心的孔隙压力扩散影响，两年多后，断裂带上CEA发布震源位置处孔隙压力增大在40 kPa以上，USGS发布震源位置处孔隙压力增大在30 kPa以上(图 4).

图 3

Fig. 3

图 3 紫坪铺水库蓄水后孔隙压力的时空演化 黄色曲线表示龙门山断裂带的位置，箭头表示水库加载位置.Fig. 3 Spatial evolution of pore pressure after impoundment of ZR The yellow curve denotes the location of the LSF.

图 4

Fig. 4

图 4 模拟得到的孔隙压力随时间变化曲线，右方小图显示点在模型中位置Fig. 4 Evolution of pore pressure with time predicted by FEM model.(a)Observation locations and (b)predicted pore pressure changes at the corresponding positions

图 4显示了模型中6个位置的孔隙压力随时间的变化曲线，6个位置分别为水库正下方5 km、10 km和14 km深处，以及CEA、USGS和陈九辉等(2009)发布的汶川地震震源位置.总的来说，各处孔隙压力都随时间增长.由于地壳中渗透率随深度减小，则距离库区越远，位于地壳越深处，孔隙压力的增长越缓慢.

在震源位置，孔隙压力随时间呈现瞬时增大，短时微弱减小，然后持续增长的微起伏过程，这是不排水孔隙压力和扩散孔隙压力共同作用的结果.因为水库重力加载为弹性瞬时加载，引起孔隙体积减小，此时孔隙中流体来不及流出(不排水反应)，导致瞬时孔隙压力急剧升高.紧接着进入排水阶段，流体从弹性孔隙介质中流出，孔隙压力些许下降，下降时间取决于介质中孔隙压力梯度和渗透率.之后，水库产生的孔隙压水头将会驱动孔隙压力持续向周围扩散，导致周围地区孔隙压力又逐渐升高.其影响与孔隙压力源、渗透率，以及观测点与水库距离相关.

从定量分析结果看，紫坪铺水库正下方的孔隙 压力增加随深度增加而减小，如：蓄水10个月后库 底5 km深处孔隙压力达到164 kPa(图 4中浅蓝色线)，10 km深处达到100 kPa(图 4中深蓝色线)，14 km深处达到75 kPa(图 4中紫色线).同深度上，水库正下方孔隙压力总是比周围地壳内孔隙压力高(图 4中紫色线和棕色线)，如：震前时间库底14 km深处孔隙压力约为101 kPa(紫色线)，而CEA发布的震源处为41 kPa(棕色线).震前USGS发布的震源处孔隙压力增加31 kPa，陈九辉等(2009)发布的震源处增加35 kPa.在地壳越深处、距离水库越远，孔隙压力的增长越缓慢，但会一直保持随时间缓慢增长的主控趋势.

库仑应力的计算会受到断层倾角的直接影响.龙门山断裂带的几何模型由龙门山地区的地震反射剖面简化得到，上部断层倾角相对陡，在底部断层近似水平并入低速层(Burchfiel et al., 2008；Jia et al., 2010).断裂带几何模型的其他约束还包括，根据Li等(2013)的断层钻探数据，由钻探与断层破裂面的相遇位置和地表破裂出露位置控制，得到近地表断层角度约为61°；USGS发布的震源破裂面角度35°；以及CEA、USGS和陈九辉等(2009)发布的震源深度.由此，龙门山断层形态可用一个连续函数

y= 2a π arctan x b 近似描述.其中参数为：a=31.6，b=9.7.断层形迹如图 5中的黄色曲线，

在近地表角度为61°，14 km深处角度为 35°，19 km深处角度约为 22°，逐渐于27 km深处趋于水平.

在此基础上，根据相应深度上的断层倾角，利用有限元模型计算了龙门山断裂带上剪切应力和正应力分布.正应力和孔隙压力的叠加得到有效正应力.水库蓄水后，库底以及水库周围地壳内的孔隙压力随时间增长，使得龙门山断裂带上有效正应力随时间增长(图 5)，影响范围逐渐扩大，有利于断层滑动.而剪应力变化主要源自水库自重加载，随时间变化不大，只水库西侧小范围剪切应力为正.

结合剪切应力和有效正应力，可计算断层周围的库仑应力(ΔCFS)演化(见图 5).ΔCFS分布是一个复杂但有一定规律的图像.水库西侧ΔCFS基本为正，对断层活动起促进作用；而水库东侧为负，对断层活动起抑制作用.水库蓄水初始，ΔCFS增加局限在水库下方及西侧很小区域.龙门山断裂带上地壳浅部ΔCFS增量为正，但深部在汶川地震震源深度的ΔCFS增量为负，说明水库的重力加载弹性响应对龙门山断裂带深处的影响是卸载作用.但2个月之后，随着孔隙压力水头的扩散，龙门山断裂带上的ΔCFS正值区域逐渐扩大，逐渐向断层深处影响，而负值范围不断缩小.直到震前，整条断层都受到正值ΔCFS的影响.结合孔隙压力扩散的特性，ΔCFS随时间增长的整体趋势不变，某一位置的ΔCFS值即使初始时刻为负，在未来某时刻仍可能由负转正，并不断增大.因此，水库蓄水作用的总趋势是向更大范围促进断层活动的方向发展.

图 5

Fig. 5

图 5 龙门山断裂带附近剪切应力、有效正应力和库仑应力分布及随时间演化 计算所用断层倾角采用相同深度上拟合的龙门山断裂倾角. 图中自上到下小图，分别表示水库蓄水、水后2个月、1年及2.7年(汶川震前)的应力表示.Fig. 5 Shear，normal， and Coulomb stress distributions on fault plane of LSF The stresses are plotted using the rainbow legend， and magenta color indicates the stresses below the negative limit. The fault dip angle of the LSF is adopted for the stress calculation on a hypothetical fault plane at the same depth. From top to bottom the figs display the stresses distributions at the ZR impoundment，2 month later，1year later and 2.7 year later(before WE)respectively.

在Lei等(2010)和Ge等(2009)的研究中，以震 源处10~50 kPa的ΔCFS增长作为判断紫坪铺水库对汶川地震触发作用的门槛.CEA发布的震源位置位于本研究拟合连续断层面上，ΔCFS增长超过25 kPa.而USGS和陈九辉等(2009)给出的震源位置位于断层面之下，ΔCFS增长仅为2 kPa左右.但若考虑到震源位置的误差，在深度或者横向上移动1~2 km，则ΔCFS结果可能稍有差异.若以CEA发布震源位置为准，则初始破裂直接触发了断层上更大范围的活动，引发汶川地震；若以USGS和陈九辉等(2009)给出的震源位置为准，则初始破裂很可能是个区域构造小震，而断层面在长期构造应力加载累积下已经积累了一定应变能，再加上紫坪铺 水库蓄水对断层整体的显著加载作用(在深度13 km 以上ΔCFS增加超过30 kPa，深度5 km以上超过 100 kPa)，则此构造小震通过改变局部应力场，也可能触发附近断层面上的失稳，导致大震发生.则无论哪种情况，从整体来看，紫坪铺水库蓄水对龙门山断裂带的发震危险性起促进作用.

断层的几何形态是计算断层面上的ΔCFS分布的关键参数，然而目前对于龙门山断层在地下深处的展布仍然没有定论.虽然本文以钻探、地质调查、震源机制等资料为约束，给出了较合理的连续断层形态，但与实际情况仍可能有一定偏差.前人在计算ΔCFS时，有的以USGS给出的主震破裂面角度35°和震源深度为根据，如：Deng 等(2010)，Gahalaut K和Gahalaut V K(2010)，孙玉军等(2012)；有的基于地震探测和地质调查结果，假设一个断层分布进行计算，如Ge等(2009)，雷兴林等(2008)，Lei(2011).但较小的断层角度调整都可能影响库仑应力的大小甚至性质(Zhou and Deng， 2011)，在断层展布形态有疑问的情况下，在确定深度仅对确定断层倾角计算的ΔCFS的讨论可能失于片面.此外，震源深度是另一个影响ΔCFS计算的关键参数.对于汶川地震的主震位置，不同的研究者给出了不同结果，并没有统一结论，并且地震学给出的主震深度本身也是具有一定误差范围内的结果.前人的研究中，对震源位置误差可能带来的影响也未进行深入讨论.因此，要对紫坪铺水库问题有更全面的认识，必须对断层形态不确定性和震源位置误差的影响进行更深入讨论.

本文计算了震前震源处及其附近相对于各种可能断层倾角(15°~75°)的ΔCFS曲线图和等值线分布图(如图 6a—6d).图 6a，6b给出在水库附近地壳内14 km和19 km深度剖面上相对于15°~75°断层倾角的ΔCFS分布.两条剖面图均显示，在剖面所在深度上的某点，若其与水库之间的水平距离在20 km范围内，则断层倾角越大，ΔCFS越大.在14 km深度剖面上，ΔCFS值最大处落在距离水库小于25 km，断层倾角大于30°的范围内；在19 km深度剖面上，ΔCFS值最大处落在距离水库小于30 km，断层倾角大于35°范围内.图 6c给出在横向距离水库7 km的垂直剖面上，相对于15°~75°断层倾角的ΔCFS分布.结果显示，随着深度的增加，需要相对于更高的断层倾角才能得到ΔCFS正值.观察ΔCFS正值范围，在地下6 km处，只要断层倾角大于20°，则水库蓄水影响为促进断层活动；而在地下21 km深度，断层倾角需大于40°，水库蓄水促进断层活动.整体来看，随着孔隙压力的扩散，ΔCFS正值范围将会不断扩大.因此，紫坪铺水库对于龙门山断裂带的影响以促进断层活动为主要趋势.

图 6

Fig. 6

图 6 相对于所有可能断层倾角的ΔCFS等值线分布及曲线图 (a)地壳内14 km 深处水平剖面；(b)地壳内19 km 深处水平剖面；(c)距离水库7 km处垂向剖面；(d)CEA、USGS和陈九辉等(2009)给出的震源处的库仑应力随断层倾角变化曲线；(d)中插图给出(a—c)中剖面和点在模型中的位置.Fig. 6 Coulomb stress variation vs. fault dip angle in hypocenter region (a)A horizontal profile on the depth of 14 km;(b)A horizontal profile on the depth of 19 km;(c)A vertical profile located 7 km from the reservoir;(d)Coulomb stress variations on the hypocenters given by CEA，USGS， and Chen et al.(2009); Inset diagram in(d)provides the locations of the profiles(a—c) and the hypocenters discussed in the model.

图 6d给出在CEA、USGS和陈九辉等(2009)给出的汶川地震震源处，相对于各种可能断层倾角(15°~75°)的ΔCFS变化曲线.均为单调递增曲线，显示断层倾角越大，ΔCFS越大，对断层发震的促进作用越强.按照CEA发布的震源位置，若断层倾角在30°~50°范围内，在震前ΔCFS加载会达到2~ 25 kPa.若按USGS或陈九辉等(2009)给出的震源 位置，则相对于30°~37°的断层倾角，ΔCFS为负；相对于37°~50°的断层倾角，ΔCFS为0~16 kPa.根据前人研究(King et al., 1994；Stein，1999；Cochran et al., 2004)，5 kPa量值的ΔCFS加载可以显著提高区域地震活动性概率，甚至触发地震；而当ΔCFS增加超过10 kPa时这一效应更为明显.则在三个发布的震源位置，都可能具备触发大震的应力条件.

本研究基于完全耦合孔隙弹性理论，利用二维有限元模型模拟了水库蓄水造成的区域孔隙压力场和弹性应力场的演化过程.由于完全耦合孔隙弹性理论计算的复杂性，前人的相关研究多忽略或简化了孔隙压力场和弹性骨架应力场之间的相互作用，或分别计算两场，然后相加，简化为非耦合问题，如：雷兴林等(2008)，Lei(2011)、Ge等(2009)和Deng等(2010)；或只考虑孔隙压力对弹性骨架的影响，不考虑骨架应力场改变对孔隙压力场的作用，简化为单向耦合问题，如：Gahalaut K和Gahalaut V K(2010)、孙玉军等(2012).而本研究完全耦合两场进行计算分析，更客观地反映了实际物理作用.周斌等(2010)虽然基于完全耦合孔隙弹性理论进行计算，但没有计算研究汶川震源区的应力场演化情况.

在前人研究中，Gahalaut K和Gahalaut V K(2010)和Deng等(2010)的研究认为紫坪铺水库蓄水对促进汶川地震的影响微乎其微，甚至起相反作用.对Deng等(2010)的工作，Ge 等(2011)进行了检验模拟计算，认为其在计算时仅考虑了扩散孔隙压力的贡献，并未包括不排水孔隙压力部分的贡献，导致所得孔隙压力偏低，进而得到的ΔCFS值偏低.不排水孔隙压力作用的特征是在水库加载瞬间导致区域孔隙压力瞬间增高，由此在库区附近的ΔCFS值也瞬间升高到一定值，而不会从零开始演化(见图 4).因此，Ge等(2011)认为Deng等(2010)的工作中的此处瑕疵，影响了他们对紫坪铺水库对汶川地震影响的结论.而Zhou和Deng(2011)认为Ge等(2009)计算所用的震源深度断层倾角45°偏高，因而高估了ΔCFS.然而Zhou和Deng(2011)没有提及，实际上紫坪铺水库垂直于汶川地震地表破裂形迹方向的剖面宽度只有2 km左右，而在Ge等(2009)的计算中水库加载宽度在5 km左右，这样 的设置使水库加载增大了1.5倍，进一步使得ΔCFS 计算结果偏高.这可能是影响更大的因素.此外，在地层的不同深度使用同样的液体扩散率，也会导致高估ΔCFS增长.Gahalaut K和Gahalaut V K(2010)文中的孔隙压力分布与其他研究结果相比明显偏低，而文章中并没有提及所用水库加载区域，由此可推知水库加载设定可能是导致其结论与其他相关研究不同的主要原因.即便如此，其结果显示在4 km深处，龙门山断裂带上CFS增量为正值，有利于断层失稳.

本研究在客观选择模型参数的基础上，根据大陆内活动构造地壳中液导率随深度的变化，在合理范围内进行了大量检验计算.结果显示，即使对应更低的渗透率，震源区ΔCFS随时间增加的趋势不变.只要水库蓄水形成的孔隙压力水头存在，孔隙压力扩散过程就不会停止，从而导致附近龙门山断裂带上ΔCFS不断增长，有利于断层失稳.如果水库蓄水导致的孔隙压力增加促进了龙门山断裂带上汶川地震的发生，那么同样的机制也会影响到水库周围的地震活动性.为此，我们分析了水库蓄水至汶川地震发生前水库周围的小震活动性.

2004年7月，在水库建设之前1年，为了检测紫坪铺水库周围的小震活动，相关部门建设了一个包括7个短周期地震仪的水库专用台网.地震仪之间的平均空间距离为10 km左右(台站位置见图 1).利用此台网结合四川地震局的区域地震台网数据，马文涛等(2011)利用双差定位法重新定位了紫坪铺水库周围的小震，小震的平面和剖面位置见图 1和图 7.2007年9月和2008年2月记录到两个地 震群，一个在汶川地震地表破裂的西南端下方，另一个位于紫坪铺水库东北约10 km左右.从垂直剖面看，紫坪铺西南的小震分布基本勾勒出了龙门山断裂带的浅部形态.在紫坪铺水库东南3 km、深度13~15 km处的小震群可能显示了断层的薄弱带所在.大部分的小震发生于水库蓄水后2～2.5年间，说明在此期间孔隙压力改变已经影响到龙门山断裂带地区，并引发了先存薄弱面上的地震.

图 7

Fig. 7

图 7 2004年9月—2008年4月间紫坪铺水库地区小震分布剖面图Fig. 7 Spatial distribution of the micro earthquakes observed in the ZR region during September 2004—April 2008 Events are colored based on their occurrence times.(a)，(b)， and (c)show the events in map view，in a depth profile along fault projecting N49°E， and in a depth profile across fault projecting N229°E，respectively. The red star shows the hypocenter location given by USGS.

为分析随时间变化的孔隙压力扩散效应，以紫坪铺水库蓄水为时间原点，将距离紫坪铺水库20 km内的小震以每6个月为时间段分组，将每组地震按震中与紫坪铺水库的距离排序，统计每组累积的小震数，见图 8.图 8中，每个时间段的地震曲线上，斜率的变化揭示了不同距离上小震发生密度的演化过程.例如，表示紫坪铺水库蓄水之前0～0.5年和0.5～1年的两条曲线B1和B2，最大斜率变化位于距离紫坪铺水库5~10 km范围内，说明小震发生最密集的地方在距离水库中心5~10 km范围内.同时，在水库蓄水后，每个时间段内的小震发生最密集处都在距离水库中心约10 km 以外，以时间为序，分别为：0～0.5年：9.5~12 km处(见曲线C1)；0.5～1年：13~15 km处(见C2)；1～1.5年：13~16 km处(见C3)；1.5～2年：14~17 km处(见C4)；以及2.5～2.7年(汶川地震前)：15~20 km 处(见C5).这些曲线揭示出，在紫坪铺水库蓄水后、汶川地震发生前，小震具有明显的以紫坪铺水库为中心随时间向周围扩散的趋势，而这一扩散趋势与紫坪铺水库蓄水引起的孔隙压力向周围扩散的趋势一致.

图 8

Fig. 8

图 8 紫坪铺水库地区震中与水库距离-小震累计发生个数曲线图 每种颜色的曲线代表 6个月的一个时间段，其中黑色线段代表汶川地震前的时间段仅含2个月.Fig. 8 Cumulative histogram of earthquakes occurred along distance from the ZR Colors of the curves denote the time ranges at 6 month intervals，except the black curve which represents data of the last two months prior to the WE.

此外，卢显等(2010)对紫坪铺水库台网记录到的汶川震前的小震进行重新定位，并比较紫坪铺水库蓄水前、后的小震分布，结果显示：(1)蓄水后比蓄水前小震明显增多，蓄水前小震震源深度主要分布在5～12 km 范围内，而蓄水后大部分小震震源深度都明显集中在5 km左右，小部分地震震源深度分布在10 km以下.(2)小震主要集中分布在龙门山地区的先存断裂带上，其中大部分地震位于北川—映秀断裂并收敛于紫坪铺水库的两端.此工作用另一种分析方法得到与本文中小震分析基本一致的结论，即：受紫坪铺水库蓄水影响，孔隙压力扩散，由浅至深改变了区域地壳内的应力状态，在先存薄弱面——龙门山断裂带上引发大量小地震，并改变了断层整体的应力状态.

综上所述，基于完全耦合孔隙弹性数值模拟研究和小震活动性分析，可得以下结论：

(1)紫坪铺水库的蓄水破坏了原有的区域孔隙压力平衡，形成恒定存在的孔隙压力源，孔隙压力水头逐渐向周围地区扩散，并影响到相邻的龙门山断裂带上，从浅到深使得ΔCFS不断增长.即使在初期ΔCFS变化为负的地区，在未来某时刻仍可能由负转正，持续增长.

(2)水库蓄水最初，龙门山断裂带浅部ΔCFS为正，深部地壳内的ΔCFS较小.之后孔隙压力扩散，导致断裂带上有效正应力不断增大，进而断裂带上ΔCFS持续增长并向深部发展.以至于在2年多后汶川地震前，整条断层都受到正值ΔCFS加载，深度 5 km以上加载超过 100 kPa，13 km以上超过30 kPa，而震源区ΔCFS增长值超过2~25 kPa，对汶川地震的发生有较强的促进乃至触发作用.

(3)紫坪铺水库蓄水使得龙门山断裂带整体都处于显著ΔCFS加载作用下.若初始破裂位于主震断层面上，则可直接触发断层上更大范围的活动，引发汶川地震；若初始破裂不发生在主震断层面上，仅为常规区域构造小震，则此小震的发生通过改变局部应力场，也可能触发附近断层面上失稳，导致大震发生.无论怎样，从整体来看，紫坪铺水库蓄水对龙门山断裂带的发震危险性起促进作用.

(4)地壳中ΔCFS计算受到断层几何形态的影响，计算结果显示，三个发布的汶川地震震源处及其附近地区，都是断层倾角越大，断层面上的ΔCFS值越大，对断层发震的促进作用越强.随着计算点位深度的增加，需要相对于更高的断层倾角才能得到ΔCFS正值.

(5)紫坪铺水库蓄水后，周围小震明显增多，并且小震随时间演化的模式与区域孔隙压力的演化模式一致，显示紫坪铺水库蓄水导致的区域孔隙压力改变，促进了汶川地震的发生.