水库诱发地震（RIS）是一类特殊的地震活动，主要是由于水库蓄水或排水过程引发的在一定时间内库区及其周边不大区域范围内所发生的地震活动^[1].迄今为止，全球已报道的RIS确切震例约有130例，其中，大于6.0级的有4例，5.0~5.9级的有15例，其余占85%以上的均小于5.0级.虽然大部分RIS的震级不高，但由于其震中位置一般邻近重要的水利工程设施，且震源浅、震中烈度高，往往具有很大的破坏性，可造成大坝及附近建筑物的破坏和人员伤亡，甚至会引起滑坡、坍塌等严重的次生灾害危及下游安全.因此，RIS不仅是水利水电工程研究的重要内容，也是地震学、区域构造稳定性和环境工程地质研究的重要内容之一.

RIS的发生与库水加卸载及渗透过程有着密切的联系，这一点可从印度地理研究所与德国波恩大学在印度Koyna-Warna库区的深井水位观测成果中得到启示^{[2, 3]}.近50年来，一些学者分别从库水荷载作用、孔隙压力扩散作用和润滑作用等方面探讨了RIS发生的力学机制，得到了一些有益的认识^[4~12]，但对于水库水位变化过程中，诱发地震活动在不同时、空尺度上对库体荷载及水库附加水头压力扩散的动态响应机制问题，目前尚少见报道.RIS实质上是在库水加卸载及渗透过程中，地下岩石介质变形与流体渗流耦合作用造成的岩体失稳问题，需要在一个综合考虑构造应力场、渗流场与岩体变形、失稳相互关系的地质力学框架体系内进行讨论.在此方面，雷兴林等^[12]在半无限弹性介质的假设下，利用格林函数和孔隙压力扩散方程分别求解应力场和流体压力场，定量讨论了断层库仑应力变化对库体荷载及附加水头压力扩散的静态响应问题.这一方法虽然便捷，但在处理复杂介质域时却有很大的局限性.RIS受库区构造应力环境、地质构造和水文地质结构条件的控制^{[13, 14]}，前人大多采用简化的地质模型，忽略了发震库区地质构造及水文地质结构复杂性的影响.因此，在库区地质调查的基础上，建立切合实际的地质构造及水文地质结构模型，基于渗流场-应力场耦合理论进行有限元数值模拟，并与监测结果对比分析，依此探讨RIS时空演化与库水加卸载及渗透过程的动态响应机制这一重要的科学问题将是一条有效的途径.

本文以紫坪铺水库为例，在对库区地质构造及水文地质结构条件、水库蓄水后地震活动时空演化特征等进行详细研究的基础上，利用有限元方法计算了水库蓄水过程中弹性附加应力场、有效附加应力场、孔隙压力和断层稳定性的动态变化，讨论了RIS时空演化与库水加卸载及渗透过程的动态响应关系.本文的成果不仅有助于丰富构造流体动力学的理论体系和拓宽RIS孕育发生机理研究的思路，同时对于区域稳定性评价及RIS预测也有重要的指导意义.

在库水加卸载及渗透过程中，地下岩石介质变形与流体渗流耦合作用造成的岩体失稳问题即为RIS.由于RIS的震源深度一般在地下10km以内^[15]，在此深度范围内，岩石体多处于上地壳脆性变形域，可在弹性力学范畴内讨论其应力与变形响应问题.为方便模型的求解，并考虑到研究问题的尺度与表征单元体（Representative Elementary Volume，REV）的关系，将岩石介质简化为多孔连续介质.同时假设：岩石介质变形为稳态的小位移变形；孔隙岩石介质中的流体运动符合Darcy渗流定律.按照弹性力学的约定，拉应力为正.

（1）连续性方程

由质量守恒定律可导出岩石骨架的连续性方程为

(1)

式中，v_s为岩石骨架的运动速度，ρ_s为岩石颗粒密度，φ为孔隙度.

若假设岩石颗粒的密度不变，在岩石骨架为稳态小位移变形的条件下，（1）式可进一步简化为^[16]

(2)

式中，ε_kk为体应变，k=1，2，3.

（2）本构方程

本构方程用张量形式表示为^[17]

(3)

式中，σ′_ij为有效应力张量，ε_ij为应变张量，D_ijkl为弹性系数张量.

（3）平衡方程

基于连续介质力学理论和多孔介质有效应力定律，当忽略惯性力作用时，岩石骨架变形平衡方程可用位移的形式表示为^{[8, 18~20]}

(4)

式中，u_i为岩石骨架位移分量，G为剪切模量，λ为拉梅（Lame）弹性常数，δ_ij为Kroneker符号，α为Biot系数，P为孔隙流体压力，f_i为体积力分量.

（1）连续性方程

对于饱和流体渗流，若忽略源汇相影响，则流体渗流连续性方程为

(5)

式中，ρ_w为流体的密度，v_s为岩石骨架的绝对运动速度，v_rw为流体相对于骨架的运动速度.

（2）流体渗流方程

饱和流体的运动方程为

(6)

式中，k为岩石的等效渗透率，μ为流体的动力黏度，φP为流体的压力梯度，g为重力加速度，φD=(0, 0, 1)为流体位势梯度，将（6）式代入（5）式，便可得到以流体位势表示的渗流方程^[21]

(7)

式中，β_p为流体压缩系数.

根据Mohr-Coulomb强度理论，若作用于断层面上的剪切应力为τ，则τ越接近于断面的抗剪强度τ_n（τ_n=c-σ_ntgφ），断层就越容易发生破裂失稳.可以定义描述断面趋于破裂失稳程度的库仑破裂应力（Coulomb Failure Stress，CFS）为^{[22, 23]}

(8)

式中，|τ|为断层面上的切向剪应力，σ_n为断层面上的法向正应力，K为断层静摩擦系数（K=tgφ，φ为内摩擦角），c为内聚力，P为流体压力.若c和K不随时间变化，则CFS的变化为

(9)

式中，Δτ为有效应力变化张量在断层滑动方向上的投影分量，若与滑动方向一致取正，反之取负，Δσ′_n为有效应力变化张量在断层面法向上的投影分量，ΔP为流体压力的变化量.

紫坪铺水库是岷江上的一座大型水利枢纽，位于四川省成都都江堰市境内岷江上游的麻溪乡，距成都市60余公里，是一座以灌溉和供水为主，兼有发电、防洪、环境保护、旅游等综合效益的大型水利枢纽工程.该水库混凝土面板堆石坝最大坝高156m，为国内仅有的几座同类型高坝之一.

