四川大岗山水库蓄水对地震活动影响的数值模拟研究

朱家正， 孙玉军

中国地质科学院地球深部探测中心, 北京 100037

0 引言

川西地区位于我国南北地震带的南段,是构造形态复杂、地壳运动剧烈、地震活动频繁的地区.在深大断裂的控制下,块体之间的水平运动和垂直差异运动都表现得十分明显,而强烈地震大多数就发生在这些深大断裂带上(李碧雄等, 2014).同时,我国西南地区水力资源丰富,近些年在该地区修建了多个大型水电站,如向家坝、溪洛渡、白鹤滩、乌东德等大型水电站,大型水电工程建设所面临的水库地震问题是重要的科学问题.一方面我国西南地区是构造活动十分强烈的区域,构造地震活跃；另一方面,修建水库后,大型水电工程蓄水造成的附加载荷和渗流作用一般都会改变地区原有的应力状态.水库地震一般震级较小,但是其造成的附加应力改变可能与该区域初始的构造应力状态叠加,从而造成大的构造地震发生.如2008年发生在龙门山断裂带上的汶川地震,毫无疑问,如此大震级的地震发生是构造运动的直接结果(Chen, 2009),但是后续的研究表明建立在龙门山断裂带上的紫坪铺水库对断层的库仑应力影响为正值(Klose, 2012;周斌等, 2010;孙玉军等, 2012;雷兴林等, 2008;张贝和石耀霖, 2010),即对断层滑动有一定的促发作用,可能促使断层滑动提前发生(Ge et al., 2009).水库蓄水可以诱发地震是一个已经被证实的科学问题(Gough and Gough, 1970;Gupta, 2002;Gupta and Chadha, 1995;Talwani, 1997;夏其发, 1993;孙玉军等, 2012).我国最著名的水库地震是1962年发生在广东的新丰江水库地震,震级达到了6.1级(丁原章等, 1983;程惠红等, 2012).因此在构造活动强烈的区域修建大型水库等水电建设工程,其对该区域构造活动特别是大型断裂的影响,是未来值得深入考虑和探究的科学问题.该问题不仅涉及到大型水电工程的安全稳定性问题,同时对建设区域的防震减灾也有重要的参考价值.

四川大岗山水电站修建于2013年,属于大渡河上游修建的一系列水电站之一.由于龙门山断裂、鲜水河断裂、安宁河断裂等一系列大型断裂在此区域交汇,地质背景复杂,构造运动活跃,属于地震高烈度区.同时,大岗山水库属于高山峡谷型水库,坝高达210 m,对于这类特殊性水库,其蓄水造成的附加灾害对该区重要的大型断裂影响需要重点关注(程万正, 2013).

1 研究背景

1.1 构造地质背景

大岗山水库库区位于青藏高原东南缘川滇南北向构造带北端,川滇地块、华南地块、巴颜喀拉地块的边缘交界处(图1),是地壳运动剧烈、构造形态复杂、地震活动十分频繁的地区(张培震等, 2013；赵静等, 2018).龙门山断裂带、鲜水河断裂带以及安宁河断裂带在此区域组成Y型构造带.大凉山断裂、小金河断裂、玉农希断裂等主要断裂也在这一区域附近交汇,其中鲜水河断裂南段的磨西断裂以及大渡河断裂都直接通过水库蓄水区域(阮祥等, 2017).龙门山断裂带是一条长约500 km、宽约30～50 km、沿NE-SW方向展布的巨大断裂带,按照由西向东的顺序,龙门山断裂带主要包含龙门山后山断裂(茂县—汶川断裂)、中央断裂(映秀—北川断裂)和山前断裂(安县—灌县断裂,亦称彭县—灌县断裂、江油—灌县断裂).这些断裂(北西段)都以逆冲滑动为主,兼具一定的右旋走滑分量(Xu et al., 2008;邓起东等, 1994;陈国光等, 2007).从地震剖面上看,逆冲面呈铲状几何形态,上地壳向西陡倾,中地壳向西逐渐倾斜(Xu et al., 2008 ).张培震等的研究结果显示,龙门山断裂西南段存在(4±2.5)mm·a-1的缩短速率,以及(7.5±2)mm·a-1左旋走滑速率(张培震等, 2003).赵静等采用DEFNODE负位错反演程序估算了龙门山断裂西南端的闭锁程度与变形状态,并综合GPS反演结果和跨断层水准结果,认为目前龙门山断裂带西南端在大部分段落处于强闭锁状态,且依然有发生大地震的可能性(赵静等, 2018).鲜水河断裂带北起甘孜东谷附近,向南经过炉霍、道孚、康定一线,至石棉县安顺场一带逐渐减弱消失(郭长宝等, 2015).晚新生代以来,鲜水河断裂表现出强烈的左旋走滑运动,是松潘—甘孜造山带内部一条大型走滑断裂,横切了松潘—甘孜造山带的主体,系造山运动后期陆内变形的产物,晚新生代以来的位移总规模在60 km左右(许志琴等, 2007).以乾宁惠远寺一带为界,鲜水河断裂带北西段和南东段的活动速率有所差异,北西段活动速率为10～15 mm·a-1(闻学泽等, 1989),南东段小于10 mm·a-1,一般为 5 mm·a-1左右(郭长宝等, 2015).李铁明等通过重力、GPS资料计算结果,得出鲜水河断裂磨西段先进整体呈现左旋走滑,速率为4.41 mm·a-1(李铁明等, 2019).

