水库蓄水与断层带流体孔隙压时－空分布的数值模拟——以紫坪铺水库为例

杨　彧杨晓松　段庆宝中国地震局地质研究所，地震动力学国家重点实验室　北京　1000292) HPT Laboratory，Department of Earth Sciences，Utrecht University，3584 CD Utrecht，The Netherlands

水库蓄水与断层带流体孔隙压时－空分布的数值模拟——以紫坪铺水库为例

杨彧1，2)杨晓松1) *段庆宝1)
1)中国地震局地质研究所，地震动力学国家重点实验室北京100029
2) HPT Laboratory，Department of Earth Sciences，Utrecht University，3584 CD Utrecht，The Netherlands

摘要通过对汶川地震破裂带上的断层岩及围岩的渗透率测量，获得了研究区岩石的渗透率随深度变化的规律。依据所得到的渗透率实验数据，数值模拟了紫坪铺水库蓄水所产生的流体孔隙压的时空分布。模拟结果显示，渗透率是否随深度变化，其流体孔隙压的模拟结果也相差很大。由于断层带具有高渗的特点，因此与水库底部存在水力联系的断层带及其宽度对流体渗透的模拟结果有明显的影响。断层带越宽，下渗作用越明显，流体孔隙压的分布与断层带的产状越接近。若断层带宽度为100～300m，2008年汶川MS8.0地震发生时，其震源区的流体孔隙压为0.1～0.15MPa，与库体载荷在断层面上产生的附加正应力相当。该模拟结果可以作为进一步分析汶川地震的发生是否与紫坪铺水库蓄水存在关联的重要约束。

关键词紫坪铺水库蓄水汶川地震流体渗透数值模拟流体孔隙压时空分布渗透率实验

〔基金项目〕地震动力学国家重点实验室自主研究课题(LED2010A03)、地震行业科研专项(201508018)和国家自然科学基金(41372202)共同资助。

0　引言

汶川地震是否由紫坪铺水库蓄水所诱发不仅是一个科学家所关注的问题(雷兴林等，2008; Ge et al.，2009;陈颙，2009;周斌等，2010; Deng et al.，2010; Gahalaut et al.，2010; Lei，2011;刘远征等，2014)，也已经成为社会各界热议的话题。西方媒体，诸如《The Telegraph》和《The New York Times》，《Science》新闻也热衷对此进行报道和评论。水库诱发地震是一类特殊的地震活动，主要指由水库蓄水或排水过程在一定时间段内在库区及其周边区域所引发的地震活动。在以往有关水库诱发地震的研究中，水库蓄水能够引发中小地震的认识已被普遍接受。但水库蓄水是否能够引起大地震一直存在争论(陈颙，2009)。汶川地震与紫坪铺水库存在关联的主要证据来自以下4个方面: 1)紫坪铺水库位于汶川地震震中东约10km; 2)在汶川地震发生之前2年多，紫坪铺水库开始蓄水，并且最大蓄水高度达到116m; 3)水库蓄水所产生的库仑应力增量与大地震应力降在同一量级; 4)水磨地震群的震源机制与汶川地震的相吻

合(胡先明等，2009)。地震史料统计表明，水库诱发的最大地震的震级还未曾超过6.5级(杨清源等，1996)。地表水能否渗透到10～20km的震源深度也存在疑问(Bell et al.，1978)。目前，数值模拟是研究水库蓄水对库区及周边地区几十千米深度范围内所造成的应力扰动的重要手段(Bell et al.，1978;雷兴林等，2008; Ge et al.，2009;周斌等，2010; Deng et al.，2010; Lei，2011)。考虑到水库蓄水引起的流体孔隙压对水扩散系数的敏感性，断层带和围岩渗透参数的选择合理与否，直接影响模拟结果，甚至有可能导致截然不同的结论。例如，Ge等(2009)假设断层带和其围岩的水扩散系数是均匀的，2维模拟计算得到震源区的流体孔隙压变化在0.01MPa量级;依据相同的假设，Deng等(2010)利用3维模型得到了更低的流体孔隙压变化，认为紫坪铺水库蓄水本身不可能诱发汶川地震。与此相反，Lei(2011)的模拟结果显示，紫坪铺水库蓄水导致的流体孔隙压扰动幅度比Ge等(2009)和Deng等(2010)的结果高1～2个数量级。纵观现有的模拟结果发现，模型参数的选择对模拟结果的影响往往被忽视。断层带和围岩通常被简化为均匀介质，而实际上其渗透率、孔隙度和弹性模量均随深度增加而发生了显著的变化。造成这种状况的主要原因是目前缺乏可靠的断层带水扩散系数随深度的变化实验数据。