在大地构造上，紫坪铺水库位于青藏高原东缘的龙门山造山带的中段.库区及邻近区域（103.2°N~103.9°N；30.7°E~31.3°E，本文研究范围）被主干断裂分隔，可划分为茂汶韧性剪切带、中央推覆构造带、龙门山前缘拆离带、前陆扩展变形带和川西前陆盆地5个地质构造单元，它们在物质组成、构造层次和变形样式等方面具有明显的差异^[24~30].紫坪铺库体座落于龙门山前缘拆离带内，邻近中央推覆构造带和前陆扩展变形带（图 1）.

图 1

Fig. 1

图 1 紫坪铺库区地质构造略图(据灌县幅1:20万区域地质图修改[33]) 1.主干断裂; 2.次级断裂; 3.逆冲断裂; 4.飞来峰; 5.紫坪铺水库坝址; 6.水系; Ⅰ.茂汶韧性剪切带; Ⅱ.中央推覆构造带; Ⅲ.龙门山前缘拆离带; Ⅳ.前陆扩展变形带; Ⅴ.川西前陆盆地 Fig. 1 Structural outline map of the Zipingpu reservior and its adjacent areas (Modified from 1:200000 areal geological map of the Guanxian[33]) 1.Main fault; 2.Sub-fault; 3.Reverse thrust; 4.Outlier of overthrust mass; 5.Zipingpu Reservior Dam; 6.Water system; Ⅰ.Mao-Wen ductile shear belt; Ⅱ.Central napped structure belt; Ⅲ.Longmenshan frontal detachment belt; Ⅳ.Expanded foreland distortion belt; Ⅴ.Chuanxi foreland basin

中央推覆构造带又称后山带，主要由前震旦纪基底杂岩和震旦系-泥盆系浅变质岩组成.该带以基底卷入冲断变形和出露彭灌杂岩体为典型特征.杂岩体是伴随着龙门山造山带的逆冲推覆作用而形成的构造侵位体^{[31, 32]}，其内部发育着一系列走向大致NE、倾向NW的近叠瓦状断裂系统.龙门山前缘拆离带又称前山带，主要由泥盆系-中下三叠统碳酸盐岩和上三叠统须家河组碎屑岩组成.该带的推覆构造体主要发育在上三叠统小塘子组、须家河组砂泥岩和含煤系地层中，表现为一系列逆冲推覆片体呈叠瓦状排列，向深部汇聚到通济场断裂主滑脱面上.该带另外一个显著的特点是彭灌飞来峰群的发育^{[24, 26~28]}，紫坪铺水库大坝就建在懒板凳-白石飞来峰之上.前陆扩展变形带又称山前带，是四川盆地在逆冲作用向前陆扩展过程中卷入变形的部分，主要由侏罗系-古近系陆相红层组成，构造变形较弱，主要为岩层产状向盆地内部迅速变缓的单斜构造（图 2）.

图 2

Fig. 2

图 2 紫坪铺库区二维地质构造剖面(位置见图 1) 蓝色虚线框部分据川S-6地震解释剖面[30]; 绿色虚线框部分据CX-NW-03-167和CX-NW-05-26地震解释剖面（中石化内部资料）；紫色虚线框部分据Inlrne740解释剖面（中石化内部资料）深部据松潘-中江大地电磁测深二维反演结果[34]. Fig. 2 The two-dimension geologic structure profile across the Zipingpu reservior and its adjacent areas (location in Fig. 1) Section framed by the blue broken line s based on Chuan-S-6 seismic interpretion profile[30]; Section framed by the green broken line s based on CX-NW-03-167 and CX-NW-05-26 seismic mterpretion profiles (restricted data from SINOPEC); Section framed by the purple broken line s based on Inline740 seismic interpretion profile (restricted data from SINOPEC); Deep section s based on the MT profile from Songpan to Zhongjiang[34].

根据区域构造沉积演化历史、地层沉积环境及地层发育特征等，可将研究区地层划分为基底和盖层两大部分^[35].基底地层形成于中-新元古代时期的晋宁运动和澄江运动，在龙门山地区称为黄水河群，主要为一套海底喷发火山岩建造和浅海泥页岩、碎屑岩建造.盖层主要发育两大地层序列，即震旦系-中三叠统海相碎屑岩及碳酸盐岩序列和上三叠统-第四系陆相碎屑岩序列（表 1）.

表 1(Table 1)

表 1 研究区地层序列及岩体渗透稳定性简表(地层序列部分据文献[35], 有修改) Table 1 Stratigraphic sequence and rockmass seepage stability of the study area (Stratigraphic sequence modified from Ref.[35])

表 1 研究区地层序列及岩体渗透稳定性简表(地层序列部分据文献[35], 有修改) Table 1 Stratigraphic sequence and rockmass seepage stability of the study area (Stratigraphic sequence modified from Ref.[35])

在RIS评价和预测研究中，明确库区的岩体结构特征及在渗流作用下岩体的稳定性问题至关重要.为了对岩体特性作进一步概化，建立起定量评价的地质模型和力学模型，依据谷德振等对岩体结构的分类标准^{[14, 36]}，将研究区的岩体结构划分为整体结构、层状结构和散体结构三种类型，其中，整体结构和层状结构分类见表 1，散体结构主要分布于主干断裂碎裂岩带.岩体渗透稳定性的分类，暂依据RIS震例关于震中岩体性质的统计结果^[14]，分为高、中、低三类（表 1）.

龙门山地区岩石圈结构具层圈性，最重要的滑脱界面有15~20km深处的低速层、寒武系泥质粉砂岩层、志留系千枚岩和泥页岩层以及中下三叠统富膏盐岩层.这一构造特征决定了龙门山冲断带深部的多级滑脱和浅层的多层次推覆^{[34, 37~40]}.

在地表露头观察，研究区内发育有茂县-汶川断裂（WMF）、北川-映秀断裂（BYF）、通济场断裂（TJCF）、安县-灌县断裂（AGF）和广元-大邑断裂（GDF）5条主干断裂（图 1）.自北西向南东，主干断裂发育时间由早到晚，变形行为韧性递减，脆性增强，反映了其形成深度也在逐渐变浅.