图1 大岗山水库区域构造背景(a)及库区地形图(b)(a) 中橙色实线代表二级地块边界,蓝色实线代表一级地块边界,黑色实线代表主要断裂,黄色方框是研究区域所在位置,黄色五角星代表水库大坝位置; (b) 中蓝色填充区域为水库主要蓄水区域,黄色方块代表大坝所在位置.Fig.1 Regional tectonic background (a) and topographic map (b) of Dagangshan ReservoirThe solid orange lines in the figure (a) represent the boundaries of the secondary blocks. The solid blue lines represent the boundaries of the primary blocks. The solid black lines represent the main faults. The yellow box represents the location of the study area, and the yellow pentacle represents the location of the reservoir dam. The blue filled area in the figure (b) is the main water storage area of the reservoir, and the yellow square represents the location of the dam.

1.2 水库蓄水前后地震活动背景

大岗山水库在2013年建成后,大致分两阶段使得蓄水水位达到预定水位.于2014年3月开启第一阶段蓄水,水位达到1000 m；于2014年11月开启第二阶段蓄水,于2015年5月达到预定水位,水库水位由2013年蓄水前的960 m水位上升到1130 m的正常蓄水位.库区的地震数据表明(图2),水库蓄水后研究区域内地震发生频率随着水位的变化发生了明显的增加,特别是伴随着水位的两次大幅度变化,库区周边地震频度以及震级均发生增大.本研究将101.9°E—102.4°E, 29.2°N—29.8°N作为研究范围,以国家地震台网的地震数据为依据,绘出蓄水前后三年内库区地震平面分布图(图3).可以看到,在水库蓄水之前,地震主要分布在库区西侧鲜水河断裂,多为ML3级以下微震,在鲜水河断裂西侧呈条带状分布；2015年5月水库经过第二阶段蓄水达到正常水位后,该区域的地震分布仍然主要分布在库区西北的鲜水河断裂磨西段,但丛集性更强,且出现ML3～4级以上地震8次,ML4～5级地震3次,活动性变强.但在库区东部的龙门山断裂南段,在蓄水达到正常水位后几乎不再出现微震,库区南侧的大凉山断裂附近的地震活动性也有一定的降低.该现象是否与水库蓄水有直接关系值得深入研究.

图2 研究区域地震活动性与水库水位关系Fig.2 The relationship between seismic activity and reservoir water level in the study area

图3 水库蓄水前后研究区域地震活动情况(a) 2011年11月—2014年11月库区地震分布; (b) 2014年11月—2017年11月库区地震分布. 黄色五角星标识为水库所在位置.a, b, c, d分别代表龙门山断裂、鲜水河断裂、大凉山断裂以及锦屏山—小金河断裂.○代表震中位置,绿色为ML0～2级,黄色为ML2～3级,橙色代表ML3～4级,红色代表ML4级以上地震.红色椭圆所圈出区域为蓄水后地震活动性增强区域,蓝色椭圆所圈出的区域为蓄水后地震活动性降低区域.Fig.3 Seismic activity in the study area before and after the reservoir storage(a) The distribution of earthquakes in the reservoir area from November 2011 to November 2014; (b) The distribution of earthquakes in the reservoir area from November 2014 to November 2017. The yellow star indicates the location of the dam. a, b, c, d represent Longmenshan fault, Xianshuihe fault, Daliangshan fault and Jinpingshan-Xiaojinhe fault respectively. ○ represents the location of the epicenter. Green represents magnitude ML0～2, yellow represents magnitude ML2～3, orange represents magnitude ML3～4, and red represents above magnitude ML4. The red ellipse circles the area of increased seismicity after water storage, and the blue ellipse circles the area of decreased seismicity after water storage.

