三峡水库蓄水后仙女山断裂北段数值模拟及其孕震机理研究

吴建超　陈蜀俊　陈俊华　蔡永建

摘要：基于线弹性理论有限元方法，建立三峡库首区的三维有限元模型，模拟计算了库首区尤其是仙女山断裂北段在蓄水后的全位移场和构造应力场变化。结果表明：三峡水库175m水位蓄水，造成仙女山断裂和九畹溪断裂北段处于全位移场归州低值区和新滩高值区的梯度带上：同时，仙女山断裂和九畹溪断裂北段深部3组主应力值较蓄水前陡增。这样的地段全位移差值变化明显，3组主应力值增高，剪应力和剪应变能易于集中，容易加速仙女山断裂和九畹溪断裂北端失稳而诱发水库地震。

关键词：三峡水库；仙女山断裂；数值模拟；孕震机理；水库地震

中图分类号：P315.7 文献标识码：A 文章编号：1000-0666（2016）02-0218-06

0 引言

三峡水库自2003年6月蓄水以来，在库区周边发生了数千次地震，绝大多数为M≤2.0非构造成因的诱发地震（车用太等，2009）。然而，伴随着三峡水库175m试验性蓄水进程，2008年11月22日、2014年3月27日和3月30日，仙女山断裂北段附近分别发生了M_S4.1、M_S4.2和M_S4.5地震，引发了社会各界对于三峡水库是否诱发了地震的讨论。

笔者参加了2014年3月30日秭归M_S4.5地震的应急调查工作，该次地震震中烈度为Ⅴ度，有感范围较大。在郭家坝镇头道河村、郭家坝村等地出现Ⅵ度烈度异常区，且造成了一定程度的房屋破坏和财产损失。根据地震监测记录波形分析，认为此次地震为构造型水库地震，与仙女山微地块的活动有关。

值得注意的是，蓄水后三峡库区的地震活动多发生在仙女山断裂和九畹溪断裂的围限区附近，属于这两条断裂的北端部。断层端部构造应力场及孕震机理一直是国内外研究的热点。断层的末端扩展受阻，往往会在断层的端部产生应力集中而引发地震（Obara et al，2009；孔礼健等，2009）。而三峡水库蓄水后的形变观测结果显示，仙女山断裂北段地块的倾斜场发生了改变，由蓄水前的正断层运动转为蓄水后的逆断运动性质（张燕等，2005；邢灿飞等，2006）。因此，综合考虑蓄水效应，研究蓄水后仙女山断裂北段构造应力场及其孕震机理具有重要意义。

本文应用三维有限元方法模拟计算了三峡水库蓄水对仙女山断裂北段构造应力场的影响，并对其孕震机理进行了分析和讨论，以期对三峡水库诱发地震机理作进一步探索。

1 地震活动概况

三峡水库自2003年6月蓄水以来，仙女山断裂北段附近共发生了41次M_S≥2.0地震，其中M_S4.0～4.9地震3次，M_S3.0～3.9地震3次，M_S2.0～2.9地震35次，且大部分位于长江南岸，呈密集面状分布特征，如图1所示，地震目录取自湖北省地震局台网中心。

车用太等（2009）认为，2008年11月22日胡家坪M_S4.1地震是在水库水体荷载与库水下渗的共同作用下沿仙女山断裂发生的构造型水库诱发地震；吴建超等（2012）认为胡家坪M_S4.1地震发生在等效应力高值区与低值区的梯度带上，为等效应力变化最为剧烈、最易产生应变能积累和应力集中的区域；陈俊华等（2013）从震群活动时序特征、震源体大小、震源机制及地震波形的频谱特征等方面进行研究，认为2012年10月31日屈原镇M_S3.2地震为水库触发型构造地震；徐长朋（2012）对三峡库区119次小震的震源机制解进行研究，认为仙女山断裂北段地震震源机制以走滑为主。

上述研究多针对蓄水所诱发的地震，未对蓄水后仙女山断裂上的位移场和深部应力场变化进行综合研究。因此，有必要通过对蓄水后仙女山断裂北段的孕震动力学环境进行数值模拟，从上地壳深部探究其孕震机理。

2 研究区三维有限元模型

本文构建的三维有限元模型范围为110km×85km，深度取至上地壳底部埋深12km处。魏贵春等（2013）利用双差定位法对三峡数字遥测地震台网记录到的秭归附近地区73个M_S≥1.0的地震进行了重新定位，结果显示，绝大多数地震事件的震源深度均在7km左右。这一深度刚好位于本文三维有限元模型的中部，计算过程中避免了受边界条件的影响，结果相对准确。因此，本文构建的三维有限元模型是合理与可靠的。图2为考虑地形地貌、断裂和构造层建立的三维模型，模型以（110°00′00″N，31°56′35″E）为坐标原点，X轴正方向为SE33°，Y轴正方向为NE33°，z轴垂直于水平面，向上为正。