通过对过龙门山中段前缘地带的CX-NW-03-167、CX-NW-05-26和Inline740等石油地震剖面与钻井资料的联合分析表明，北川-映秀断裂以东的沉积盖层构造在垂向上被三叠系雷口坡组和嘉陵江组主滑脱面隔开，其下属于印支期前沉积的地层，主要发育早期的正断层；后者为印支期间及以后沉积的地层，构造变形样式更为复杂，除发育有一系列倾向NW的叠瓦状逆冲推覆构造外，还发育有双重构造、反冲断层、三角带、滑脱构造等多种构造样式（图 2）.

断裂的渗透结构对孔隙压力变化下断裂的力学响应具有重要的影响.Jonathan S C等^[41]在对走滑断裂渗透性研究的基础上，将断裂渗透结构划分为隔水的断层核和导水的断层破碎带两部分，并依据断层核和断层破碎带的性质及其相互关系，划分为局部导水、局部阻水、散状导水和复合型4种类型.实际的断裂渗透结构情况要复杂得多，并非此四类就能概括.断裂的力学性质、运动性质、变形行为特征、新活动性及断裂与围岩的组合形式等因素都对渗透结构有着不同程度的影响.

研究区域现今处于以挤压为主的构造应力环境，主干断裂均不同程度地表现出走滑逆断性质.近年来，国内外学者对分布在褶皱冲断带内的滑脱构造层的研究发现，在构造挤压环境中，往往滑脱层之上的地层变形强烈，滑脱层之下地层变形相对较弱^{[40, 42~45]}.根据地表露头调查，研究区内北川-映秀、通济场和安县-灌县三大主干断裂的上盘地层构造变形相对强烈，影响带构造岩劈理化、碎裂化程度高，胶结差，羽状裂隙发育.在岩性、地貌有利地带出露大泉，表明它们具有较好的渗透条件^[46].

此外，断裂的新构造活动性是控制其渗透性能的另一个重要因素.李明诚等^{[47, 48]}认为，无论是张性正断层，还是压性逆断层，只要它们处在活动中，就具有通道性.根据构造地貌、线性影像、地质剖面揭露的断裂活动性以及所获得的年代学证据^{[49, 50]}，研究区内的茂县-汶川断裂、通济场断裂和安县-灌县断裂均为晚更新世活动断裂.北川-映秀断裂是2008年汶川M_S8.0级地震的主要发震断裂之一，为全新世活动断裂.广元-大邑断裂至今尚缺少新年代学的直接证据，在位于断裂下盘的花园场竹瓦铺一带，第四系厚达500m，而位于断层上盘的相应地层厚度尚不足100m.从这一现象推测，该断裂在较近时期可能仍在持续活动.

基于以上分析认为，研究区内的5大主干断裂都不同程度地具有使地表水体向深部渗流的通道性，但其渗透结构并不同于Jonathan S C分类，是一种断裂上盘破碎带导水、下盘地层及断层核阻水的新的渗透结构类型，我们将其称为“下阻上导型”.主干断裂导水性能的差异主要是由断裂与围岩组合形式的不同造成的.横跨研究区从北西到南东，无论是主干断裂还是围岩其力学性质都是脆性逐渐增强，反映在导水性上是影响带逐渐变宽，渗透性能逐渐变好.

研究区内发育的NW向张扭性断裂规模小，延伸不远，破碎带窄，多出现在碳酸盐岩飞来峰内和飞来峰底座的三叠系须家河组煤系地层中.这些小断裂属于Jonathan S C分类中的散状导水类型，导水性好，往往切割了NE向的压性断裂，起到了沟通不同含水层的作用.

历史地震记载，研究区惟一一次4级以上地震是1787年12月发生在都江堰附近的4³/4级地震，考虑到中国西部地区的地震监测台网普遍建立较晚，估计该区应有不少4级以上地震被漏记.自1970年四川及邻区测震台网正式观测以来，研究区M_L2.5级以上地震记录相对完整、可靠^[51]，可作为地震活动性分析的起始震级.1970年1月1日至2004年11月30日期间，该区共发生M_L2.5~2.9级地震94次，M_L3.0~3.9级地震34次，M_L4.0级以上地震3次，分别是1983年3月19日漩口M_L4.0级地震，1986年11月10日映秀M_L4.1级地震和1993年12月30聚源M_L4.4级地震.在这131个M_L≥2.5级的地震中，仅75个有震源深度定位数据，对其进行统计分析显示，60%以上的震源深度分布在6~15km范围内，最深可达34km.可能是由于可靠地震资料的时间较短，该时段震中空间分布非常分散，也未明显表现出地震频度和强度随时间的周期性变化规律.

为了便于分析水库蓄水后地震活动的时空演化规律，我们利用双差定位法对紫坪铺水库台网（包括7个子台和1个中继站）在2004年8月16日至2008年5月10日期间记录到的1569个小震进行了精定位.精定位后，到时残差的均方根平均值为0.048s，震源位置的估算误差在水平方向上平均为0.63km，在垂直方向上平均为0.51km，得到了研究区较为精确的地震参数^[52].

相比于水库蓄水之前地震震中分布比较分散，在2004年8月份坝前水位有明显升高以后，至2008年4月30日期间，地震活动图像上明显呈现出条带状分布和丛集分布的特点（图 3a）.震中大多分布于库岸周边10km以内的范围（Ⅳ区除外），小震活动平行于主构造线呈NE向优势分布，在这一总体特征下，库体覆盖区（Ⅱ区）的地震活动却非常少，几乎为空白.此外，小震活动密集分布在3个小区域：第1个是图 3a中的Ⅰ区，该区位于水库西南侧的龙门山前山带，其东南和西北两侧分别被岷江的支流水系寿溪河和古溪沟所夹持；第2个是图 3a中的Ⅲ区，该区位于水库东北侧的龙门山前山带和北川-映秀断裂地表露头带，邻近岷江的支流水系龙溪河；第3个是图 3a中的Ⅳ区，该区位于坝址下游9~18km的都江堰市幸福乡、中兴镇和聚源镇一带，横跨广元-大邑隐伏断裂.在这3个地震丛集区中，小震分布最为密集的是水库西南侧的Ⅰ区.

图 3

Fig. 3