为了定量分析大岗山水库蓄水对该区域断层的影响,本文依据构造地质特征及精细的DEM数据建立了三维孔隙弹性有限元数值模型(图4),通过计算大岗山水库蓄水过程造成库区断层的库仑应力变化,定量分析了大岗山水库蓄水对地震活动性的影响.

图4 三维孔隙弹性有限元模型(模型中不同颜色代表不同材料参数)Fig.4 Three-dimensional poroelastic finite element model (Different colors in the model represent different material parameters)

2 研究方法

2.1 孔隙弹性体基本方程

孔隙弹性体的一般本构关系为

(1)

水的渗流扩散方程(暂不考虑饱和孔隙水不排水效应的影响)为(Bell and Nur,1978;Roeloffs, 1988)

(2)

(3)

其中c为水力学扩散系数,单位m2·s-1,k是介质的渗透率,单位m2,η是水的流体动力学黏滞性系数,单位Pa·s,泊松比、剪切模量、Skempton系数对结果的影响主要是通过改变扩散系数，本文将c作为计算参数.c的取值范围主要参考Talwani等(2007).

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

几何方程：

(4)

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

(5)

本文采用孙玉军等开发的三维孔隙弹性有限元模拟程序进行计算(孙玉军等, 2012).流体渗流和孔隙弹性耦合的研究方法不仅可以用于研究水库地震，同时在研究页岩气开采过程中压裂诱发地震时也被广泛应用(Hemam et.al., 2021; Wong et.al., 2021; Yang et.al., 2021).

2.2 建立模型

三维模型的建立考虑地层的垂向分层,模型表面层采用实际的地表地形数据,本研究采用SRTM30 m数字高程数据(https:∥www.usgs.gov/products/maps/topo-maps).设定模型的深度为25 km,模型范围为101.9°E—102.4°E,29.2°N—29.8°N,采用三棱柱单元进行网格剖分,并对库区部分进行网格加密,模型总结点数702259,总单元数1399412.数值模型可以考虑实际的分层性,但是由于缺乏实际渗流参数,本研究只区分了断层和非断层区域,对其设置了不同的模型参数.水平层状网格中,对蓄水区域进行加密,分辨率为60 m,非蓄水区域网格大小为900 m.垂向上网格分为15层,深度为25 km,底部四层网格无断层穿过,断层产状倾向北西西,倾角83°.

2.3 边界条件和计算参数

在研究区域范围内,对水库蓄水位以上的其他区域,给定力的边界和孔隙压力边界均为0,对水位以下地区给定力的边界和孔隙压力边界条件为

p(x,y,z,t)=ρωgh(x,y,z,t),

(6)

ρω为水的密度,g=9.8 m·s-2为重力加速度,h(x,y,z,t)为t时刻水库水位减去不同地点的地形高程值,即实际水深,p(x,y,z,t)为对应的水库蓄水造成的库底水压力,对于力的边界,模型侧边界和底边界均为滑动边界条件,对于孔隙压力,侧面给定为等渗流梯度,底面固定孔隙压力为0.

2.4 库仑应力计算

库仑应力的变化量定义为(Harris,1998;石耀霖和曹建玲, 2010)

Δcfs=Δτ+μ(Δσn+Δp),

(7)

Δτ为所考虑的断层面上剪应力的变化量,μ是摩擦系数,Δσn为断层面上的正应力改变量(拉伸为正),Δp为孔隙压力改变量.本文以coulomb-s代表库仑应力变化量.

同时,为了比较不同参数对计算结果的影响,分别建立不同的模型进行比较(表1).

表1 各个模型中计算参数的选取Table 1 Selection of calculation parameters in each model