水库蓄水导致的库区深部地壳的应力场扰动来自2个方面: 1)库体荷载引起的地壳变形所产生的附加应力场; 2)水渗透产生的流体孔隙压变化所诱发的应力场扰动。由于二者的时间响应速率不同，作为一级近似，可以将二者分别处理。计算前者时，将地壳考虑为等效弹性介质，其附加应力场与时间无关。当断层模型选定后，其附加应力场的计算相对简单，不同学者的计算结果出入不大。计算后者时，通常将库区地壳和断层视为孔隙介质。流体渗透引起的应力扰动依赖于流体孔隙压的时－空分布状态。龙门山地区处于以水平挤压为主的构造应力环境，汶川地震断层带具有走滑－逆冲的属性，由库体荷载引起的附加应力场倾向使断层的稳定性增加。因此，理论上讲，如果紫坪铺水库蓄水能够引起龙门山发震断层带稳定性下降，则源于流体孔隙压变化。为此，本文以渗透率实验结果为基础，通过模拟计算探讨断层带及其围岩的渗透率、断层带空间分布、断层带宽度等3个参数的变化对库区及周边地区流体孔隙压时－空分布产生的影响。

1　地质背景及地震时空分布

研究显示龙门山地区大致分为基底和盖层2个部分。基底地层主要为中、新元古代形成的海底喷发火山岩建造、浅海泥页岩和碎屑岩建造。盖层则是由三叠纪的海相碎屑岩及碳酸盐岩和陆相碎屑岩组成。龙门山断裂带由汶川－茂县断裂(后山断裂带)、映秀－北川断裂(中央断裂带)、安县－灌县断裂(前山断裂带)和山前隐伏断裂等组成。汶川地震过程中，沿映秀－北川断裂、安县－灌县断裂形成了2条同震地表主破裂带(徐锡伟等，2008)，其破裂长度分别约为280km和85km。其中前者在南段分叉为2条地表破裂。

紫坪铺水库是岷江上的一座大型水利枢纽，坐落于四川省都江堰市麻溪乡，距成都市大约60多千米，东北边可顺支流向北延伸到原龙池镇镇址，西南边距水磨镇约2km，西北则可达映秀镇下游。水库从2005年开始蓄水，水位由海拔759.5m升高至最高水位为876m，形成约116m的水头落差。2008年5月12日汶川地震时，紫坪铺水库的水位为828m(雷兴林等，2008)。详细研究表明，沿映秀－北川断裂带发生的地表破裂在都江堰段分为3支，其中一支

从水磨镇延至深溪沟(任俊杰等，2012)，地表破裂贯穿紫坪铺水库，地表垂直位移量达到0.7m，并造成横跨紫坪铺水库的庙子坪大桥第10跨桥身跌落和桥身位移。该地表破裂的位置对应于中摊铺断裂带。

在距水库约10km范围内，早在蓄水之前已经存在一些小地震活动，地震月频次＜25，最高震级＜3。蓄水之后，虽然在地震震级上没有出现明显的增大现象，但地震次数明显增多。在2008年2月达到了汶川8.0级地震前的最高频次(125次/月)。这些地震以震群形式出现。主要的几个地震震群分别为水磨震群、深溪沟震群以及都江堰震群。从震群发生的时间序列来看，蓄水后初期首先在水磨北发生了一系列地震，形成了紫坪铺水库西南侧的水磨震群，继而又在2个月之后形成了深溪沟震群(位于紫坪铺水库的东北边)，大约在2008年之后形成了都江堰震群(位于紫坪铺水库的东南边)。另外，在水库水位快速上涨或维持在高水位期间，部分地震表现出了对于水位变化的快速反应，同时也存在着对水位响应的滞后现象(马文涛等，2011)。

2　模型、参数及数值模拟方法

在MATLAB软件平台上，作者编写了用于计算流体扩散所引起的库区深部岩石和断层带流体孔隙压对于水库蓄水的响应的有限元程序。MATLAB软件平台可以方便灵活地划分计算网格，并能够对介质性质明显变化的位置实施网格加密处理。孔隙压时空分布的模拟计算涉及3方面内容: 1)介质模型的建立; 2)流体扩散方程及边界条件的确定; 3)计算参数值的选定。

图1 水库模型网格图Fig.1 Mesh model for numerical simulation.

2.1几何模型设置

为了计算水库蓄水造成的整个库区深部地壳的流体孔隙压扰动，本文建立了基于2维地质剖面有限元数值模拟的地质模型(图1)。模型以紫坪铺水库为中心，东西长60km，深30km。2维切面上共划分了57450个三角形单元，总节点数＞50000个。有限元网格的平均分辨率为40m。通过网格加密处理，使得断层带及附近几百米范围内的分辨率优于10m。基

本模型中包含了4条主要的断裂带，从西向东依次为汶川－茂县断裂、映秀－北川断裂，中滩铺断裂、安县－灌县断裂。地震测深资料显示，在研究区15～20km深度上，存在近水平的拆离断层带，上述4条断层向下收敛于该水平断裂带(徐锡伟等，2008;刘保金等，2008)。