模型按上地壳岩石层的物性参数可划分为沉积盖层和基底构造层。其中，沉积盖层厚约7km，主要为花岗岩、灰岩和砂岩；基底构造层层厚约5km，主要为花岗岩。依据陈蜀俊（2005）研究成果，断裂带的参数分别取为：弹性模量E=1.07×10¹⁰Pa，泊松比v=0.24，密度ρ=2570kg·m^-3。有限元网格划分以满足计算精度为要求，对仙女山断裂北段附近的网格再做细化。本次计算模型单元数共计18862个、节点数共计22498个。

参照陈蜀俊等（2005，2006）的研究成果和三峡地区地震加密观测资料（马文涛等，2010），本文选用的加载方案为：主压应力方向为NE-SW向，大小为790MPa，主张应力方向为NW-SE向，大小为540MPa。有限元模型采用三维线弹性模型，除作了约束的部位外，每个节点有3个自由度。

3 计算结果分析与讨论

三峡水库蓄水提升了长江上游水位，并淹没了长江部分支流河段。对库水影响范围内的水文地质条件与岩石物性参数进行研究，是蓄水背景下构造应力场数值模拟的基础。根据车用太等（1998）、刘远征等（2010）对水库诱发地震成因机理的研究成果，本文着重考虑了三峡水库蓄水后，水体荷载作用、库水下渗作用和孔隙水压作用的影响，并计算得到了蓄水后三峡库首区的全位移场和构造应力场。

3.1 全位移场

为更具体地分析三峡水库蓄水对仙女山断裂北段全位移场的影响，将175m水位与蓄水前全位移值相减，得到了三峡水库175m水位相对于蓄水前的全位移差值图（图3），与杜瑞林等（2004）给出的蓄水后GPS形变场观测结果较为一致。

蓄水造成的垂直形变远大于水平形变。因此，垂直位移对全位移作出的贡献最大。由图3可知，175m水位相对于蓄水前的全位移差值分布不均匀，出现了归州和官渡口等几个低值区，最大值出现在巴东附近，而最小值出现在归州一带。极值区大部分位于长江南岸，紧邻极值点的等值线梯度大，其后逐渐变缓，极值区等值线长轴方向与断裂带走向基本一致。值得注意的是，仙女山断裂和九畹溪断裂北段刚好处于归州低值区和新滩高值区的梯度带上。这样的地段全位移变化明显，易于剪应力和垂直剪应变能的集中，断裂容易失稳错动而引发地震。

将模拟结果与观测结果进行对比分析可知，数值模拟获得全位移场与GPS观测形变场所反映的研究区垂直形变特征及趋势基本吻合，全位移场的极低值区在归州-新滩-带出现（杜瑞林等，2004）。虽然上述模拟结果与GPS观测结果所反映的运动趋势整体上基本一致，但在局部区域还是存在差异，例如模型的边界等。这可能是由于有限元计算中对于模型边界的加载和围限作用所致，而其对处于模型中部的仙女山断裂附近则影响较小。因此，本文关心的仙女山断裂北段的数值模拟结果是可信的。

3.2 构造应力场

蓄水前后构造应力场的变化，直接反映了由蓄水引起的构造应力的变化，是蓄水引起库首区孕震环境变化的有力证据之一。为了详细了解三峡水库在蓄水进程中产生的主应力场变化，尤其是对仙女山断裂深部应力变化进行详细的剖析和研究，本文对模拟结果做了剖面输出。以坐标点（39400，0，-8000）作为切面点，以向量{-1，1，0}作为切面的法向向量，切出了跨仙女山断裂和九畹溪断裂的剖面线（图1），蓄水前后跨仙女山断裂和九畹溪断裂的最大主应力结果如图4所示。

由图4可以看出，伴随着蓄水进程，跨仙女山断裂北段最大主应力较蓄水前出现明显变化。蓄水前，深部最大主应力在郭家坝至新滩一带没有明显变化，呈舒缓波状分布，应力值自上而下缓慢递增。而蓄水至135m水位后，最大主应力等值线在郭家坝和新滩间，即仙女山断裂和九畹溪断裂深部出现了较大变化，具体表现为两组断裂带深部的最大主应力陡增，呈“V”字形，断裂带内部最大主应力均大于断裂带外部的应力大小。在175m水位下模拟计算得到的最大主应力云图基本同135m水位下的类似，只是在模型顶部有应力值的变化。