鲜水河断裂南东段的磨西断裂主要走向为160°(李大虎等, 2015),在靠近库区的位置取其走向为166°,根据程佳、徐锡伟等对鲜水河南端水平和垂直滑动速率的总结,将断层滑动角参数设置为21°(徐锡伟等, 2003;程佳等, 2012).龙门山断裂南段大渡河断裂的断层参数参考前人的区域应力场及震源机制解研究来选取(易桂喜等, 2013;李君等, 2019;杜瑶等, 2016;闻学泽等, 2008).计算库仑应力所采用的接收断层参数选取如表2所示.

表2 计算库仑应力采用的各断层参数Table 2 The parameters of each fault used to calculate the change of coulomb stress (coulomb-s)

3 计算结果及分析

3.1 不同模型参数对结果影响

为讨论不同模型参数对计算结果影响,在水库蓄水后库区西北侧选取参考点1坐标为(102.12°,29.53°,-5 km),分别对比三组模型在该点处的孔隙压力以及库仑应力变化,库仑应力的计算采用鲜水河断裂磨西段的断层参数(表2).图5中coulomb-s代表库仑应力值.

图5a、5b为模型1计算得到结果.模型1中通过保持其他参数不变而只改变了断层的弹性模量以探究其对结果影响.从图中可以看出,在2013年12月水位第一次上涨时水库蓄水就对该参考点位置造成正的库仑应力影响,但该蓄水过程中孔隙压力几乎没有变化,说明此时水库蓄水的弹性载荷效应在起主要作用,随着水库在2014年12月水位快速增长110 m,三种不同模型参数对应的库仑应力变化量也增加明显,同时参考点的库仑应力随断层弹性模量的增大而减小,但由于三个模型并未改变渗流状态,因此三种模型所对应的孔隙压力基本一样(图5b),随着后期孔隙压力的增加,水库蓄水过程造成了正的库仑应力增加.

图5c、5d为模型2计算结果.模型2只改变了断层的扩散系数而保持其他系数不变.可以从结果中看出孔隙压力的结果在2015年4月出现明显改变；而库仑应力的三条结果曲线在2014年12月水位大幅度改变之前基本保持一致,计算结果显示此时水的渗流还未达到该深度.在2014年12月之后库仑应力随着孔隙压力的增大而逐渐变大,同时断层的扩散系数越大,在参考点处造成的库仑应力变化量则越小.

图5e、5f为保持断层的扩散系数不变,而改变非断层区域的扩散系数.结果显示非断层区域的扩散系数对结果影响较为明显,在第二阶段蓄水之后,当扩散系数c=1 m2·s-1时,参考点处的库仑应力和孔隙压力值都快速增长,而扩散系数低于该值的另外两个模型显示,库仑应力和孔隙压力变化较为平缓.

图5g、5h为改变非断层区域的排水泊松比所得到的计算结果.计算结果显示，在保持扩散系数不变的情况下，改变非断层区域的排水泊松比对孔隙压力的结果没有影响.水库蓄水初期库仑应力随着排水泊松比的增大而减小，随着时间的增加三种计算结果在一点交会后，两者关系转为正相关.该结果可以通过孔隙弹性的本构关系来理解，也即单纯的弹性介质下，库仑应力随排水泊松比增大而减小，但是根据孔隙弹性的本构关系，孔隙压力的增加，高排水泊松比更有利于库仑应力的增大.

图5 不同模型计算得到的库仑应力变化量以及孔隙压力随时间变化(a、b) 模型1结果; (c、d) 模型2结果; (e、f) 模型3结果. (g、h) 模型4结果; (b、h)中红色、绿色和粉色线重合.Fig.5 The change of Coulomb stress and pore pressure calculated by different models vary with time(a,b) The results of Model 1; (c,d) The results of Model 2; (e,f) The results of Model 3. The red, green and pink lines in (b) coincide.

以上几种模型中都可以看到在水库水位变化引起库仑应力变化主要受到两种效应影响,一种是蓄水造成弹性载荷变化引起的库仑应力变化,另一种是水的渗流造成孔隙压力变化而引起的库仑应力变化.但这两种效应在时间上的表现有所不同,蓄水造成弹性载荷变化所引起的库仑应力变化是瞬时的,如第一阶段蓄水之后库仑应力都随之快速增加,此时断层弹性模量对计算结果有较大影响,但此时孔隙压力并未明显改变.而第二种效应中孔隙压力改变所造成的库仑应力变化相对于水库蓄水过程有所滞后,如第二阶段蓄水之后的库仑应力增长都随着孔隙压力的变化而变化,说明水库达到设计水位后,主要是水的渗流作用影响库仑应力,而对此影响较大的参数是断层或非断层区域的扩散系数,扩散系数越大,孔隙压力越大,从而对参考点造成的库仑应力改变越大.