模型的下边界和左右2个侧边界设定为Neumann边界(即在边界处，指定微分方程解的导数或偏导数之值) ;模型的上边界(包括水库两侧的上边界和水库的侧边界和底界)设定为Dilichlet边界(即在边界处，指定微分方程解的值)。计算时间步长的设定主要考虑模拟结果的精细程度和计算量。紫坪铺水库从第1次快速蓄水到汶川地震发生，历经约29个月，为能够精细地刻画孔隙压的时－空变化，时间步长设定为194.4h，约为8d。

2.2流体渗流和扩散方程

流体孔隙压由流体扩散方程描述(Biot，1941; Rice等，1976; Bell等，1978) :

式(1)、(2)中，P为孔隙压(Pa)，k为渗透率(m2)，η为水的黏性系数(Pa·s)，ν为泊松比，νu为非排水泊松比，B为Skempton系数。

2.3模拟参数

式(1)和式(2)中的参数中，除了渗透率外，其余的参数变化范围均不大。花岗岩的泊松比ν=0.26～0.27，石英砂岩的泊松比ν=0.25，页岩的泊松比ν=0.26～0.27，碳酸岩的泊松比ν=0.28(Christensen，1996)。非排水泊松比νu=0.29～0.31。η=10－3Pa·s。Skempton系数B= 0.7～1.0。参考紫坪铺水库库区的岩性，本文取ν=0.265;νu=0.30;η=10－3Pa·s，B=0.9。

岩石的渗透率变化幅度可达数个量级，因此合理确定渗透率参数对于准确模拟库区地壳的流体孔隙压至关重要。一般而言，结晶岩系的渗透率比沉积岩低得多;在沉积岩中，页岩和泥岩的渗透率远比砂岩低;而完整岩石的渗透性明显低于未胶结的断层岩，如断层角砾岩。介质的渗透性依赖于介质中孔隙度和孔隙的连通性。随着压力的增加，岩石中的孔隙被挤压而逐渐减小，尤其是联通孔隙的喉管可能被阻断，因此岩石的渗透性随着围压的增高(深度的增加)而显著下降，其变化可达1～4个数量级(20km范围)。为简便模拟计算又不失合理性，模拟计算中赋予围岩和断层带不同的渗透率值。围岩的渗透率依据对库区岩石渗透率的实验结果取值。断层带的渗透率则依据采自汶川地震破裂带上的断层岩的实验结果取值。

实验结果表明，断层岩和围岩的渗透率均随着深度增加而幂次减小(图2)。考虑到不同岩石的渗透率可存在较大变化，故在模拟中除了选用k0基本参数外，还分别用10k0和0.1k0作为地表渗透率(分别记为和)进行模拟计算，以便了解渗透率的不确定性对模拟结果产生的影响。具体的渗透率见表1。

由于断层带和围岩的渗透率相差悬殊，水库底部存在断层有可能对库区的渗透性结构产生重要的影响。因此在模拟中，需要对比研究库区底部存在断层与否对流体孔隙压分布的作用。断层带的渗透率较围岩可高出数个量级(随深度增加二者的差别也逐渐增大)，断层的宽度变化将会导致模型的渗透性结构变化，因此本文还考虑了断层带宽度对流体孔隙压结构的可能影响。地表观察显示，断层带宽度可＞100m。围陷波探测显示(Li et al.，2009)，中央断裂

带的断层宽度为170～230m。假设其余断层宽度不大于中央断裂带，模拟中选择100m和300m作为断层带宽度的上、下限是合理的并可以涵盖实际断层带宽度。

表1 断层岩及围岩的渗透率取值Table 1　 Permeabilities of fault and host rocks used in numerical simulation

紫坪铺水库于2005年9月30日快速蓄水，水位从海拔约750m迅速达到820m，并于2006年10月达到876m的最高水位。之后水库又经历了较大幅度的泄水—蓄水—泄水过程。至汶川地震发生时，其水位维持在821m左右。在近29个月的时间间隔中，平均水位约为840m。由于孔隙流体渗透缓慢，对水库小幅度的加载和卸载波动并不敏感，因此模拟过程中假设在29个月的时间内水库平均水位始终保持在840m，对应的水库的平均蓄水高度为90m。

2.4模型

图2 渗透率随深度变化图Fig.2 Depth dependence of permeabilities of fault rocks and its host ones.

表2 模型分类表Table 2 Category of models

如前所述，模型中变化幅度最大的参数为渗透率，而围岩与断层带的渗透率相差最为显著。因此本文着重讨论: 1)断层带和围岩渗透率随深度变化对流体孔隙压时－空分布的影响; 2)水库底部存在规模性断层与否对模拟结果的影响; 3)断层的宽度对模拟结果的影响。模型编号(表2)及其各模型的基本参数列于表3。

(1)模型A－C－E。围岩和断层带的渗透率均随深度的增加而幂次衰减，其所遵循的变化规律依据渗透率实验结果。断层带渗透率围岩渗透率水库底部存在规模性断裂带——中滩铺断裂。所有断裂带宽度为100m。为考虑渗透率变化带来的影响，计入和模型衍生出3个子模型，分别记为: A0－C－E、Amax－C－E和Amin－C－E，其中A0－C－E为基本子模型。