图 3 紫坪铺库区不同蓄水时段地震震中分布图 浅灰色区域代表飞来峰构造；浅蓝色区域代表此时段受水库回水影响库段的大致位置. Fig. 3 Earthquake epicenter distribution in Zipingpu reservior and adjacent areas during different water impounding time intervals Light grey area express outlier of overthrust mass; Light blue area express reservoir segment being influenced by banked-up water in that time interval.

为了进一步分析在空间分布上地震活动随水库蓄、放水过程演化的规律，依据2004年1月1日至2008年4月30日期间的坝前水位变化资料，并考虑到精定位数据的时间范围，将库水加卸载过程划分为5个不同的时段，即：2004年8月16日~2005年9月29日为低幅加卸载期；2005年9月30日~2006年11月30日为第1次高幅、快速加载期；2006年12月1日~2007年4月15日为第1次高幅、快速卸载期；2007年4月16日~2007年12月15日为第2次高幅、快速加载期；2007年12月16日~2008年4月30日为第2次高幅、快速卸载期（图 4）.同时，利用精定位数据绘制了不同时段的震中分布图（图 3b~3f）.

图 4

Fig. 4

图 4 不同库水加卸载阶段的划分 Fig. 4 The classification of different period of reservoir water load-unload process

在2004年8月开始的低幅加卸载时段（阶段Ⅰ），坝前水位从海拔高度约754.5m上升到约773.5m后又有回落.据库区海拔可知，此时水体大致分布于图 3b中浅蓝色区域所覆盖的库段.该时段小震活动主要集中在库体的北西侧、北侧和北东侧距离库岸约2km以外、10km以内的区域，沿“彭灌”飞来峰构造和前山断裂破碎带呈NE向条带状分布，而在水库的东南侧及库体覆盖区则很少有地震发生.此外，在距离坝址约25km处的红路附近发生的小震群应是本底自然地震活动，与水库蓄水无关.

2005年9月30日起的阶段Ⅱ，坝前水位从海拔高度约759.5m迅速升高到820m以上，至2006年秋汛末首次抬升至正常蓄水位875m，并持续近2个月.在水库达到875m正常蓄水位时，回水刚好到达岷江与北川-映秀断裂地表露头带的交界处，此时水体大致分布于图 3c中浅蓝色区域所覆盖的范围.随着水库水位的不断升高，回水影响范围的逐步扩大，该时段地震活动表现出了三个显著的特点：一是Ⅲ区的地震活动沿袭了上一时段小震条带的展布方向，向远离库岸的虹口附近迁移；二是小震活动向水库西南侧的Ⅰ区迁移；三是Ⅳ区开始发生了一些响应性的小震活动.尽管此时库水位已有大幅抬升，库体影响范围不断扩大，但在库体覆盖区（Ⅱ区）及水库东南侧的近库岸区地震数目却依然是寥寥无几.

在第1次高幅、快速卸载时段（阶段Ⅲ），坝前水位从海拔高度约874.5m下降到约819m.据库区海拔可知，当坝前水位为820m时，回水可影响到图 3d中浅蓝色区域所示的库段.因此，此时段受水库回水影响的库体范围从图 3c中的浅蓝色区域逐渐收缩到图 3d中浅蓝色区域.相比于前两个时段，该时段除在Ⅰ区小震分布相对集中外，地震活动频度和强度都明显的降低.Ⅰ区的地震活动表现出向寿溪河和古溪沟所夹持的中间部位收缩的趋势.

在第2次高幅、快速加载时段（阶段Ⅳ），库水位从海拔高度约819m上升到约873m，水库水体上升到图 3e中浅蓝色区域所覆盖的库段.该时段小震活动主要集中在水库西南侧的Ⅰ区，Ⅲ区的地震活动明显减弱.

在第2次高幅、快速卸载时段（阶段Ⅴ），库水位从海拔高度约873m下降到约821.5m，回水影响情况与阶段Ⅲ相似.该时段地震活动主要集中在Ⅰ区和Ⅳ区，在水库快速蓄水以来时隔近29个月后，在Ⅳ区发生了都江堰小震群活动（最大震级为M_L3.7级）.

水库蓄水后的震源深度统计结果显示，研究区范围内约占70%的小震震源深度分布在4~10km，10~25km也有少量地震发生，并且震源深度分布有一定的连续性.相对比水库蓄水前震源深度分布的特征，蓄水后小震震源分布似乎显示出从深部向浅部迁移的趋势，这一变化可能与水库蓄水有关.为便于寻找震源深度分布与地下断裂构造展布的关系，我们横跨几个地震丛集区分别作了震源深度剖面A-A′、B-B′、C-C′、D-D′（剖面位置见图 3a），将剖面两侧各4km范围内的震源位置垂直投影到剖面上，同时，将主要断裂构造也投影到剖面上（图 5）.

图 5

Fig. 5

图 5 震源深度分布剖面图 1.水库位置；2.地震；3.逆断层；4.正断层；5.滑脱层；6.推测 Fig. 5 The profile map of earthquake depth 1. Reservoir location; 2. Earthquake; 3. Thrust fault; 4. Normal fault; 5. Decollement zone; 6. Hypothetical fault

A-A′剖面和B-B′剖面分别垂直于主要构造线方向横跨了Ⅰ区和Ⅲ区，在这两个剖面上表现出了较为一致的震源分布特征，即：震源深度集中在地下4~10km范围内，10km以下虽也有地震发生，但其比例不超过总数的8%；沿通济场断裂在地腹的延伸方向，震源相对密集地分布在与断裂邻近的上、下盘地层中，而在北川-映秀断裂和安县-灌县断裂的陡倾段则很少有地震发生；在邻近通济场断裂和安县-灌县断裂的深部汇聚处，由这两条断裂和反冲断层构成的三角带震源分布最为密集；在通济场断裂和安县-灌县断裂地表露头带之间向下约6km以上的区域几乎没有地震发生.

C-C′剖面平行于龙门山前山带主构造线方向，大致跨越了库体覆盖区及Ⅰ、Ⅲ地震丛集区.该剖面中左、右两侧的震源深度集中分布在4~10km范围内，与A-A′剖面和B-B′剖面相符.在第1次低幅加卸载时段，Ⅲ区中南段的震源深度主要分布在8km左右；而在随后的第1次高幅、快速加载时段，中南段地震活动减弱，且北段出现了震源深度集中在5km以内的浅层地震活动.此外，剖面中段的库体覆盖区以下区域几乎为地震空白区.