从孔隙弹性介质材料来看,水库蓄水过程所产生的这两种效应,可能会互相叠加增强,也有可能会互相抵消减弱.为了更详细的分析两种效应的作用,我们在计算库仑应力的时候可以分别对其进行考虑,即考虑水库渗流造成的库仑应力改变(弹性和渗流的综合效应)和不考虑水库渗流造成的库仑应力变化(弹性效应).

3.2 水库蓄水对不同区域地震活动性影响

综合上述三种模型,以及前人的研究成果(Talwani et al., 2007;陈建业等, 2011)选取模型5为最终的模型计算参数.在库区西侧的磨西断裂选取参考点1(102.12°E,29.53°N,-5 km),大渡河断裂选取参考点2(102.20°E,29.57°N,-5 km)(图3a),对水库蓄水后的地震活动性增强和减弱区域进行库仑应力和孔隙压力进行分析.对点1库仑应力计算选取磨西断裂的断层参数,对点2的库仑应力计算采用大渡河断裂的断层参数(表2).图中coulomb-e为不考虑孔隙压力影响计算得到的库仑应力变化量.

对比两图发现,图6a中在水库蓄水后不含孔隙压力的库仑应力coulomb-e(弹性效应)为正值,并且随水位增加而迅速增长,使断层更加危险；而图6b的结果与之相反,不含孔隙压力的库仑应力coulomb-e(弹性效应)为负值,随水位增长而快速下降,对断层起到抑制作用.这是由于水库蓄水区域跨越了龙门山和鲜水河两大断裂的南端,位置分别处于磨西断裂的下盘,以及大渡河断裂的上盘.因此库水对断层施加的弹性载荷使得磨西断裂的主要区域地震活动性增加,使得大渡河断裂的主要控制区域地震活动性减弱.但是,随着水位保持高位,水不断向断层中渗流,孔隙压力随之增大,使得库区两侧的库仑应力coulomb-s都转变为正值,断层趋向于更加危险.

图6 地震活动性增强区域与减弱区域选定点的库仑应力和孔隙压力变化(a) 地震活动性增强区域点1结果; (b) 地震活动性减弱区域点2结果.Fig.6 Changes in Coulomb stress and pore pressure at selected points in the areas of enhanced and weakened seismic activity(a) The result of point 1 in the seismic activity enhancement area; (b) The result of point 2 in the seismic activity weakened area.

3.3 对磨西断裂影响

将磨西断裂参数作为计算库仑应力的参数,得到2017年1月的库仑应力变化量与孔隙压力的平面分布(图7).孔隙压力增大的区域对应水库主蓄水区域,且随着深度的增加,孔隙压力值逐渐减小,如在深度5 km孔隙压力最大值为63 kPa,而在10 km和15 km孔隙压力最大值分别为22 kPa和10 kPa；同时由于断层的扩散系数一般较大,孔隙压力在断层区域一般较高.库仑应力变化量正值区域除了在水库库区范围外,主要分布在库区的北西和南东两个方向,同时随深度增加逐渐减小,在5 km深度处库仑应力变化量正极大值可达45 kPa,在10 km和15 km库仑应力变化量正极大值分别为17 kPa和7 kPa.库仑应力变化量的负值区域主要分布在库区东北侧和西南侧,在5 km深度负的极大值可达-6 kPa,在10 km和15 km负的极大值分别为-4 kPa和-2 kPa.深度为10 km以上时,水的渗流造成的孔隙压力较小,库仑应力变化量的区域性分布更为明显,库区西南方向出现了库仑应力变化量的负值区域,在水库蓄水后的地震增加区域(库区北西侧)看到明显的2 kPa左右的正值区域,说明水库蓄水对该区域地震的发生具有促进作用.