表3 模型编号及其相应的渗透率Table 3 Models and corresponding permeabilities

(2)模型A－D－E。围岩和断层带渗透率以及断层宽度与A0－C－E基本子模型相同，但水库底部及3km范围内不存在规模性断裂带。

(3)模型A－C－F。所有断裂带宽度为300m，其余参数与基本模型A0－C－E相同。

(4)模型B－C－E。围岩和断层带的渗透率均为常数，即在整个几何模型范围内，断层和围岩的渗透率不随深度发生变化。该模型可以划分为4个子模型，即B1－C－E、B2－C－E、B3－CE和B4－C－E。前3个子模型中的断层岩渗透率均为2.0×10－18m2(约相当于4km深度的断层带渗透率) ;围岩的渗透率分别为2.0×10－20m2、2.0×10－21m2和2.0×10－22m2(分别对应约2km、7km和20km深度的围岩渗透率)。作为对比，第4个子模型的断层岩和围岩渗透率分别为2.0×10－19m2和2.0×10－21m2。水库底部存在规模性断裂带——中滩铺断裂，其宽度为100m。

3　模拟结果

对于水库荷载之后孔隙压扩散的数值模拟结果揭示出了基本的渗透压时－空分布特征。流体主要沿着和水库有水力联系的断裂带渗透。伴随着水体的渗透过程，断层带及其附近的孔

隙压明显增高，而远离断层带的围岩以及与水库没有直接水力联系的断层带上的孔隙压力的变化十分微弱。此外，时间对孔隙压的演化有重要的影响。随着时间的推移，某一点的孔隙压会逐渐增加;孔隙压等值线沿着断层带由浅逐渐向深处延伸。

3.1时间对模拟结果的影响

图3显示了A0－C－E基础模型自水库开始蓄水至汶川地震发生约29个月里，库区孔隙压的时空变化。从渗透压分布的整体变化规律可以看出，水的渗透方向主要沿水库底部断层带向下渗透。渗透压等值线分布形态呈现出以垂直向下渗透为主，略偏向西侧的特点。该特征显示出水体的渗透与中摊铺和灌县－安县断裂带的产状密切相关，水的渗透方向完全取决于高渗透率区域的几何形态。远离水库底部的汶川－茂县断裂和映秀－北川断裂对水体的下渗影响甚微。图3－1～4分别显示了在水库蓄水3个月、6个月、12月和29个月(汶川地震发生前夕)时的孔隙压空间分布情况。模拟计算结果显示在汶川地震震源区，在上述4个时间点上，由水渗透而产生的孔隙压分别约为0.05MPa、0.07MPa、0.10MPa和0.12MPa。另外，在水库蓄水1a左右，水库下方的近水平的拆离断层带上的孔隙压能够达到0.1MPa。在整个水库蓄水期间，汶川－茂县断层上所承受的流体孔隙压均为最低，灌县－安县断裂带和中摊铺断裂带上承受的流体孔隙压最高，并且在不同的深度上，其流体孔隙压变化显著。例如中摊铺断裂带上的流体孔隙压随深度变化逐渐降低，自水库底部的0.9MPa降至15km深处的0.05MPa(3个月，图3－1)至0.12MPa(29个月，图3－4)。灌县－安县断裂带上所承受的流体孔隙压在10km左右的深度上达到最大;向浅部和深部其流体孔隙压均逐渐减小。映秀－北川断裂带上的流体孔隙压相对均匀，随着时间推移，其断裂带流体孔隙压自浅部向深部逐渐增加到约0.1MPa(图3)。

3.2围岩和断层渗透率对模拟结果的影响

对比基本模型(A0－C－E)，将围岩和断层带的渗透率同时增大和减小一个数量级(分别对应Amax－C－E和Amin－C－E子模型)来考察渗透率对模拟结果的影响程度。图4和图5分别为Amax－C－E和Amin－C－E模型在水库蓄水29个月时的模拟结果。与基本模型A0－C－E的模拟结果(图3－4)对比可以看出二者形态相似，但模型Amax－C－E的流体渗透压等值线分布向深部延伸更加明显，显示出水体更倾向于向下渗透(图4)。0.05MPa的等值线延伸到30km以下，但流体渗透压＞0.1MPa的区域则明显小于基本模型(A0－C－E)的模拟结果。与此相反，模型Amin－C－E的流体渗透压等值线分布近似球形，但略向西偏(图5)。其渗透压的分布状态(蓄水29个月)与模型A0－C－E在水库蓄水3个月时的渗透压分布状态(图3－1)接近，表明降低断层带和围岩的渗透率会显著地制约流体向下进行渗透。由此可见，渗透率的大小显著地影响水的渗透行为和流体渗透压的分布状态。

图4 Amax－C－E模型结果图Fig.4 The simulation result of model Amax－C－E.断层标号同图1

图5 Amin－C－E模型结果图Fig.5 The simulation result of model Amin－C－E.断层标号同图1