D-D′剖面垂直于主要构造线方向，横跨了龙门山后山带、前山带、山前带和川西前陆盆地4个地质构造单元，并且过Ⅳ区地震丛集区.该剖面中通济场断裂以西部分除地震发生个数较少外，其震源深度分布特征与A-A′剖面和B-B′剖面类似.剖面右侧的震源分布主要是2008年2~4月份都江堰幸福乡、中兴镇和聚源镇一带的小震群活动所致，震源深度主要集中在5~16km范围内，最深已达25km.从该震群震源深度分布图像与断裂构造的位置关系来看，小震群发生与广元-大邑断裂和原地断层的活动有关.

为了分析库区地震活动对库水加卸载过程的时间响应规律，我们针对以上3个地震丛集区，分别绘制了水库坝前水位变化与各分区地震发生频次的关系图（图 6）.

图 6

Fig. 6

图 6 紫坪铺水库水位变化与各分区地震频次关系图 (a)Ⅰ区; (b)Ⅲ区; (c)Ⅳ区.曲线表示水位线; 方柱表示地震每月发生的频次. Fig. 6 The relationship between the change of water level and rate of earthquake occurrence in different areas (a) Area Ⅰ; (b) Area Ⅲ; (c) Area Ⅳ. Curves express water level; Square columns express the number of earthquake occurrence per month.

在2005年9月30日之前，由于水库回水尚未影响到邻近Ⅰ区的库段，该区的小震活动基本代表了背景性的自然地震活动，地震数目较少.随着库水位的快速抬升，库水影响范围的逐渐扩大，该区地震数目明显地增多，表现出了伴随着库水位的高低变化，小震发生频度呈同步变化的趋势.在2006年9~10月和2007年9~12月两个蓄水高峰期，地震活动频度亦达到了峰值（图 6a）.Ⅰ区的小震活动显然与水库蓄水有关，按照SimpsonDW等的分类^[53]，可能属于“快速响应型”水库诱发地震.

Ⅲ区邻近坝址区，受坝前水位变化的影响较大.该区在库水低幅加卸载的初期就迅速发生了响应性的小震活动，出现了一个高峰期.在第1次高幅、快速加载时段，伴随着水位的快速抬升，该区又出现了一个小震活动的高峰期，此后便逐渐地减弱（图 6b）.

在水库快速蓄水之前，Ⅳ区小震发生的数目很少.2005年9月30日之后，伴随着水位的快速抬升，该区小震发生的频度似有微弱的抬升，到2008年2~4月期间的水库卸载阶段，忽然发生了都江堰小震群活动（图 6c）.如果该区的小震活动与水库蓄水有关的话，依据SimpsonDW等的分类，可能属于“滞后响应型”水库诱发地震.此方面在后文中将做进一步的讨论.

有限单元法在对复杂地下介质域和边界条件的处理方面，以其灵活的几何单元划分、边界条件处理及对不同性质介质的处理技术见长.周斌根据Galerkin有限元理论，由岩石介质变形与流体渗流耦合数学模型控制性方程，推导出“弱积分”形式的Galerkin有限元公式，并利用有限元程序自动生成软件FEPG编写了计算程序^[21].此程序可离散求解流体压力场、位移场、有效应力场和应变场的分布.本文重点计算了在库体荷载作用下产生的弹性附加应力场，以及在库体荷载和附加水头压力扩散耦合作用下产生的有效附加应力场，并在不同断层上分别设置了观察点（图 7），观察水库动态加卸载过程中的孔隙压力响应及断层稳定性变化.

图 7

Fig. 7

图 7 有限单元划分及观察点位置示意图 Fig. 7 The sketch map of finite element division and observation sites

考察水体载荷造成地下岩体内应力的变化，本应考虑实际的地形及地质构造情况建立三维模型.但由于受资料的限制，本文将其简化成平面应变问题，选用图 2所示的二维地质剖面作为有限元数值模拟实验的地质模型，进行了有限单元划分，共12513个单元，6434个节点（图 7）.

在确定地质体力学参数与物性参数时，主要参考了中石化南方勘探开发分公司在龙门山地区及四川盆地开展的岩石力学参数与物性参数测试成果，以及岩石力学手册中研究区附近的映秀湾水电站和四川“红层”地层的实测数据^[54]，同时考虑到岩体与岩石之间力学性质和渗透性的差异，对参数进行了适当的调整（表 2）.

表 2(Table 2)

表 2 地质体力学参数与物性参数 Table 2 The mechanic parameters and physical properties of geologic bodies

表 2 地质体力学参数与物性参数 Table 2 The mechanic parameters and physical properties of geologic bodies

依据2005年3月31日~2008年4月14日期间坝前水位变化资料，将相对于2005年3月31日基准水位（海拔高度约754.5m）的水位变化用分段直线加载函数在库体处施加动水头边界和动荷载边界.模型的底边界设置为水平滑动，两端点固定.计算时共设置了38个计算步，每步历时30天.

在水库蓄水过程中，地表水体荷载的作用会在库底岩体内产生弹性附加应力的响应，从而影响到断层稳定性的变化.在不同的构造应力环境中，弹性附加应力对断层稳定性所起的作用截然不同.在拉张构造环境中，最大主应力方向近于垂直，最小主应力方向近于水平，弹性附加应力与初始构造应力场叠加的结果，会使莫尔圆半径变大，增加逼近外包络线的机会，促使断层破裂趋势增大；在剪切构造环境中，最大和最小主应力方向都近水平，中间主应力方向垂直，弹性附加应力对断层的稳定性几乎无影响；在挤压构造环境中，最大主应力方向近于水平，最小主应力方向近于垂直，弹性附加应力的叠加会使莫尔圆的半径变小，并远离外包络线，对断层破裂起到一定的抑制作用.地表水体荷载的作用除了会产生弹性附加应力响应以外，在水库快速加载时，岩石的压缩变形还会造成孔隙压的瞬时升高，随着时间的推移、孔隙水的排出，这一响应将逐渐被消弱.