图7 2017年1月模型计算得到的深度为5 km (a, b),10 km (c, d),15 km (e, f)的库仑应力变化量以及孔隙压力分布(单位：Pa)○ 代表第二阶段蓄水之后距离剖面垂直距离范围2 km内的地震震中位置在平面上的投影,灰色代表震级ML2以下微震,黄色代表ML2～3级地震,橙色代表ML3～4级地震,红色代表ML4级以上地震.Fig.7 The coulomb stress and pore pressure distribution at the depth of 5 km (a, b),10 km (c, d), and 15 km (e, f) in January, 2017 (unit: Pa)○ Represents the projection of the earthquake epicenter position on the plane within 2 kilometers of the vertical distance from the profile after the second stage of water storage. Gray represents seismic magnitude below ML2, yellow represents magnitude ML2～3, orange represents magnitude ML3～4, and red represents magnitude above ML4.

以参考点(102.12°,29.53°)为中心位置对计算结果做纵向剖面(图8),孔隙压力主要分布在水库区域且有随着断层在地下的延展更快速地扩散(图8c),库仑应力变化量的分布与图7对应,在库区南北两侧形成负值区域,中心位置为正值；库区东侧以及南侧的5～10 km深处存在库仑应力变化量的负值区域(图8b、8d),西侧沿断层的延展分布库仑应力变化量正值区域.

图8 模型4计算得到的2017年1月不同剖面孔隙压力(a，c)和库仑应力变化量(b, d)分布图剖面A-A′为102.15°E的纵切面,剖面B-B′为29.53°N的横切面. 剖面位置已在图3a中显示.Fig.8 The pore pressure (a, c) and Coulomb stress (b, d) distribution of different profiles calculated by Model 4 in January 2017 (The positions of the two profiles are shown in Fig.3a)

从应力触发区域和地震活动的对应关系上,可以看到水库蓄水对10 km和15 km深度范围的地震在空间位置上具有明显的对应关系：在10 km深度处的应力影响区边缘出现了较为集中的ML2～3级地震分布,15 km深度范围内的ML3～5级的地震也集中分布在西北侧应力触发区域的边缘,说明水库蓄水对该区域地震发生具有促进作用,而不同深部的应力影区域的地震活动均受到抑制.

3.4 对大渡河断裂影响

以大渡河断裂参数(表2)作为计算依据得出2017年1月库仑应力变化量(图9a, 9c)及不考虑水渗流的库仑应力变化量(图9b, 9d)平面分布图.不考虑水渗流的库仑应力变化量在库区东侧存在一个明显的负值区域,在5 km深度靠近水库大坝位置的库仑应力变化量负极大值达到-12 kPa,考虑水渗流的情况下该值增加至-8 kPa左右,且与采用磨西断裂参数时一样，在蓄水区域出现正值区,5 km深度处库仑应力变化量正值极大值达到50 kPa,10 km深度处的库仑应力变化量分布与采用磨西断裂得到的结果相比,库仑应力变化量(考虑水渗流)负的极值增大1 kPa左右,达到-5.6 kPa.对比采用磨西断裂参数计算得到的结果可以看出,尽管采用的两条断层参数有所改变,但水库蓄水所造成的库仑应力变化量正、负区域(即应力触发和应力影区域)并未明显改变,除了造成库区的库仑应力变化量为正值外,在库区北西和南东为库仑应力变化量正值区(应力触发区域),库区北西侧相对较大；而在库区北东和南西则为库仑应力变化量负值区(应力影区域),库区北东负值更明显.该影响与水库蓄水后地震活动性比较一致.

图9 模型4计算得到的2017年1月深度为5 km(a, b)和10 km (c, d)的库仑应力变化量分布 (a, c)和(b, d)分别为考虑水渗流和不考虑水渗流得到了库仑应力变化量分布(单位：Pa). ○代表第二阶段蓄水之后距离剖面垂直距离范围2 km内的地震震中位置在平面上的投影,灰色为ML2级以下,黄色为ML2～3级,橙色代表ML3～4级,红色代表ML4级以上地震.蓝色三角形代表石棉地震震中投影.Fig.9 The distribution of the change of coulomb stress at the depth of 5 km (a, b) and 10 km (c, d) in January, 2017 for the model 4 (unit: Pa)(a, c) and (b, d) are the change of Coulomb stress with and without considering pore-pressure, respectively. ○ Represents the projection of the earthquake epicenter position on the plane within 2 kilometers of the vertical distance from the profile after the second stage of water storage. Gray represents microseisms below magnitude ML2, yellow represents magnitude ML2～3 earthquakes, orange represents magnitude ML3～4 earthquakes, and red represents earthquakes above magnitude ML4. The blue triangle represents the asbestos earthquake epicenter projection.

