0 引言

随着地震的发生，断层应力得到释放，由此产生的静态库仑应力变化将对后续破裂事件造成影响。区域静态库仑应力增加，将会促进附近断层后续地震的发生，而区域静态库仑应力减小，将减缓附近断层后续地震的发生，因此地震导致的静态库仑应力变化会影响其附近区域的地震活动性（刘方斌等，2018）。地震造成断层面上和邻区库仑应力的变化，被普遍用来研究其对周围断层的应力触发作用及其与后续地震之间的联系。前人研究表明，如果断层上的库仑应力变化超过0.1 bar，就可能触发后续地震并造成附近区域地震活动变化（Harris，1998）。

珊溪水库地处浙江省温州市，2000年5月开始蓄水，2002年7月28日库区发生M_L 3.5地震。此后，地震活动不断发生，至今已发生地震8 000多次，其中最大地震为2006年珊溪水库M_L 4.6地震。在2006年和2014年，库区先后发生震群型地震，最大震群为2006年珊溪水库M_L 4.6震群。据浙江省地震台网记录，2006年2月4日至11日珊溪水库发生M_L 4.0以上地震8次，截至2009年12月10日共发生地震1 063次。

前人对珊溪水库震群在地震精定位、发震构造、断裂运动特征等方面进行了深入研究，如：钟羽云等（2011）针对2006年珊溪水库M_L 4.6震群进行了精确定位并讨论了其地震序列的发震构造，此次地震震群发生在珊溪水库双溪—焦溪垟断裂东南端，震中较为集中，震源深度为5 km以上；杨福平等（2019）利用大地电磁测深法开展了珊溪水库地区三维地震深部构造探测，讨论了珊溪水库地震的深部构造环境；龚俊等（2021）利用GPS资料联合测震反演了2014年珊溪水库震群发震构造——双溪—焦溪垟断裂的运动学特征。然而，对于2006年M_L 4.6地震震群静态库仑应力变化及其对后续地震的影响分析较少。本文依据库仑破裂准则，计算此次M_L 4.6震群的静态库仑应力变化，并讨论其对后续地震活动性的影响，同时利用Dieterich（1994）提出的速率—状态摩擦定律，分析并预测研究区地震活动性变化。

1 区域地质构造

珊溪水库地处华南褶皱系东南褶皱带温州—临海坳陷。该坳陷区发育有NW向和NE向陡倾角逆断裂和逆走滑断裂。规模较大的NW向双溪一焦溪垟断裂为区域控制断裂，走向310°，倾向SW，由F_11-1、F_11-2、F_11-3三条近乎平行分支断裂组成，库区地震活动集中分布在该断裂附近（图 1）。其中：①F_11-1分支：长约10 km，倾角60°—70°，断裂带内挤压构造透镜体、断层泥及片理发育。该分支断裂形成于燕山晚期，破碎带宽约1 m，胶结良好，不利于库水下渗。②F_11-2分支：长度小于5 km，断层面近直立，局部倾向NE，倾角65°—70°，具正断性质，野外露头调查显示破碎带带宽10—20 cm，剖面出露地层显示为青灰色厚层状凝灰质砂岩。③F_11-3分支：倾角较陡，在水库附近展布长约8 km，在杜山村南水库边见断层出露，发育于上侏罗统灰绿色中厚层状凝灰岩中，破碎带节理较为发育，宽约5 m，断面上见近水平向擦痕，断层性质以走滑为主，兼逆断性质。该分支断裂在3个分支断层中破碎带规模最大、胶结程度相对较差、淹没段最长，破碎带内竖向节理发育，利于库水下渗（马志江等，2016）。

2 水库蓄水及对地震活动影响

珊溪水库位于浙江省飞云江上游河段，水库绝对高度156.8 m，坝长308 m，设计最大水位154.75 m。在水库大坝建设之前，库区周边25 km范围内无历史地震记载。水库于2000年5月12日开始蓄水，至2002年5月17日，水位达最高值133.8 m，同年7月28日，首次记录到库区地震，震级M_L 3.5，后地震活动逐年发生，最大地震为2006年2月9日M_L 4.6地震，此次地震最大影响烈度达Ⅵ度，对珊溪水库周边造成较大的经济损失和社会影响。据杨福平等（2021）的研究，研究区内侏罗系地层主要表现为中高阻特通性，易导致库水向深部下渗，三维大地电磁测深发现，F_11-3断裂附近出现“低阻通道”，导致地下水向深部渗透，断层面上的正应力被降低，对地震的发生有一定促进作用。