紫坪铺库区处于以挤压为主的构造应力环境中，单纯的库体荷载作用会有利于库底断层的稳定.由不同计算时间步的弹性附加应力场图像可以看出，弹性效应主要位于库底之下的高附加应力区.图 8为水库蓄水接近历史最高水位时的弹性附加应力场云图（此时水位约874.5m，相对基准水位面升高了120m）.结果显示，由库底向远处，弹性附加应力逐渐衰减，垂向分量（σ_yy）向深部的传播明显受陡立断层结构面的限制，集中分布在库底正下方的区域（图 8a）.不同观测点断层库仑应力动态变化曲线（图 9）亦表明，距离库底越远，库体荷载对断层稳定的影响程度就越低，在距离较远的汶川-茂县断裂上的S1观察点上，ΔCFS的最大值仅为-0.00066MPa，而在较近的S5观测点上，ΔCFS可达-0.3MPa.局部构造条件对弹性附加应力场分布也有一定程度的影响，图 8b中黑色虚线所围限的区域为安县-灌县断裂拐点与反向冲断层交汇的部位，此构造部位在垂向附加应力为较高的压应力状态，且剪应力分量数值较低的情况下，水平附加应力却出现了较高的张应力状态.在此情况下，附加应力场与初始构造应力场叠加的结果会促使该区域更趋稳定.在震源深度分布剖面上该区域几乎为地震空白区（图 5），就进一步证实了库体荷载对该区稳定起到了主导性的作用.从附加剪应力分布图像（图 8c）可以看出，在地下6km以上的区域，以库底向深部延伸的垂线为界，左侧的剪应力为负值，右侧为正值，这将更有利于右侧NW倾向高角度逆断层的稳定，S1和S5两个观察点的断层库仑应力变化情况也证实了这一点（图 9）.

图 8

Fig. 8

图 8 相对水位变化120 m时的弹性附加应力场云图 Fig. 8 The additional elastic stress field induced by water level change of 120 m

图 9

Fig. 9

图 9 不同观察点孔隙压力及断层库仑应力变化曲线 a孔隙压力变化；b库体荷载与孔隙压力扩散耦合作用下断层库仑应力变化；c库体荷载作用下断层库仑应力变化. Fig. 9 The change of pore pressure and Coulomb failure stress at different sites a Pore pressure change; b Coulomb failure stress change affected by the coupling process both from reservoir body gravity and pore pressure diffusion; c Coulomb failure stress change affected by loading reservoir body gravity.

由于仅从地表观察无法准确判断深部岩体和断层的渗透性能，因而对于地表水体和深部流体之间联系的通道、强弱以及地表水体向深部循环所能影响的范围等问题一直是困扰着水文地质界的一个难题.近年来，一些学者对含油气盆地流体输导体系的研究发现，断裂构造带是地下流体运移的主要通道^[55]，其输导性能又与断裂渗透结构密切相关.由上文的分析可知，库区范围内的主干断裂可能属于“下阻上导型”的渗透结构，不同程度地具有地表水体向深部渗流的通道性.

我们对水库附加水头压力扩散的数值模拟结果表明，在水库水位变化过程中，附加水头压力主要沿与水库有水力联系的断裂带和裂隙发育的碳酸盐体进行传导，贯穿到库底和飞来峰底座的次级断裂成为库水下渗的优势通道.由水库蓄、放水过程中不同观察点的孔隙压动态变化曲线可以看出，随观测点埋藏深度的不同及与水库水力联系强弱的差异，孔隙压对附加水头压力变化的响应亦表现出不同的特征（图 9）.在与水库有直接水力联系的断裂带的浅部，如S3、S4和S8观察点，孔隙压响应迅速、变化幅度高；在有水力联系的断裂带的深部（S7点）或水力联系弱的断裂带上（S2、S6点），孔隙压响应并不明显，经过长时间以后，基本维持在0.2~0.3MPa的水平上不再升高；在与水库无明显水力联系的断裂带上（S1点），孔隙压几乎没有变化.各观察点孔隙压变化的另一个重要特征是滞后于附加水头的变化，距离库底越远，滞后现象表现得就越为突出.

相比于单纯在库体荷载作用下产生的弹性附加应力场分布图像，在库体荷载和水库附加水头压力扩散的共同作用下，有效附加应力场的水平向分量（σ′_xx）和垂向分量（σ′_yy）发生了明显的变化，而剪切分量（τ′_xy）变化不大.图 10是2005年3月31日坝前水位开始抬升起，60天、570天、990天和1140天后有效附加应力水平分量和垂直分量的变化云图.可以看出，随着水库蓄水时间的延长，孔隙压扩散对有效附加应力场的影响将逐渐扩大到地下10km以上的范围.在4个典型的区域（图 10中Ⅰ、Ⅱ、Ⅲ、Ⅳ），有效附加应力状态表现出不同的特征：（1）Ⅰ区位于库底正下方、通济场滑脱断裂面以上的区域，岩体破碎程度高，次级断裂发育且与水库有直接的水力联系.经过长时间以后，该区的水平向和垂向有效附加应力分量逐渐由压应力状态转变为张应力状态.到2008年4月14日，即坝前水位开始抬升后的1140天，水平向有效附加应力分量的平均值达到0.2MPa左右，垂向有效附加应力分量的平均值达到0.5MPa左右，由于此时附加剪切分量并无明显变化，这将促使该区逐渐由稳定状态向不稳定状态转变.经过长时间以后，该区可能会诱发一定数量的微震活动.（2）Ⅱ区位于安县-灌县断裂的上盘，发育有三叠系须家河组陆相碎屑岩和雷口坡组-嘉陵江组膏岩，岩性软弱，渗透条件差.该区除中部的反向逆冲断层受附加水头压力扩散影响较大外，其他部分的有效附加应力状态与单纯在库体荷载作用下的情况无明显差异.为考察其孔隙压响应及岩体稳定性，我们设置了SA观察点（图 7），并假设过该点的岩体中发育有一条倾向NW、倾角45°的断层.实验结果表明，该点的孔隙压响应极其微弱，最高值仅为0.034MPa，无论是仅在库体荷载作用下，还是在库体荷载与孔隙压扩散的共同作用下，假设断层的库仑应力均为幅度很大的负变化，且相差不大，进一步证实了该区岩体稳定性主要受控于库体荷载的作用（图 10）.在位于该区边界附近反冲断层上的S4观察点，尽管孔隙压响应值较高，但由于其邻近库底，在孔隙压力促使断层失稳与库体荷载抑制断层失稳两种不同作用的相互牵制下，ΔCFS几乎为零.（3）Ⅲ区位于安县-灌县断裂与通济场断裂在地腹的交汇区域，主要发育泥盆系-下三叠统飞仙关组海相碳酸盐岩，岩性坚硬，裂隙相对发育，渗透条件好.在2005年9月30日~2006年11月30日坝前水位高幅、快速上升时段，无论是有效附加应力场的水平向分量还是垂向分量都快速地从压应力状态转变为张应力状态，表明了附加水头压力扩散对该区有较大的影响.我们在该区设置了SB观察点（图 7），并同样假设过该点的岩体中发育有一条倾向NW、倾角45°的断层.实验结果表明，单纯在库体荷载作用下，假设断层的ΔCFS为负值，但由于孔隙压的快速升高，在水库快速蓄水90天后，ΔCFS转变为正值，并随时间的延续不断升高，最高达0.37MPa.在位于该区西缘通济场断裂上的S3和S8观测点，ΔCFS亦呈现出相似的变化.从震源深度分布剖面（图 5）可以看出，水库蓄水后的小震活动主要集中在这一区域，说明附加水头压力扩散作用对该区的诱发地震活动起到了主导性的作用.（4）Ⅳ区位于广元-大邑断裂和深部原地断层交汇的区域.由于该区距离库底较远，在库体荷载作用下产生的弹性附加应力场本身很微弱，有效附加应力场的变化主要受控于孔隙压扩散过程（图 10）.在S6观测点上，孔隙压持续升高，经过长时间以后基本维持在0.3MPa的水平上，受孔隙压变化的影响，ΔCFS亦不断增高，最高达0.1MPa左右（图 9）.