由于以往的模拟均假设库区岩石的渗透率为常数，因此有必要对比研究渗透率为常数与其随深度变化2种情况对库区深部地壳流体渗透压分布状态的影响。模型B1－C－E、B2－C－E 和B3－C－E分别描述了4种常数围岩渗透率所对应的流体渗透压的分布(表3)。均匀渗透率模型的基本特征是在不同深度上断层带和围岩的渗透率差别均相同，其流体孔隙压向断层两侧的扩展非常明显，尤其是在地壳浅部。当kf/kc为100时，其流体孔隙压分布近似半球形(图6)。当kf/kc为104时，其流体渗透压呈现出宽阔的带状分布，其产状与断层带一致(图8)。

3.3模型结构对模拟结果的影响

模型结构对于模拟结果的影响主要从两方面讨论。对比基本模型A0－C－E，本文考察了水库下方断层带宽度为300m(模型A0－C－F)和水库下方有无规模性断层带(模型A0－D－E) 2种情况。

模型A0－C－F中的断层带宽度为300m。其模拟结果与A0－C－E存在非常大的差异，主要表现在2个方面:其一，流体渗透压的分布受中摊铺和灌县－安县断裂带产状控制，渗透压等值线明显西偏，显示水体下渗主要受到具有高渗透特征的断层带控制;其二，在震源区深度上(15km)，绝大多数区域的流体渗透压均＞0.1MPa;在水库下方及偏西的区域，流体渗透压达

到了0.15MPa(图9)。然而，当水库下方不存在规模性断层时，水体渗透以垂直向下为主，流体渗透压在垂直方向上梯度最小。而两侧的断层对其渗透影响有限。水库蓄水29个月后，水库下方15km的深度上流体渗透压略＞0.1MPa(图10)，其值与基本模型A0－C－E在12月的状态相似。可见水库底部规模性断层的存在与否和断层带宽度对于水库蓄水所引起的流体孔渗透压分布有重要的影响。

图6 B1－C－E模型结果图Fig.6　 The simulation result of model B1－C－E.断层标号同图1

图7 B2－C－E模型结果图Fig.7　 The simulation result of model B2－C－E.断层标号同图1

图8 B3－C－E模型结果图Fig.8　 The simulation result of model B3－C－E.断层标号同图1

4　讨论

在水库加载及渗透过程中，引起断层带上应力变化的原因主要有2个: 1)库区荷载作用导致断层面上应力变化，对于以逆冲运动为主的汶川地震断层，库区荷载具有抑制滑动(增加稳定性)的功效(陈颙，2009) ; 2)水库蓄水形成的水头落差驱动水体向深部渗流造成的流体孔隙压变化所诱发的应力场扰动。流体孔隙压的增加会降低断层面上的有效正应力。此外，水的弱化作用会降低断层带的摩擦系数，从而使断层带的稳定性下降。对于汶川地震的发震断层，了解断层带上流体渗透压的变化是判断紫坪铺水库蓄水是否会实质性地降低断层带的剪破裂强度，最终导致灾难性结果的基础。

不同的模型参数对模拟结果的影响程度是本文重点考虑的问题。本文将渗透率随深度的变化纳入到模拟计算中，即模拟过程中不再将渗透率作为常数处理。断层岩和围岩的渗透率随深度的变化均依据实验结果。尽管渗透率随深度的变化规律是由实验确定的，但考虑到不同的岩石渗透率变化范围较大，本文还将实验得到的渗透率放大和缩小1个量级进行模拟计算。由于水的渗透不仅与渗透率相关，还依赖于渗透通道的空间尺度;虽然断层带的渗透率要远高于围岩，但断层带宽度相对于围岩尺度很小;因此只有同时考虑这2个参量才可能对断层带尺度如何影响流体渗透压的分布问题有深刻的了解。

图9 A0－C－F模型结果图Fig．9　 The simulation result of model A0－C－F．断层标号同图1

图10 A0－D－E模型结果图Fig.10 The simulation result of model A0－D－E.断层标号同图1

4.1渗透率不均匀性对模拟结果的影响

渗透率参数是否随深度变化对模拟结果有巨大的影响。当考虑渗透率随深度增加逐渐减小的真实情况时，流体渗透孔隙压的分布呈近似向下展布的水滴形状，即在垂直向下的20°～30°范围内压力梯度最小，随着角度的增加，其压力梯度逐渐增加(表明水体以向下渗透为主)。受水库下方断层产状的影响，其分布略向西偏。

常数渗透率模拟结果显示出截然不同的孔隙压分布特征。在近地表侧向梯度很小，与垂向梯度接近或可比，即形成近似半圆的孔隙压分布形状(图6)，显示出均匀渗透的特征。该结果与前人结果一致(Bell et al.，1978; Ge et al.，2009;周斌等，2010; Deng et al.，2010; Lei，2011)。其次，断层和围岩渗透率的比值(kf/kc)控制了流体孔隙压的分布状态。即kf/kc越小(均匀性越高)，侧向渗透能力越强。例如，当kf/kc= 102时，流体渗透压等值线近乎半圆形(图6) ;相反当kf/kc= 103和104时，流体渗透压等值线偏离半圆形分布越来越明显(图7，8)。但从图2可以看出，围岩和断层带的渗透率之比并非是常数，而是随着深度的增加逐渐增大。在近地表，二者相差不到一个数量级;而在20km的深度上，二者的差别可达4个数量级。由于在浅处，断层带和围岩的渗透性相差不大，流体从断层带向两侧围岩方向的渗透远比图8显示的要强，因此，图8显示的流体孔隙压分布形态严重地夸大了流体沿断层的渗透作用。