珊溪水库地震震级—时间及水位对应关系见图 2。由图 2可见，在水库蓄水初期，水位的快速变化与地震的发生有一定相关性，可能是受到水库蓄水的弹性载荷效应和库水渗流引起的空隙压力增加的影响，但随着蓄水时间的增长，库区岩层孔隙压趋向稳定，水库水位与地震发生的相关性降低，后续地震的发生可能与构造应力变化有关。2006年珊溪水库M_L 4.6震群发生前，水库蓄水时间已达4年，库区小震活动比较频繁，2005年9月11日库水位由140 m开始下降，至2006年此次震群发生之前，水位下降至124 m。水库蓄水过程造成的附加应力对区域断层的定量影响尚需进一步研究。

3 库仑应力变化计算

为了探究珊溪水库8个M_L≥4.0地震产生的静态库伦应力变化对后续地震的影响，计算其产生的静态库伦应力变化。

3.1 静态库仑应力

静态库仑应力触发是指，地震断层的同震位错引起地壳应力场永久变化后地震活动性变化。库仑破裂准则下静态库仑应力变化是，基于各向同性介质弹性半空间模型，根据发震断层同震位错量计算接收断层的库仑应力变化。静态库仑应力变化可用下式表达

$ \Delta \sigma_f=\Delta \tau_{\mathrm{s}}+\mu^{\prime} \Delta \sigma_{\mathrm{n}} $

(1)

式中，Δτ_s为断层滑动方向上的剪应力变化，当Δτ_s与接收断层的滑动方向一致时为正，反之为负；Δσ_n为断层面上的正应力变化，使断层两盘分离为正，挤压为负；μ^′为有效摩擦系数。若Δσ_f为正，则应力变化会促进目标断层的滑动，反之则抑制目标断层滑动。以往研究表明，当Δσ_f达到0.1 bar时，会触发后续地震破裂事件（Harris，1998）。

3.2 库仑破裂应力计算方法

依据库仑破裂准则，当岩石接近破裂时，库仑破裂应力为

$ F_{\mathrm{CFS}}=\tau+\mu(\sigma-p)-S $

(2)

式中，τ为主地震产生的应力张量在断层面上的剪切应力大小；μ为摩擦系数；p为孔隙压力，以沿断层滑动方向为正；σ为断层面上的正应力，以张应力为正；S为岩石内聚力。一般比较关注的是库仑破裂应力的变化，当μ和S不随时间变化时，库仑破裂应力的变化ΔF_CFS可由下式计算

$ \Delta F_{\mathrm{CFS}}=\Delta \tau+\mu(\Delta \sigma-\Delta p) $

(3)

为近似计算孔隙压力变化对库仑破裂应力变化的影响，假设介质是均匀且各向同性的，若静态应力发生变化，而流体自由流动前，空隙压力变化如下

$ \Delta p=-\beta^{\prime} \Delta \sigma_{k k} / 3 $

(4)

式中，β^′为Skempton’s系数，一般取值0.7—1.0；σ_kk为应力张量的对角元素之和。假设介质在断层带外是均匀且各向同性的，而在更具韧性的断层带内部介质是均匀且各向同性的，则Δσ_kk/3 = Δσ，代入式（3），可得

$ \Delta F_{\mathrm{CFS}}=\Delta \tau+\mu^{\prime} \Delta \sigma $

(5)

其中，μ^′ = μ(1 - β^′)，为有效摩擦系数，通常包括孔隙流体的影响以及断层带的材料特性。通过对比式（3）—（5），可知μ^′是σ_kk和Δσ的函数，则

$ \mu^{\prime}=\mu\left(1-\frac{\beta^{\prime}}{3} \frac{\Delta \sigma_{k k}}{\Delta \sigma}\right) $

(6)

有效摩擦系数μ^′通常取值为0—0.75，本文据Stein等（1992）的研究，取平均值0.4。数值实验表明，取不同μ^′值计算得到的库仑破裂应力的空间分布相差不大，但应力大小有一定差别（万永革等，2000）。

3.3 计算参数的选取

基于钟羽云等（2017）列出的2006年珊溪水库M_L≥4.0地震序列震源机制解（表 1），据库仑应力变化计算公式（3），计算8个M_L≥4.0地震的库仑应力变化，计算参数的选取参考浙江省地震局珊溪水库库区断层野外考察结果，如下：发震断层走向322°，倾角81°，滑动角3.75°，等效摩擦系数取值0.4。由于地震未造成地表破裂，计算库仑应力变化所需断层破裂长度、宽度及断层位错等信息，可由经验公式所得。

据Wells等（1994）提出的矩震级与破裂长度L和宽度W之间的关系，有

$ M_{\mathrm{W}}=4.49+1.49 \text{lg} L $

(7)

$ M_{\mathrm{W}}=4.37+1.95 \text{lg} W $

(8)

平均位错可由地震矩经验公式得到，公式如下

$ D=M_0 /(\mu A) $

(9)

式中，M₀为标准地震矩，μ为剪切模量，A为断面面积。地震震级M_S与M₀有以下关系