图 10

Fig. 10

图 10 有效附加应力场变化云图 (a)60天后；(b)570天后；(c)990天后；(d)1140天后. Fig. 10 The cloud map of additional effective stress field (a) 60 days later; (b) 570 days later; (c) 990 days later; (d) 1140 days later.

RIS的发生与库水加卸载及渗透过程中库底断层稳定性的变化密切相关，在库水加卸载及渗透过程中，引起断层库仑应力发生变化的原因主要有两个：一是由库体荷载作用产生的断层面法向应力变化Δσ_n和沿滑动方向的剪应力变化Δτ，这一力学响应是瞬时的，在不同的构造应力环境中对断层稳定性的贡献会有所差异；二是由地表水体向深部渗流造成的断层面上孔隙压的改变ΔP，这一变化对剪应力并无直接影响（暂不考虑水的弱化效应），但孔隙压的增加会降低断层面上的法向有效应力，从而降低了断层的抗剪强度，促使断层破裂失稳.同时，由于孔隙压力部分抵消了围压的影响，使岩体更易于产生张性破坏.从数值模拟结果可以看出，在以挤压为主的构造应力环境中，库体荷载作用的结果一般会使库底断层更趋向稳定，而附加水头压力扩散的效应则是促使断层趋向失稳.虽然地下岩石体对库体载荷的弹性应力响应是瞬时的，但在水库水位的动态变化过程中，这一响应也是一个动态变化的过程；同时，岩体孔隙压对水库附加水头变化的响应每时每刻也都在发生着变化.正是这个矛盾双方相互制约与平衡动态过程，控制了断层库仑应力变化的取向，从而决定了RIS时空演化的特征.这可能就是RIS时空演化对库水加卸载及渗透过程动态响应的力学机制.

从紫坪铺水库蓄水后小震活动的特征来看，亦表明了RIS时空演化受控于ΔCFS的动态变化过程.下面，就图 3中的1个稳定区和3个地震丛集区进行分别的讨论.在Ⅱ区，库体正下方6km以上范围的ΔCFS总体上在负值区间变化（图 10中S4、S5和SA观察点），库体荷载作用对岩体和断层的稳定起到了决定性的作用，因而该区小震活动非常之少.在距离水库较远的Ⅳ区，单纯库体荷载作用对断层稳定性的影响已非常微弱，ΔCFS主要受孔隙压变化的控制，在2008年2~4月期间，ΔCFS上升到0.1MPa左右（图 9中S6观察点），该区相继发生了都江堰震群活动，依此判断该震群活动可能与水库蓄水有一定的关系.水库东北侧的Ⅲ区远离库体覆盖区，库体荷载作用对断层稳定性的影响已不占据主导地位.由于数值实验选用的地质模型剖面穿过了库体覆盖区，而没有穿过Ⅲ区，所以对于Ⅲ区而言，数模结果肯定夸大了库体荷载作用的影响，其真实的情况会有所不同.同时由于Ⅲ区邻近坝址，受坝前水位变化的影响较大，加之，NW向张扭性断裂的沟通使该区深部的断裂与地表水体之间建立了良好的流体运移通道，所以即使在水库水位低幅变化时段，附加水头压力的变化亦可快速传导到该区的深部，引起断层库仑应力的变化，从而诱发小震活动，地震精定位分析结果已证实了这一点.在坝前水位快速抬升过程中，该区小震活动向远离库岸的虹口方向迁移的现象就进一步说明了孔隙压力扩散对该区诱发地震活动的促进作用.水库西南侧的Ⅰ区的情况与Ⅲ区类似，不再赘述.

断层上积累的应力在地震时部分释放而产生应力降，一次地震所产生的应力降一般为几兆帕到几十兆帕^[12]，相比而言，由水库蓄水造成断层库仑应力的变化量已微不足道.这也是许多学者强调的库区应变能积累已达到临界状态是诱发地震活动前提条件的主要依据.全世界已建大中型水库约有1万多座，但已诱发水库地震的只占水库总数的1%左右^[56]，在库区地质条件上，发震水库与未发震水库并未表现出显著的差别^[13]，这也可能主要是由于库区应变积累程度的差异造成的.

岩石渗流力学实验证实，脆性岩石在进入宏观破裂之前的扩容阶段，体积快速膨胀，内部产生大量的微裂隙，渗透性能得到有效地提高^[57~59].在诱发地震活动发生之前，孕震区岩石体可能已经进入扩容阶段，由扩容破碎作用产生的贯通微裂隙改善了脆性岩石体的孔渗条件，提高了断裂带的输导性能.此时，更利于地表水体向震源区的渗透，促使水库附加水头压力的快速扩散，从而诱发库区地震活动.同时，诱发地震（包括微小地震）的发生也会改善脆性岩体的渗透性，促进孔隙压的扩散.从另一个角度来讲，如果孕震区脆性岩石体的应变能积累程度普遍较低，岩石体尚处于弹性变形阶段，则其渗透性不会随围压的升高而得到有效地改善，甚至会有所降低^[59]，从而制约了水库附加水头压力的扩散，不利于诱发库区地震活动.对塑性岩石体而言，其渗透性随应变的增加而逐渐减小^[59]，孔隙压扩散效率降低，从而不利于诱发地震活动.因此，诱发地震活动在很大程度上依赖于岩性条件，同时，孔隙压扩散与诱发地震活动之间亦存在着相互促进与制约的辩证关系.