本文的研究重点是仙女山断裂北段深部构造应力场的变化，考虑到该断裂北段地震的震源深度多为7km左右。因此，取Z=-7km作为输出层，得到了蓄水至175m水位后，库首区7km深度处跨仙女山断裂的主应力场的结果。从图5可以看出，175m水位蓄水后，在仙女山断裂和九畹溪断裂带上均出现了3组主应力大小陡增的现象，且均大于断裂带两侧的主应力值。而两条断裂带中间所挟持的区域则明显地出现了主应力减小，总体呈现“M”型。

上述数值模拟结果集中表明，三峡水库蓄水给包括仙女山断裂和九畹溪断裂的整个库首区的应力状态和位移场带来了明显变化。

4 应力场变化对孕震环境的影响

大地构造上，仙女山断裂北段位于黄陵背斜的西南缘。黄陵背斜为一古老的刚性地质体，在地质历史中长期表现为隆起，对局部应力作用产生强大的反作用力。因此，黄陵背斜周缘的隆起区存在垂直向上的作用力，从而导致上部地壳引张（高士钧，1992）。

在三峡水库蓄水进程中，水体荷载会在库底岩体内产生弹性附加应力的响应，进而将影响到周边断层稳定性。在不同的构造应力环境中，附加应力对断层稳定性所起的作用是截然不同的。在拉张应力环境中，最大主应力方向近于垂直，最小主应力方向近于水平，水体荷载与初始构造应力场叠加，将会使莫尔圆半径变大，增加逼近外包络线的机会，进而加剧断层破裂而发生地震（周斌等，2010）。

该理论用于复杂的三维断面并不简单，而用相近的办法和思路，也可获得相似的认识。即蓄水导致断裂带上构造应力场发生变化，并叠加在原有的构造应力场上，形成了对构造应力场及孕震环境的扰动和影响，将产生软化、促滑、诱震作用，使原本相对稳定的断裂发生活动或原本就存在的地震危险趋势增加，甚至使未来地震震源机制和时、空、强度特征发生改变。

仙女山断裂北段主应力的模拟计算结果表明，蓄水后断裂带上的应力加载效应显著。一般而言，加载效应有利于局部应力积累，可能加速地震发生或提高局部地震危险度。卸载效应有利于区域应力释放，可能延迟地震发生或降低区域地震能量。

张秋文等（2005）研究认为，几个兆帕的应力场变化，虽不足以直接导致地震，却可以使地震发生的时间提前、发震概率提高。王绍晋等（2005）研究发现，澜沧江漫湾水库蓄水前3.0≤M_S≤4.0的历史地震震中均在库中心10km以外，震源机制解平均应力场主压应力方位为SSE或近SN向，蓄水后震中向库心汇集，震源机制解表明断层面错动具有较大的正倾滑动分量，主应力方位则变为近EW向，方位向东偏转了50°以上。

综上所述，三峡水库175m水位蓄水后，仙女山断裂北段浅部构造应力场的动态变化、调整蕴含着足以引起重视的地震动力学意义。仙女山断裂和九畹溪断裂的北段之间全位移变化明显，3组主应力大小陡增，应力加载效应明显，容易加速仙女山断裂北端部失稳而引发地震。

5 结论

本文模拟计算结果虽然在局部区域和模型边界上与观测结果存在一定的差异，但总体变化趋势基本一致。仙女山断裂北段深部构造应力场的变化，能够综合反映蓄水效应，而构造应力场局部应力状态的调整又与地震的孕育存在直接联系。因此，本文通过三维数值模拟研究蓄水后库首区构造应力场的变化，将有助于探讨蓄水后断裂带深部的孕震机理及其地震动力学意义。

通过对计算结果的分析讨论，我们获得了以下主要认识：

（1）利用三维有限元方法，计算得到了三峡库区175m水位相对于蓄水前的全位移差值结果，与蓄水后GPS形变场观测结果较为一致。沿归州、新滩和周坪三地出现了全位移场低值区，反映了175m水位蓄水造成的上地壳表层全位移场变化的基本状况。

（2）伴随三峡水库135m和175m水位的蓄水进程，仙女山断裂和九畹溪断裂浅部最大主应力较蓄水前出现明显变化，表现为两组断裂带深部的最大主应力值陡增，且均大于断裂带两侧的数值。

（3）三峡水库蓄水导致了仙女山断裂北段浅部构造应力场和全位移场的动态变化、调整。蓄水后，地震多发生于仙女山断裂和九畹溪断裂的北端附近，应力易于在这两组断裂端部积累。

（4）仙女山断裂北段在三峡地区的位置重要，175m水位蓄水以来，库首区M≥2.0地震多发生于这个区域。综合地层岩性、断裂构造、渗透条件等，运用水-岩耦合数值模拟手段，探究该区域地震的成因机理是我们下一步研究的目标。