$ M_{\mathrm{W}}=\frac{2}{3} \lg M_0-10.73 $

(10)

$ \text{lg} M_0=1.5 M_{\mathrm{S}}+16.1 $

(11)

根据滑动角r，可计算出断层的水平位错h和垂直位错v，则有

$ h=D \cos r $

(12)

$ v=D \sin r $

(13)

3.4 结果分析

选取珊溪水库M_L 4.6震群及之后至2009年12月发生的1 063次地震，以珊溪水库2006年2月8个M_L≥4.0地震的发震断层作为源断层（地震断层参数见表 1），目标断层的几何参数选用构造应力驱使下的最优断层面解，设置泊松比ν = 0.25，杨氏模量Ε = 8×10⁴ MPa，利用Shinji Toda开发的Coulomb3.3程序（Stein et al，1997）计算8个M_L≥4.0地震的静态库仑应力变化，讨论此次震群的静态库仑应力变化及其对后续地震的影响。若大多数后续地震位于库仑应力增加区，则表明后续地震的发生受到震群触发作用，否则震群对后续地震产生延迟作用。计算结果见图 3。

图 3(a)给出珊溪水库M_L 4.6震群中8个M_L 4以上地震主发震断层库仑应力变化，其中紫色五角星为8个地震的震中分布位置，图 3(b)给出该震群库仑应力变化后的后续地震分布，图中黑色加号为后续地震。由图 3可见，库仑应力增加区呈NW和NE两个共轭方向，F_11-2断层两端和垂直方向的库仑应力增加，8个M_L 4以上地震的震中位于库仑应力影区，说明地震发生后震中附近区域库仑应力降低，并导致周边区域库仑应力增大。库仑应力变化结果表明，研究区库仑应力正值区域约占整个库仑应力变化场的45%，结合珊溪水库2006年2月9日至2009年12月的地震分布，在库仑应力增加区域所发生地震占比76.3%，且震中位置较为集中，说明后续地震在库仑应力影区受到抑制作用，分布相对稀疏。后续地震总体呈NW向条带展布，与双溪—焦溪垟断裂走向基本一致，其余地震随机分布在其他区域，根据库仑应力变化结果推测，此次珊溪水库震群受断裂展布影响较大。

由于震群发生在水库库区，库水对断层带的软化、弱化作用是触发地震的可能影响因素之一。为此，设置不同的断层带有效摩擦系数予以分析。图 4给出有效摩擦系数分别为0.0、0.2、0.4、0.6时珊溪水库库区库仑应力变化特征。

由图 4可见，在不同有效摩擦系数条件下，库仑应力增加区域和降低区域的总体分布相似，差异体现在沿发震断层区域的应力影响范围大小上，位于库仑应力增加区的余震次数占余震总数的75.6%—78.7%。有效摩擦系数的取值大小对库仑应力场与余震分布之间的关系影响不大，地震主要发生在区域应力积累区域，双溪—焦溪垟断裂构造活动对地震的发生起到控制作用，并且在水库蓄水加卸载过程中，可能由于库水沿破碎程度较高的F_11-3断层渗流而导致有效摩擦系数降低，同时孔隙压力的变化进一步对地震的发生起到促进作用。

为了增加研究的可信度，将珊溪水库库区库仑应力变化场应力增加区和下降区划分为AA′、BB′、CC′、DD′四个剖面，如图 4(c)所示，计算4个剖面的库仑应力变化，结果见图 5。由图 5可见，沿主发震断裂，即双溪—焦溪垟断裂方向的AA′剖面两端，库仑应力增加，后续地震大部分发生在库仑应力增加区域，说明库仑应力变化对后续地震的发生起到一定触发作用。沿近NS向BB′剖面，库仑应力降低，说明M_L 4.0以上地震的发生降低了这一区域的库仑应力，后续地震发生在该区域的可能性较低，表明库仑应力变化对该区域后续地震的发生起到一定延迟作用。

4 库仑应力变化对地震活动性的影响 4.1 背景地震活动

地震活动率是指一定区域内一定时期的地震活动特性。为了分析珊溪水库地区地震活动率是否受到2006年8个M_L 4以上地震产生的静态库仑应力变化的影响，引入Dieterich（1994）速率—状态摩擦定律，分析珊溪水库区域地震活动速率情况（Hainzl et al，2009）。在一定条件下，应力变化幅度可以用库仑应力变化代替。Dieterich模型可以用来预测一次地震发生之后，由于应力变化对周边断裂影响而造成的地震活动速率和地震发生概率的变化，其计算公式如下

$ R(t)=\frac{r}{\left[\exp \left(\frac{-\Delta \mathrm{CFS}}{A \sigma}\right)-1\right] \exp \left(\frac{-t}{t_{\mathrm{a}}}\right)+1} $

(14)