4.2断层带宽度对模拟结果的影响

对比图3－4和图10(渗透率值相同，但前者对应的断层带宽度为100m，后者为300m)可以看出，断层带宽度的增加导致水更容易沿断层带向下渗透。最直观的现象是流体孔隙压影响范围更大。此外，流体孔隙压的分布受到断层带产状的控制更加明显，即流体更多的是沿着断层带进行渗透。例如，当断层带宽度达到300m时，蓄水29个月后，汶川地震震源区的流体孔隙压可达到0.15MPa。而如果断层带宽度只有100m时，其流体孔隙压约为0.1MPa。

4.3紫坪铺水库蓄水诱发地震问题的探讨

水库蓄水所产生的附加质量(库体荷载)在断层面上引起的附加应力和水体通过缓慢渗透而引起孔隙压的变化共同影响着断层带的稳定性。为讨论水库蓄水可能产生的作用，可以考察水库蓄水而引起的断层带上的库伦应力变化(ΔS)。

式(3)、(4)中，Δσ和Δτ分别为垂直于断层面上的正应力和沿断层面上的剪应力变化，挤压状态下Δσ为正值，μ为静摩擦系数。负Δσ和正的Δτ将降低断层面的稳定性，即有促进断层失稳的作用;反之亦反。ΔP使有效正应力减小，其作用总是促使断层向不稳定的方向发展。前人对紫坪铺水库的库体荷载所引起的库伦附加应力已经做了计算和讨论，其库体荷载在地震断层面上引起0.1～0.2MPa的附加正应力和负的剪切应力(雷兴林等，2008;周斌等，2010)，并且随着距库底的距离增加，其附加应力逐渐减小。因此，如果仅仅考虑库体荷载，则蓄水增加了汶川地震断层的稳定性。模拟结果显示，当假设断层的宽度为300m时，映秀－北川断层在10～15km深度上的流体孔隙压力ΔP可达到0.15MPa，如若断层宽度较窄(100m)，其ΔP约为0.1MPa。因此，在蓄水29个月后，其ΔP值与Δσ大体相当，表明流体孔隙压基本抵消了由库体荷载所引起的附加正应力。从应力扰动的角度分析，在29个月里，紫坪铺水库蓄水不会导致震源区断层带稳定性的显著下降。但由于流体的下渗会随着时间的推移而逐渐积累，因此，维持更长时间的水头压力将导致断层深部的孔隙渗透压渐渐增加，并最终有可能超过库体荷载所引起的附加正应力。

多数模型(A－C－E、A－C－F和A－D－F)的模拟结果显示在库区下方5km深度范围内，其流体孔隙压均可达到0.4～0.5MPa。该孔隙压远远超过了库体荷载所引起的附加正应力，由此推断，紫坪铺水库蓄水会影响库区浅部断层的稳定性。在紫坪铺水库蓄水后，库区周围的地震活动明显增强，例如出现水磨震群、都江堰震群(周斌等，2010)等反映出断层稳定性对流体孔隙压的响应。

不同的学着对于流体渗透压的计算结果相差甚大。周斌等(2010)的计算显示，在约10km的深度上，与水库有水力联系的断层带其最大的流体孔隙压可达近1MPa，非断层带的ΔP也可达到0.3MPa。Gahalaut等(2010)则认为流体孔隙压的数值在0.4～0.5kPa。Ge等(2009)给出约0.05MPa的模拟结果。前人在模拟计算中所选取的水扩散系数(c值)也相差悬殊，雷兴林等(2008)选择c值为0.15～0.30m2/s; Ge等(2009)在计算中选用0.1m2/s的参数值;周斌等(2010)选择的围岩和断层岩的c值约为0.001～100m2/s(对应的渗透率为10－17～1.5×10－12m2) ; Gahalaut等(2010)所选的c值为1～10m2/s。此外，以上的模拟计算过程中均假设水扩散系数为不随深度变化的常数。因此，模拟计算结果相差甚远不足为奇。本文所用的水扩散系数是深度的函数，其值是依据渗透率实验结果推导得出的，所以其参数更为接近紫坪铺水库库区的实际情况。其模拟得到的流体孔隙压分布与水库下部的断裂带的产状及宽度存在明显的相关性。模拟结果表明，当考虑渗透率随深度变化时，其流体孔隙压的时－空分布与常数渗透率模拟结果截然不同。水库蓄水驱动水体向断层带渗透时，水不仅沿断层带向下渗透，同时也会沿断层带的走向渗透迁移。水沿断层带走向渗透的效应导致2维模拟得到的渗透压时空分布应该比实际的渗透压高。换言之，本文模拟给出的渗透压时空分布是上限值。