4 讨论与结论

本文采用三维孔隙弹性模型,研究大岗山水库蓄水前后水位变化对库区地震及断层的影响.根据地震台网数据，已知经过水位变化的第二阶段蓄水之后,水库附近的微震活动变化明显,首先,库区西北的鲜水河断裂南段附近地震震级增大,出现了ML3级以上的地震,且数量增多,丛集性更明显,由于蓄水位置处于磨西断裂下盘,促进了断层活动性,根据计算结果,在地震丛集增强的区域靠近库区的磨西断裂一侧为库仑应力变化量正值区(应力触发区)；而在地震活动性降低的库区北东区域,库仑应力变化量为负值(应力影区).但是随着渗流的持续,库仑应力变化量负值区域也会逐渐变为正值,地震活动性增强.根据前人研究成果,鲜水河断裂带的最大剪切应变速率为40～60×10-9/a(Wang and Shen, 2020).如果取最大剪切应变速率为50×10-9/a,假设弹性模量为80 GPa的情况下,应力累积速率为4 kPa/a,而蓄水造成的库仑应力变化量增加在浅层5 km深度最大可达到45 kPa(图10),这相当于11年的构造运动产生的应力累积.

地震增强区域和减弱区域的平面方向与临近的鲜水河断裂和龙门山断裂存在一定的夹角,这是否对应了两条断裂带地下的展布状态,需要进一步讨论.由于水库位于磨西断裂的下盘且主要蓄水区域都被其直接穿过,导致无论水的渗流还是库水的弹性载荷都促使其更加活跃；而大渡河断裂的上盘由于受到水库的垂向重力载荷,导致活动性减弱,地震减少.但是随着水的不断渗流,以及水位的稳定,水库弹性载荷效应造成的对地震的抑制作用会不断被削弱,库区的库仑应力可能都处于增加状态,即整个区域的地震活动性更强.根据地震数据,该区域在蓄水之后三年的时间中都没有再发生地震,说明现实中库水的弹性载荷对该区域地震的抑制效应比预计的更强；另外一种可能是水库区域地壳结构的黏弹性效应，黏弹性效应往往也造成库仑应力分布的不均匀，也可能对该区域造成负的库仑应力，从而造成渗流效应在该区域发生作用需要更长的时间.因此，在未来的研究中对于地壳结构的黏弹性效应也应该有所考虑.水库北东侧的断裂活动受到抑制，产生闭锁效应，造成龙门山断裂西南段的应变累积速率加快，水位的快速下降有可能造成此前积累的应变能快速释放，增加大震发生的可能性.

由于模型和实际断层详细数据的限制,没有给出鲜水河和龙门山断裂的铲状结构,实际情况中断层的倾角可能会随着深度增加而减小,因此水的渗流影响会在深部随着断层扩散至更远的区域,造成西北部区域的孔隙压力增加值比模型计算得到的值更大,该区域的断层可能更加活跃.另外，在鲜水河断裂的西北区域也可能存在与断裂平行的小断裂，如果这些小断裂与鲜水河主断裂参数一致，那么目前计算得到的库仑应力也适用于这些小断裂，从而促进断裂活动.因此未来如果获取了该区域更为详细的活动断裂探测数据，则可以利用该模型开展更为详细的计算分析.

本文的主要结论如下：

(1)从孔隙弹性材料的角度看,水库蓄水所造成的库仑应力改变主要有两种效应,一种是蓄水造成弹性载荷变化引起的库仑应力变化,另一种是水的渗流造成孔隙压力变化而引起的库仑应力变化.这两种效应在时间上的表现有所不同,蓄水造成弹性载荷变化所引起的库仑应力变化是瞬时的；而孔隙压力改变所造成的库仑应力变化相对于水库蓄水过程有所滞后.

(2)依据断层实际滑动参数计算得到的库仑应力结果来看,大岗山水库蓄水过程所造成的库仑应力改变使得鲜水河断裂磨西段库仑应力增加,而使得大渡河断裂库仑应力降低.这与水库蓄水后地震台网观测得到的磨西断裂西侧地震活动性明显增强,而大渡河断裂东侧地震活动性减弱的事实一致.

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