紫坪铺水库蓄水后，小震活动主要集中发生在脆性程度高、渗透稳定性低的碳酸盐岩地层中，而在岩性较软弱、渗透稳定性高的三叠系须家河组砂泥岩和煤系地层中很少有地震发生.除上文分析的原因外，由于在构造应力场作用下，脆性岩体内能聚集更高的变形能，原位地应力高于塑性岩体，在诱发地震发生之前，脆性岩体已经进入扩容状态，有效提高了岩石体的渗透性能，促使了诱发地震活动，这也可能是一个重要的诱震因素.

（1）在前人研究的基础上，本文运用比较构造学和解析构造学的理论和方法，对研究区深浅构造组合特征进行了解析，建立了研究区二维地质构造与水文地质结构模型.研究表明：紫坪铺库区位于青藏高原东缘龙门山造山带的中段，茂县-汶川、北川-映秀、通济场、安县-灌县和广元-大邑5条主干断裂控制了本区的基本构造格架.库区附近的主干断裂都不同程度地具有使地表水体向深部渗流的通道性，可能是一种断裂上盘破碎带导水、下盘地层及断层核阻水的“下阻上导型”的渗透结构类型，这种断裂渗透结构对孔隙压力变化下断裂的力学响应具有重要的影响.在地质构造单元上，紫坪铺库体座落于龙门山前缘拆离带内，该带在垂向上被三叠系雷口坡组和嘉陵江组主滑脱面隔开，上、下两部分的沉积地层组成及构造变形程度、样式等都存在着显著的差异.同时，依据谷德振等对岩体结构的分类标准，将研究区的岩体结构分为整体结构、层状结构和散体结构三种类型，岩体渗透稳定性分为高、中、低三个类别.岩体渗透稳定性差异在很大程度上导致了诱发地震活动对于岩性条件的依赖性.

（2）紫坪铺水库蓄水后，小震活动明显呈现出条带状分布和丛集分布的特点.震中大多分布于库岸周边10km范围以内且平行于主构造线呈NE向优势分布.小震活动在水库西南侧、东北侧和坝址下游9~18km的都江堰市幸福乡、中兴镇和聚源镇一带丛集分布，在库体覆盖区地震活动非常少，几乎为空白区.随水库蓄、放水过程的变化，小震活动亦表现出一定程度的震中迁移特征.水库蓄水后，小震震源深度优势分布在4~10km范围内，在通济场断裂与安县-灌县断裂的深部汇聚区域震源分布最为密集.同时，小震活动主要集中发生在脆性程度高、渗透稳定性低的碳酸盐岩地层中，而在岩性较软弱、渗透稳定性高的三叠系须家河组砂泥岩和煤系地层中很少有地震发生.在水库蓄水后地震活动的时间响应特征上，水库西南侧和东北侧两个丛集区的小震活动可能属于“快速响应型”RIS，而都江堰小震群活动可能属于“滞后响应型”诱发地震活动.

（3）RIS的发生不仅与库底先存断层等大型结构面有关，而且还受控于断层与围岩的组合形式，以及围岩体的岩石组合、岩性变化反映的岩体力学性质与渗透性能不均匀性.基于此认识，本文将RIS诱发机制的定量化模型分为2个层次：一是以孔隙介质为载体的流体渗流对岩体变形和稳定性的影响，由流-固耦合形式的岩体变形与孔隙渗流模型进行描述，此模型可反映岩体力学性质与渗透性能的不均匀性及其对库水加卸载的变形响应上的差异；二是对断层相关的RIS定量研究则将水库附加水头压力沿断层面（区）的扩散与断层库仑应力变化联系起来.两种形式模型方法的结合将能对RIS定量研究提供一个相对宏观的力学框架.

（4）RIS的发生与库水加卸载及渗透过程中库底岩体有效应力的变化密切相关，其先决条件是孕震区构造应变能的积累已接近临界程度，在库水作用下产生的弹性附加应力场、孔隙压力场和原地应力场的叠加场与岩体变形场之间的耦合作用下，致使岩体和断层面在特定的地段提前达到临界强度或破裂失稳条件，促使累积应变能以地震的形式释放出来.在以挤压为主的构造应力环境中，库体荷载作用的结果一般会使库底断层更趋向稳定，而水库附加水头压力扩散的效应则是促使断层趋向失稳，正是这个矛盾双方相互制约与平衡的动态过程，控制了断层库仑应力变化的取向，从而决定了RIS时空演化的特征.

由于紫坪铺水库距2008年5月12日汶川M_S8.0级大地震震中的距离大约十几公里，大地震发生后，有关紫坪铺水库蓄水是否会诱发汶川大地震的问题迅速引起了国内外广泛的关注.汶川大地震的成核区域在地下约19km深处的脆韧性转换带内，对于此深度范围内断裂渗透结构的变化、流体状态及其运移特征等，目前我们还知之甚少.另外，弹性力学方法已不适用于该深度范围的研究，这也是本文把所有数值模拟观察点都布置在10km以上的原因.对于发震活断层来说，脆韧性转换带对地震孕育发生具有重要的控制作用，特别是，地下流体在转换带的作用，断层摩擦滑动与脆性破裂转化与地下流体的耦合机制，是进一步研究中需要重视的科学问题.对于有限分布的水库载荷来说，二维有限元计算结果将会夸大水库的作用，尤其是水库重量的作用，在条件具备时，有必要基于三维有限元模型开展进一步的研究.

论文写作过程中，多次向中国地震局车用太研究员请教，得到了许多有益的启示；中石化研究院许华明博士提供了宝贵的第一手资料；四川省地震局程万正研究员、广东省地震局杨马陵研究员、中石化国际勘探项目部李行船博士、中国石油大学朱秀星博士都给予了热情的帮助，在此谨表衷心的感谢.同时，感谢三位匿名审稿专家提出了宝贵的修改意见.