5　结论

依据汶川地震断层岩及围岩的渗透率实验数据，运用数值模拟的方法，对紫坪铺水库蓄水所导致的流体孔隙压时空分布进行了模拟。为对比断层结构参数、断层带和围岩的渗透率对

流体渗透压时空分布的可能影响，本文构建了8个渗透结构模型。模拟得到如下主要结果:

(1)依据断层带和围岩的渗透率随深度增加而减小规律所得到的流体孔隙压之时－空分布与前人假设渗透率为常数所得到的结果截然不同。显示出渗透率对流体渗透有极为重要的控制作用。

(2)渗透率相差1～2个数量级，基本不改变流体孔隙压的空间分布形态，但其孔隙压大小则有显著的差异，即渗透率越高，流体下渗越明显。

(3)断层带具有高渗的特点，因此水库下方的断层带宽度对流体渗透有显著的影响。断层带越宽，下渗作用越明显，流体孔隙压的分布与断层带的产状越接近。

(4)若断层带宽度为100～300m，2008汶川MS8.0地震发生时，其震源区(15km)的流体孔隙压上限为0.1～0.15MPa，与库体载荷在断层面上产生的附加正应力大体相当。该模拟结果可以作为进一步分析汶川地震的发生是否与紫坪铺水库蓄水存在关联的重要约束。

致谢感谢审稿人对本文提出的修改意见和建议。

参考文献

陈颙.2009.汶川地震是由水库蓄水引起的吗?［J］.中国科学(D辑)，39(3) : 257—259.

CHEN Yong.2009.Did the reservoir impoundment trigger the Wenchuan earthquake?［J］.Science in China(Ser D)，52(4) : 431—433.

胡先明，张永久，谢蓉华，等.2009.紫坪铺水库区小地震震源机制研究［J］.地震地质，31(4) : 676—685.doi: 10.3969/j.issn.0253－4967.2009.04.011.

HU Xian-ming，ZHANG Yong-jiu，XIE Rong-hua，et al.2009.Research on focal mechanism of small earthquakes in Zipingpu reservoir area［J］.Seismology and Geology，31(4) : 676—685(in Chinese).

雷兴林，马胜利，闻学泽，等.2008.地表水体对断层应力与地震时空分布影响的综合分析:以紫坪铺水库为例［J］.地震地质，30(4) : 1046—1064.

LEI Xing-lin，MA Sheng-li，WEN Xue-ze，et al.2008.Integrated analysis of stress and regional seismicity by surface loading: A case study of Zipingpu Reservoir［J］.Seismology and Geology，30(4) : 1046—1064(in Chinese).

刘保金，杨晓平，酆少英，等.2008.龙门山山前疑似汶川MS8.0地震地表破裂的浅层地震反射调查［J］.地震地质，30(4) : 906—916.

LIU Bao-jin，YANG Xiao-ping，FENG Shao-ying，et al.2008.Exploration of suspected surface rupture of the MS8.0 Wenchuan earthquake at frontal areas of Longmenshan using shallow seismic reflection［J］.Seismology and Geology，30(4) : 906—916(in Chinese).

刘远征，马瑾，马文涛，2014.探讨紫坪铺水库在汶川地震发生中的作用［J］.地学前缘，21: 150—160.

LIU Yuan-zheng，MA Jin，MA Wen-tao.2014.The role of the Zipingpu reservoir in the generation of the Wenchuan earthquake［J］.Earth Science Frontiers，21: 150—160(in Chinese).

马文涛，徐长朋，张新东，等.2011.紫坪铺水库与汶川地震关系的讨论［J］.地震地质，33(1) : 175—191.doi: 10.3969 /j.issn.0253－4967.2011.01.017.

MA Wen-tao，XU Chang-peng，ZHANG Xin-dong，et al.2011.Study on the relationship between the reservoir-induced seismicity at Zipingpu reservoir and the MS8.0 Wenchuan earthquake［J］.Seismology and Geology，33(1) : 175—191(in Chinese).

任俊杰，孙鑫哲，丁锐，等.2012.基于高程测量控制点和工程破坏对汶川地震都江堰段地表破裂特征的重新

界定［J］.大地测量与地球动力学，32: 37—40.

REN Jun-jie，SUN Xin-zhe，DING Rui，et al.2012.Determination of surface rupture characteristics of Wenchuan earthquake at Dujiangyan segment based on elevation controlling points and engineering damage［J］.Journal of Geodesy and Geodynamics，32: 37—40(in Chinese).

徐锡伟，闻学泽，叶建青，等.2008.汶川MS8.0地震地表破裂带及其发震构造［J］.地震地质，30(3) : 597—629.

XU Xi-wei，WEN Xue-ze，YE Jian-qing，et al.2008.The MS8.0 Wenchuan earthquake surface ruptures and its seismogenic structure［J］.Seismology and Geology，30(3) : 597—629(in Chinese).