式中，R(t)为地震活动率；r为背景地震活动率；ΔCFS为库仑应力变化量；A为摩擦本构状态参数；σ为正应力；Aσ取值为0.3 bar；t_a为应力扰动的持续时间，也就是应力变化导致的地震活动率回到背景活动率的时间。

通过计算2006年珊溪水库M_L 4.6震群发生后的地震活动率和背景地震活动率的比值来表示地震活动率变化。其中背景地震发生率为一定圆域内完整性震级（M_c）以上的地震频次，使用震级—频度经验公式来描述研究区内地震活动性差异，根据G—R关系，有

$ \lg N(M)=a-b M $

(15)

式中，N(M)是以震级M（M＞M_c）为中心的小区间在一定时期内发生地震的次数，a表示研究区内地震活动程度，b表示地震震级数的比例。

采用2006年2月—2009年12月地震目录中1 063个地震得出研究区震级—频度分布，计算得到最小完整性震级M_c = 1.2，G—R关系式为lg N = 3.507 - 0.52M。

计算中选取0.01°×0.01°的网格单元，以0.6 km为扫描半径统计样本点数，得到研究区背景地震分布，结果见图 6，可见地震活动性范围呈椭圆形展布，远离圆心地震活动逐渐减弱。与图 2对比发现，与余震分布相似，背景地震集中分布在双溪—焦溪垟断裂附近（中部）。

4.2 参数计算

根据式（14），∆CFS是库仑应力变化结果，采用对数线性拟合方法来确定余震持续时间参数。也就是说，余震持续时间以震群发生后余震发生数目随时间变化的函数关系回到背景地震水平来确定。依据珊溪水库M_L 4.6震群发生后单位时间内的余震数目，对其进行线性拟合处理，得到余震活动水平降低到背景地震水平所需时间t_a = 601天，拟合关系见图 7。

4.3 地震活动趋势

使用Coulomb3.3程序计算双溪—焦溪垟断裂F_11-1、F_11-2及F_11-3中部的库仑应力变化，分别为-0.52 bar、-0.29 bar、1.24 bar，将库仑应力变化、余震持续时间、Aσ等参数代入式（13），得到R(t)/r随时间变化的函数关系。图 8给出F_11-1、F_11-2及F_11-3断裂上地震活动趋势，可见在此次震群发生后，F_11-1、F_11-2断裂库仑应力变化小于0，受到主震导致的库仑应力降低的影响，周边断裂上地震的发生被延缓，断裂附近地震活动性小于背景地震活动性，且F_11-2断裂所受抑制作用较大。图 8(b)显示，余震持续时间后基本无地震发生；在震群活动作用下，F_11-3断裂中部库仑应力增加1.24 bar，未来一段时间内其地震活动高于背景地震活动，601天后地震活动将主要分布在该断裂附近。由图 3(a)可见，库仑应力增加的NW区域与图 9(b)中该区在2007年后地震多发生在F_11-3断裂的实际情况一致。

自2002年7月珊溪水库第一次地震（M_L 3.5地震）发生，至2006年M_L 4.6震群发生前，库区小震震源深度基本分布在5 km以下，而自该震群发生至2008年以前，库区小震震源深度在0—10 km均有分布，2009年以来震源深度大于5 km，且集中分布在双溪—焦溪垟断层的F_11-3分支断层上（图 9），而华南地震区震源深度平均为10 km（龚俊等，2021），可见珊溪库区震源深度较浅，触发地震可能为水库诱发地震。作为双溪—焦溪垟断层库区淹没段最长的分支断层，F_11-3断层的地质构造特点是，胶结成度较差，破碎带竖向节理发育，有利于库水下渗，说明后续地震不仅受到主震的触发作用，还可能受到库水下渗导致孔隙压增大的诱发作用。然而，水库蓄水过程造成的附加应力对区域断层的定量影响仍需进一步研究。

5 结束语

通过对珊溪水库2006年2月4日—11日M_L 4.6震群库仑应力变化的计算，分析库仑应力变化对后续地震活动性的影响，结果表明：

（1）珊溪水库8次M_L 4.0以上地震的静态库仑应力变化及库水下渗导致的静摩擦力变化，对后续地震的发生起到一定促进作用，约76.3%的地震发生在库仑应力变化为正的区域。

（2）有效摩擦系数的取值大小对库仑应力场与余震分布之间的关系影响不大。

（3）计算F_11-1、F_11-2及F_11-3分支断层中部的库仑应力变化，发现F_11-1、F_11-2断裂上的库仑应力变化降低，F_11-3断裂上的库仑应力变化增加，且F_11-3断层库区淹没段最长，胶结成度较差，破碎带轴竖向节理发育，有利于库水下渗导致孔隙压增大，进一步表明未来地震将更多发生在F_11-3断裂附近。