杨清源，胡毓良，陈献程.1996.国内外水库诱发地震目录［J］.地震地质，18(4) : 453—461.

YANG Qing-yuan，HU Yu-liang，CHEN Xian-cheng，et al.1996.Catalogue of reservoir induced seismic events in the world［J］.Seismology and Geology，18(4) : 453—461(in Chinese).

周斌，薛世峰，邓志辉，等.2010.水库诱发地震时空演化与库水加卸载及渗透过程的关系:以紫坪铺水库为例［J］.地球物理学报，53: 2651—2670.

ZHOU Bin，XUE Shi-feng，DENG Zhi-hui，et al.2010.Relationship between the evolution of reservoir-induced seismicity in space-time and the process of reservoir water body load-unloading and water infiltration: A case study of Zipingpu reservoir［J］.Chinese J Geophys，53: 2651—2670(in Chinese).

Bell M L，Nur A.1978.Strength changes due to reservoir-induced pore pressure and stresses and application to Lake Oroville［J］.J Geophys Res，83: 4469—4483.

Biot M A.1941.General theory of three-dimensional consolidation［J］.J Appl Phys，12: 155—164.

Christensen N I.1996.Poisson's ratio and crustal seismology［J］.J Geophys Res，101: 3139—3156.

Deng K，Zhou S，Wang R，et al.2010.Evidence that the 2008 MW7.9 Wenchuan earthquake could not have been induced by the Zipingpu reservoir［J］.Bull Seism Soc Am，10: 2805—2814.

Gahalaut K，Gahalaut V K.2010.Effect of the Zipingpu reservoir impoundment on the occurrence of the 2008 Wenchuan earthquake and local seismicity［J］.Geophysical Journal International，183: 277—285.

Ge S M，Liu M，Lu N，et al.2009.Did the Zipingpu Reservoir trigger the 2008 Wenchuan earthquake?［J］.Geophys Res Lett，36: L20315.doi: 10.1029/2009GL040349.

Lei X.2011.Possible roles of the Zipingpu Reservoir in triggering the 2008 Wenchuan earthquake［J］.J Asian Earth Sci，40: 844—854.

Li S L，Lai X L，Yao Z X，et al.2009.Fault zone structures of northern and southern portions of the main central fault generated by the 2008 Wenchuan earthquake using fault zone trapped waves［J］.Earthquake Sci，22: 417—424.

Rice J R，Cleary M P.1976.Some basic stress diffusion solution for fluid-saturated elastic porous media with compressible constitutes［J］.Rev Geophys Space Phys，14: 227—291.

NUMERICAL SIMULATION OF TEMPORAL-SPATIAL DISTRIBUTION OF PORE—FLUID PRESSURE INDUCED BY ZIPINGPU RESERVOIR IMPOUNDMENT

YANG Yu1，2)YANG Xiao-song1)DUAN Qing-bao1)
1) State Key Laboratory of Earthquake Dynamics，Institute of Geology，China Earthquake Administration，Beijing 100029，China
2) HPT Laboratory，Department of Earth Sciences，Utrecht University，3584 CD Utrecht，The Netherlands

Abstract

The depth dependence of permeabilities of the fault rocks collected from the ruptures of 2008 Wenchuan earthquake and their host rocks were experimentally determined.On the basis of the measurement data，the temporal-spatial distributions of pore-fluid pressure induced by Zipingpu reservoir impoundment were numerically simulated.Modeled results indicate that the pore-fluid pressure is closely related to the patterns of the permeabilities.In other words，the results with the depth dependent permeabilities are strikingly different from those with the constant permeabilities.Existence of a fault hydro-dynamically connecting to the bottom of reservoir has a significant influence on the simulation of fluid infiltration owing to higher permeability in fault zones.The wider the fault zone is，the more obvious of fluid infiltrating along the fault.If fault width ranges from 100m to 300m，the pore-fluid pressure in the hypocenter region is about 0.1～0.15MPa，which is comparable with the normal stress caused by the reservoir impoundment on fault zone near the hypocenter of 2008 Wenchuan earthquake.This simulation provides an important constraint on the studies associated with the potential relationship between Zipingpu reservoir impoundment and Wenchuan earthquake.

Key wordsZipingpu Reservoir impoundment，Wenchuan earthquake，fluid infiltration，numerical simulation，temporal-spatial distribution of pore-fluid pressure

作者简介〔〕杨彧，女，1987年生，2014年于中国地质大学(北京)获矿床学专业硕士学位，目前为荷兰Utrecht大学岩石物理及岩石力学专业在读博士研究生，电话: +31－0681836909，E－mail: y.yang @uu.nl。

doi:10.3969/j.issn.0253－4967.2015.02.013

通讯作者:*杨晓松，男，研究员，E－mail: xsyang@ies.ac.cn。

收稿日期〔〕2014－01－28收稿，2014－06－22改回。

中图分类号:P315.9

文献标识码:A

文章编号:0253－4967(2015) 02－0510－14