华北平原地下水开采对区域地震活动性的影响

窦甜甜, 程惠红, 周元泽, 石耀霖

中国科学院大学地球与行星科学学院，计算地球动力学重点实验室， 北京 100049

0 引言

研究表明地下水开采、水库蓄水、页岩气开发、深井注入等对地震发生存在显著影响，人类活动诱/触发地震已成为地球物理学研究的一个重要领域(Gupta, 2002; González et al., 2012; Keranen et al., 2014; Amos et al., 2014; Segall and Lu, 2015; Kundu et al., 2015; Bao and Eaton, 2016; Grigoli et al., 2018; Kim et al., 2018; Wetzler et al., 2019;Yang et al., 2021; Zhou et al.,2021).对于区域性地下水开采，大规模的地壳卸载不仅能直接扰动应力和应变积累过程(González et al., 2012; Amos et al., 2014; Gambolati and Teatini, 2015; Kundu et al., 2015; Wetzler et al., 2019; Carlson et al., 2020; Pang et al., 2020)，也会引起区域孔隙压力变化并以扩散的形式向周围传播(Rice and Cleary, 1976; Talwani and Acree, 1984; Kundu et al., 2015; Segall and Lu, 2015; Bhattacharya and Viesca, 2019).Amos等(2014)对地下水开采引起的弹性应力变化进行了研究和分析，并认为卸载降低了断层面上的正应力.González等(2012)根据弹性位错模型计算了Lorca地震的库仑应力变化，认为该地震可能是由地下水卸载引起的.Kundu等(2015)通过模拟计算得出2015年Nepal地震以及主喜马拉雅逆冲断裂区发生的地震均受到印度河—恒河平原地下水超采的影响.Segall和Lu(2015)分别计算了地下水开采引起的弹性应力和孔隙压力变化对地震活动的影响，并强调了孔弹耦合应力的作用.庞亚瑾等(2016)模拟计算并分析了华北平原地下水开采对地壳应力的影响，但未考虑孔隙压力变化.Wetzler等(2019)据观测数据发现地下水快速抽取后Kinneret湖区域发生两次异常的浅地震群，通过孔隙弹性模型认为与地下水开采有关.根据孔隙弹性耦合理论，孔隙压力变化会引起应力场改变，同时孔隙压力也受到平均正应力的影响，两者相互作用共同控制断层的稳定性(Rice and Cleary, 1976; 程惠红等, 2012; Segall and Lu, 2015; Wetzler et al., 2019; Hemami et al., 2021; Barbour and Beeler, 2021).

目前，地下水仍是人们生活和灌溉的主要来源，特别是我国人口相对密集的华北平原，自20世纪60年代以来该地区地下水存在明显的长期亏损，山前平原区(包括石家庄、邢台等地)浅层地下水水位下降达30～50 m(费宇红等, 2009).2000年后，华北地区全年水供应约74%来源于地下水(Zheng et al., 2010).苏晓莉等(2012)利用GRACE重力观测数据并与相关水文模型比较得出华北地下水水位下降速度约0.5 cm·a-1.而随着“南水北调”工程的实施，近年来该地区地下水水位下降速度有所变缓(周志博等, 2020).

华北平原作为中国大陆典型强震多发区，近几十年来持续的地下水开采造成的地壳卸载和孔隙压力减小可能会对其产生较大影响.因此，本文基于孔隙弹性耦合理论，定量计算1959—2016年间华北平原地下水开采引起的区域应力和孔隙压力变化，并根据历史地震参数计算库仑应力变化，分析并探讨地下水开采对区域地震活动性的影响.

1 计算方法及模型建立

1.1 华北地区三维模型

华北平原位于中国大陆东部，西邻山西断陷盆地带，东邻郯庐断裂带(张培震等, 2013).自晚中生代以来该地区经历太平洋深俯冲和弧后扩展作用，活动断裂发育，地震活动强烈，历史上曾多次发生8.0级以上的破坏性地震(Liu et al., 2011).近代华北地区也多次记录到7.0级以上的强震，如1966年邢台7.2级地震、1976年唐山7.8级地震(Liu et al., 2011; 王辉等, 2011) (图1).研究表明华北地区地震活动存在明显的“期幕”特征，时间上呈现活动-平静期交替出现的周期性，空间上表现为地震带之间的丛集迁移特征(Liu et al., 2011; 尹晓菲等, 2020).

图1 华北地区强震分布及1959—2016年地下水水位变化(Pang et al., 2020)XFF:夏垫—凤河营断裂; TNF:唐山—宁河断裂; HJF:河间断裂; XHF:新河断裂; TXF:汤西断裂. Fig.1 Distribution of strong earthquakes and groundwater water level change from 1959 to 2016 in NCP (Pang et al., 2020)XFF: Xiadian-Fengheying fault; TNF: Tangshan-Ninghe fault; HJF: Hejian fault; XHF: Xinhe Fault; TXF: Tangxi Fault.

中、新生代时期，西太平洋板块俯冲导致华北克拉通破裂(朱日祥等, 2012; Zheng et al., 2017; Ma et al., 2022)，发育了一系列活动断裂，以北东向的高倾角右旋走滑张扭性断层为主(邓起东等, 2003; 张培震等, 2013).区域内主要活动断裂有夏垫—凤河营断裂(XFF)(何付兵等, 2013)、唐山—宁河断裂(TNF)(李志义和虢顺民, 1979; 张素欣等, 2020)、河间断裂(HJF)(杨家亮等, 2010)、新河断裂(XHF，又称百尺口断裂)(徐杰和方仲景, 1988)和汤西断裂(TXF)(韩慕康和赵景珍, 1980).远震波形数据及流动台站地壳观测数据显示华北地区地壳厚度平均约40 km(危自根等, 2015; 王椿镛等, 2017)，其中上地壳10～12 km，中地壳8～10 km，下地壳10～15 km(王椿镛等, 2017).

地貌上，华北平原属典型平原地貌，其上覆盖第四纪含水层(Foster et al., 2004; 孟素花等, 2013).地貌构成以山前冲、洪积倾斜平原，中部冲、湖积多层叠积平原和东部冲、湖积为主夹数层海积层的滨海平原为主(孟素花等, 2013).根据水文地质特征，可分为山前冲洪积平原、中部冲积湖积平原、古黄河冲洪积平原和东部冲积海积平原，其中山前平原颗粒粗，具有较强的渗透能力，中东部平原岩性逐渐变细，多为粉土、粉质黏土和黏土.

基于以上地质背景资料建立了华北地区三维有限元模型(图2).模型采用分层分块结构，深度上包括沉积层及上、中、下地壳在内，共40 km，断层切割深度15 km，地表地形采用SRTM 90 m数字高程数据(http:∥dds.cr.usgs.gov/srtm/).根据水文地质特征，开采区沉积盖层可划分为4个单元.采用三棱柱单元进行网格剖分，整个模型节点数为378690，单元数为698572，如图2a.

图2 华北地区有限元模型(a) 三维地质模型; (b) 边界条件.Fig.2 The finite element model of NCP(a) Three-dimensional geological model; (b) Boundary conditions.

1.2 计算方法

地下水开采不仅对地壳产生卸载作用，同时引起孔隙压力变化，孔隙弹性耦合理论可以较好地分析流体与固体骨架之间的耦合作用(Talwani and Acree, 1984).根据Biot(1941)固结理论，Rice和Cleary(1976)推导得出了孔隙弹性耦合方程，经众多学者研究和发展，目前已经广泛应用于分析岩体在流体注入/抽取及水体载荷作用下岩石应力场和孔隙压力场的变化(雷兴林等, 2008; 程惠红等, 2012; 孙玉军等, 2012; Segall and Lu, 2015)，其本构关系为

(1)

(2)

在已知断层参数的情况下，库仑应力变化可以由下式计算得到：

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

(3)

其中，Δσn是断层面上正应力变化，Δτ是断层面上沿滑动方向的剪应力变化.由于大多数情况下缺少研究区域初始构造应力场，库仑应力变化成为定量分析地震危险性的常用方法，被广泛应用于地震触发问题研究(King et al., 1994; Stein et al., 1997;石耀霖和曹建玲, 2010).由孔隙弹性模型计算得到应力张量和孔隙压力变化量，可以计算特定震源参数下的库仑应力变化.若结果为正，表明断层滑动危险性增加；反之，则表明断层趋于稳定，从而了解分析地震相对更容易发生于哪些地方，何种断层.

1.3 计算参数和边界条件

根据华北地区岩石圈P波和S波速度结构(Pang et al., 2020; 石富强等, 2020)可以给出模型弹性参数.相关研究表明(Talwani et al., 2007)孔隙弹性介质的扩散系数范围一般在0.1～10 m2·s-1，本文基于华北平原沉积盖层的水文地质特征给定计算所用的扩散系数.模型具体计算参数见表1，计算总时间为1959—2016年(58年)，计算时间步长为1年.

表1 模型计算参数Table 1 Model calculation parameters

由于华北平原地下水开采主要集中于上部含水层，其深度相对模型深度较浅，且在开采时间尺度内可忽略地幔松弛效应而仅考虑弹性作用.因此，地下水的卸载可近似为对模型上边界施加垂直向上的法向力.模型弹性位移场边界条件设定为：底边界、侧边界法向固定，切向自由，上边界自由，其中开采区内施加垂直向上的拉力，大小为地下水位变化、孔隙度和水的容重之积；孔隙压力场设定为：底边界、侧边界渗流梯度为0，上边界由水体卸载量给定，如图2b.相关研究表明华北地区含水层孔隙度范围在0.03～0.21(Pang et al., 2020)，本文计算中设模型具有均匀孔隙度，大小为0.18.

2 计算结果

2.1 华北平原地下水开采引起的位移场变化

地下水开采对地壳产生卸载作用，其引起的位移变化主要受弹性卸载控制，空间上与地下水水位下降分布相对应.图3给出了华北平原地下水开采引起的不同深度上的垂直位移变化.由于水体卸载，华北地区不同深度上都呈现出一定抬升，最大抬升出现在水位下降最明显的邢台、石家庄附近.开采区内地表平均抬升量大于5 mm，最大抬升35 mm，位于邢台以北(图3a).位移变化结果随深度增加而减小，上地壳上表面(地下5 km)和中地壳上表面(地下15 km)最大垂向位移分别为28 mm(图3b)和18 mm(图3c).相较于垂直方向上的位移变化，水平方向上位移变化较小，结果在5 mm以下.开采区以外的位移场变化不大，基本可忽略.

图3 华北平原地下水开采引起的不同深度上的垂直位移变化(a) 地表,蓝色三角为GPS站点; (b)地下5 km; (c)地下15 km.Fig.3 Vertical displacement changes at different depths caused by groundwatermining in NCP (a) Surface, blue triangles are GPS sites; (b) The depth of 5 km; (c) The depth of 15 km.

根据模型计算结果，华北平原地下水开采引起的地壳隆升主要发生于邢台、石家庄和北京等地，平均位移变化速度约0.09 mm·a-1，最大0.6 mm·a-1.GPS和GRACE联合观测数据(Liu et al., 2014)显示华北地区四个基岩站点(BJFS, BJSH, JIXN, TAIN)(图3a)的抬升速度在0.4～2.0 mm·a-1.计算结果相较于观测数据较小的原因可能是现有P、S波速度结构模型对近地表分辨率不够，相关研究表明(沈伟森等, 2010)首都圈近地表(100～500 m)的平均S波速度可能低至300～800 m·s-1，较低的S波速度意味着更大的地表位移(Chen et al., 2011).实际上，地下水开采导致含水层松散沉积物孔隙压实，其造成的地表沉降可能会抵消地下水卸载引起的地壳隆起.因此华北地区大部分建立在沉积层而非基岩上的GPS站点无法观测到这种垂向抬升(赵斌等, 2014; 吕健和杨超, 2015).

2.2 华北平原地下水开采引起的应力场变化

华北地区历史地震目录显示，其震源深度集中于地下5～20 km，平均约10 km(石富强等, 2020).图4给出了华北平原地下水开采在地下10 km深度处引起的应力变化结果.从图中可以看出，地下水开采主要引起垂直方向上的正应力Δσzz变化，变化明显的区域主要为石家庄、邢台等地，最大92 kPa(图4c).水平方向上的正应力变化Δσxx和Δσyy结果相似(其中x指东，y指北)，最大分别为18 kPa和24 kPa(图4a,b)，且后者略大于前者，这可能与开采区整体上沿山前平原呈NNE向分布有关.水平方向上的剪应力变化Δσxy结果最小，仅在-3～3 kPa范围内(图4d).垂直方向上的剪应力变化Δσyz和Δσzx结果分别在-12～10 kPa和-8～16 kPa范围内(图4e,f).

图4 华北平原地下水开采引起的地下10 km处的应力变化Fig.4 Stress changes at -10 km caused by groundwater mining in NCP

地震发生是断层面上积累应力的释放过程，对于大型地震，构造应力是最主要的驱动因素.研究区域位于古老的华北克拉通(朱日祥等, 2012)之上，受太平洋板块和菲律宾板块的联合俯冲作用，区域构造应力场表现为NE-SW向压缩及NW-SE向拉张，最大和最小主应力方向接近水平，主要活动断层为高倾角走滑断层，构造应力积累速率约0.5 kPa·a-1(柳畅等, 2012; 张培震等, 2013).计算结果显示：地下水开采导致地下10 km深度处水平正应力增加约20 kPa，相当于该地区40年的构造应力积累.可见华北平原地下水开采对该地区的应力扰动不可忽略，一定程度上干扰了构造应力加载进程.Pang等(2020)建立黏弹性有限元模型也对华北平原地下水开采引起的应力变化进行了计算，其水平方向上的正应力及剪应力变化结果与本文相近，但垂直方向上正应力变化存在差异，黏弹性模型计算结果最大约70 kPa，相对本文结果偏小.考虑地壳黏弹性性质时，下地壳较低的黏滞系数会使得地壳上部呈现拉张，下部呈现挤压的状态.由于华北地区地壳黏滞系数普遍大于1021Pa·s(Pang et al., 2020)，地下水开采的几十年时间内松弛作用很小，因此本文研究中没有考虑地壳的黏性性质.

2.3 华北平原地下水开采引起的孔隙压力变化

流体抽取引起孔隙压力减小，图5给出了华北平原地下水开采导致的地下10 km处孔隙压力变化随时间的扩散分布.整体上，华北地区孔隙压力随时间减小，且空间上与地下水水位变化分布对应.在开采初期孔隙压力减小发生于石家庄、邢台和北京等水位下降明显的地区，最大减小10 kPa(图5a)；随着地下水不断开采，孔隙压力逐渐向周围扩散，2000年孔隙压力最大减小37 kPa(图5b)；2016年达56 kPa(图5c)，整个区域平均减小约8.1 kPa.断层相对岩体的扩散系数较大，孔隙压力扩散较快，大小和扩散速度取决于水位下降幅度及周围介质的扩散系数.计算结果显示该深度处新河断裂(XHF)的孔隙压力变化最大，约减小36 kPa.

图5 华北平原地下水开采引起的地下10 km深度处的孔隙压力变化(a) 1975年; (b) 2000年; (c) 2016年.Fig.5 Pore pressure changes at -10 km caused by groundwater mining in NCP(a) 1975; (b) 2000; (c) 2016.

2.4 华北平原地下水开采引起的库仑应力变化

华北地区发育有一系列破坏性活动断裂，方向多为NNE和NWW(方菲, 2020)，且大部分为高倾角走滑断层.应力和孔隙压力变化结果显示：该地区地下水卸载对应力场产生明显扰动，而开采区内的活动断裂，如新河断裂(XHF)、唐山—宁河断裂(TNF)和夏垫—凤河营断裂(XFF)等均位于应力变化明显的区域.为了定量分析华北平原地下水开采对区域地震活动性的影响，这里选取1679年三河—平谷8.0级、1937年菏泽7.0级、1966年邢台7.2级和1976年唐山7.8级地震的震源机制解作为库仑应力变化的计算参数，如表2(Pang et al., 2020; 方菲, 2020; 石富强等, 2020).

表2 华北地区强震震源参数Table 2 Source parameter of strong earthquakes in North China

图6给出了四个震源参数计算得到的地下10 km深度处的库仑应力变化结果(摩擦系数μ选取为0.6).

图6 华北平原地下水开采引起的地下10 km深度处的库仑应力变化(a) 三河—平谷地震,考虑孔压; (b) 菏泽地震,考虑孔压; (c) 邢台地震,考虑孔压; (d) 唐山地震,考虑孔压; (e) 三河—平谷地震,忽略孔压; (f) 菏泽地震,忽略孔压; (g) 邢台地震,忽略孔压; (h) 唐山地震,忽略孔压.Fig.6 Coulomb stress changes at -10 km caused by groundwater mining in NCP (a) Sanhe-Pinggu earthquake, consider pore pressure; (b) Heze earthquake, consider pore pressure; (c) Xingtai earthquake, consider pore pressure; (d) Tangshan earthquake, consider pore pressure; (e) Sanhe-Pinggu earthquake, ignore pore pressure; (f) Heze earthquake, ignore pore pressure; (g) Xingtai earthquake, ignore pore pressure; (h) Tangshan earthquake, ignore pore pressure.

为了分析孔隙压力的影响，这里分别给出考虑和忽略孔隙压力时的计算结果.从图中可以看出，四种参数的库仑应力变化结果在空间分布上基本相同，即华北平原地下水开采对该地区活动断层的影响相似.由于地下水开采主要沿NNE向分布，相较于NW向断层，NE向断层的结果稍大.

孔隙压力变化对库仑应力变化计算结果影响较大，若忽略孔隙压力变化(图6e—h)，整个华北地区的库仑应力增加，最大约16 kPa.考虑孔隙压力时(图6a—d)，地下水开采区内的库仑应力减小，最大减小20 kPa，而“开采漏斗区”边缘的库仑应力变化结果为正，尤其是石家庄西北区域.也就是说，地下水开采区内孔隙压力减小起主导作用，结果为负，断层滑动可能性降低；而开采区外围的弹性卸载作用要大于孔隙压力的影响，使库仑应力增加了约3 kPa，可能会促进该区域地震活动的发生.Amos等(2014)对San Joaquin Valley地区的计算结果显示地下水开采使河谷西部Coslinga逆冲断层系的库仑应力增加了1.0～1.7 kPa，与我们考虑孔隙压力变化时的计算结果较接近.

华北地区地震普遍发生于地下5～20 km，但部分地震仍发生于中下地壳，孕震深度可达25 km(Dong et al., 2018).通过本次计算分析，我们可以给出一定的解释：孔隙压力变化随深度衰减较快，会存在地壳深处的库仑应力变化相对浅部为正的值更大，意味着更危险.例如，地下10 km处的库仑应力变化范围在-20～4 kPa，而25 km处范围约在-7～7 kPa.并且，对于临近破裂的断层，kPa量级的应力扰动就可以触发地震(Johnson et al., 2017)，叠加深部岩石圈流变作用对应力积累的影响，华北平原地下水开采也可能对深部断层破裂产生一定的影响.

3 讨论

3.1 华北平原地下水开采对区域地震活动性的影响

地下水开采对区域应力状态的改变是一个复杂的过程，涉及区域构造应力场、岩体岩性及构造分布等因素，所有影响联合控制了地震在何时、何地发生.虽然物理计算模型经过抽象简化，但本文工作的意义在于考虑了区域地质背景和水位变化，定量计算了华北平原地下水开采造成的区域孔隙压力变化、应力变化及库仑应力变化，并探讨地下水开采对区域地震孕育过程可能的影响.

根据库仑应力变化计算结果，若仅考虑水体卸载作用引起的应力改变，华北平原地下水开采使得区域内库仑应力变化结果为正(图6e—f)，增大断层滑动可能性.然而，地下水开采还会引起区域孔隙压力减小并随时间向周围扩散，开采区内库仑应力变化结果为负，开采区边缘呈现3 kPa的增加(图6a—d).Johnson等(2017)根据GPS数据计算了加州断层面上的水文载荷及由此产生的应力变化，得到走滑型地震的剪应力幅度在-1.1～1.6 kPa，其结果增大与地震数量增加存在相关性，表明微震活动的周期性变化受水文载荷影响.同时，他们计算了5.5级以上地震事件的库仑应力变化(-0.6～0.6 kPa)，虽然应力扰动幅度较小，但仍认为水体载荷变化会促进地震成核并调节地震活动性.尤其对临近破裂的断层，低幅度的应力扰动也有触发大地震的可能性(King et al., 1994).历史地震目录显示，华北地区除1966年邢台和1976年唐山地震及其余震外，地下水大规模开采以来该地区以微、小震为主(Pang et al., 2020)，且主要沿山前平原分布，如图7所示.两次强震均发生于20世纪70年代，释放了区域内积累的构造应力.随着地下水大量开采，区域内孔隙压力持续减小，一定程度上降低了该地区强震发生可能性.从该地区地震数目来看，近几十年来微小型地震有所增加，尤其是北京地区.由于缺乏完整的地下水开采前的地震目录，不能证明地震数目的增加是由地下水开采引起的，但开采引起的弹性卸载可能对部分微震的发生存在影响.在此研究中，我们发现华北平原地下水开采引起的孔隙压力减小会抑制开采区内断层滑动，延缓强震发生.由于孔隙压力扩散效应的滞后，唐山地震发生时地下水开采引起的库仑应力变化基本为0，因此开采初期的强震活动可能更多地受构造应力的影响；而开采区外断层面上的孔隙压力减小要小于卸载作用引起的正应力增加，可能促进部分微震发生.

图7 1970—2017年华北地区3.0级以上地震活动性分布Fig.7 Seismicity distribution with M≥ 3.0 in North China from 1970 to 2017

地下水开采、水库和水压致裂等人类活动都会对区域构造应力加载过程造成一定干扰，从而影响区域地震活动性.地下水开采与水库蓄/放水引起的地壳加/卸载相比较，后者较快的水位变化使得其影响幅度相对较大，大部分库区下方的库仑应力变化可达到50～70 kPa(雷兴林等, 2008; 孙玉军等, 2012; 程惠红等, 2012)，但库仑应力变化随深度和水平距离的增大而迅速减小，影响区域范围较小，主要集中于近地表的库区附近.华北平原地下水开采则不同于水库蓄水，它对区域应力场影响幅度相对较小，但地下水开采区面积往往是水库面积的上百倍，地下水水位变化的时滞效应会大一些.

此外，为了缓解华北地区水资源问题的制约，2014年9月国家完成并正式实施“南水北调“工程(吴海峰, 2016).该工程大大缓解了华北地区水资源严重缺乏的情况，华北平原地下水损失量不断降低，浅层地下水水位下降速度减小，甚至有不同程度的回升(崔亚莉等, 2009).地下水水位下降速度的放缓减轻了水体卸载作用的影响，缓解其对构造应力积累的干扰.然而，随着部分地区地下水回灌，地下水水位恢复，区域内弹性荷载和孔隙压力将会增加，在一定程度上将会影响区域断层的活动性(Hemami et al.,2021).Pang等(2020)仅考虑地壳的黏弹性性质，也发现若地下水水位恢复会促进一些地震活动发生.

3.2 摩擦系数对库仑应力变化计算的影响

近年来，库仑应力变化被广泛应用于研究地震触发问题，特别是在初始构造应力场不明确的情况下，库仑应力变化是分析地震发生时间-地点以及主震与余震关系等应力触发问题的有效方法(King et al., 1994).由式(3)可以看出，库仑应力变化计算结果不仅受到断层几何形状和滑动参数的影响，还受到介质摩擦系数(μ)的影响.为了进一步分析μ的影响，我们以1976年唐山7.8级地震的震源机制解(表2)为例，分别选取μ为0.4和0.6进行计算，以分析其对结果的影响，如图8所示.可以看出，μ值虽然没有改变库仑应力变化的极性，但对计算结果值影响较大，最大差别约为6 kPa(图8c).实际上，μ并不是一个定值，地震发生前后会发生变化.若忽略这一变化，会出现库仑应力随摩擦系数增大而增加的矛盾(朱守彪和缪淼, 2016).目前，构造应力场的影响常常被忽略，使得结果的解释与实际情况出现较大偏差，这也是目前库仑应力变化研究触发问题的缺点之一，值得进一步深入研究.

图8 不同摩擦系数下的库仑应力变化及差值(a) μ=0.4; (b) μ=0.6; (c) 两者结果差值.Fig.8 Coulomb stress changes and difference under different friction coefficients(a) μ=0.4; (b) μ=0.6; (c) The difference between the two results.

3.3 不同断层类型对地下水开采的响应

华北地区主要断层均为高倾角走滑型，计算结果显示地下水开采引起开采区内库仑应力变化结果为负，会使得断层破裂危险性降低.但是，不同断层类型会对水体卸载作用产生不同响应，为此我们进一步分析地下水开采对典型正断层和逆断层的影响，如图9所示.图9a—c和图9d—f分别以纯逆冲断层(走向45.0°，倾角30.0°，滑动角90°)和正断层(走向45.0°，倾角60.0°，滑动角270°)参数计算了库仑应力变化、正应力和剪应力结果.结果显示：1)对于正应力Δσn，无论逆冲断层还是正断层，地下水开采都会引起断层面上的Δσn增加(图9b,e)，变化大小受倾角控制，倾角30.0°和60.0°对应的最大值分别为75 kPa和40 kPa.考虑到开采引起孔隙压力减小，断层面上的有效正应力会随倾角增大而逐渐减小；2)对于剪应力Δτ，地下水开采会使逆冲断层的Δτ增加(图9c)，正断层的Δτ减小(图9f)，倾角在45°时达到最值；3)地下水开采使低角度逆冲断层面上的有效正应力和剪应力都增大，库仑应力变化为正(图9a)，地震发生可能性增加.而对于正断层，开采区内孔隙压力减小占主导，断层面上的有效正应力减小，库仑应力变化为负；开采区以外，尤其是开采区边缘，正应力变化可能大于孔隙压力减小使得库仑应力变化出现正值，断层滑动可能性增加.

图9 不同断层类型对地下水开采的响应(a—c)以纯逆冲断层(走向45.0°,倾角30.0°,滑动角90°)为计算参数; (d—f)以正断层(走向45.0°,倾角60.0°,滑动角270°)为计算参数.Fig.9 Responses of different fault types to groundwater mining(a—c) Take the thrust faults (strike angle 45.0°, dip angle 30.0°, slip angle 90°) as the calculating parameters; (d—f) Take the normal fault (strike angle 45.0°, dip angle 60.0°, slip angle 270°) as the calculating parameters.

总之，地下水开采是对地壳的卸载过程，会使断层面上的压性正应力减小，增大逆断层滑动可能性，尤其是对低角度逆冲断层(González et al., 2012; Kundu et al., 2015)；而对于有一定角度的正断层或走滑断层，其断层面上的库仑应力变化结果受到断层角度及断层位置的影响.开采区内，断层孔隙压力减小，断层面上的有效正应力取决于于断层倾角；开采区外，弹性卸载的应力增加大于孔隙压力减小，断层可能更易滑动.

4 结论

自20世纪60年代以来，华北平原地下水水位经历了一个快速下降过程.在前人已有的研究基础上，本文基于孔隙弹性耦合理论建立了华北地区三维孔隙弹性模型，定量计算了华北平原1959—2016年间地下水持续开采引起的区域形变、应力、孔隙压力和库仑应力变化，分析其对区域地震活动性的影响.初步得出：华北平原地下水开采造成的卸载作用会引起区域内地壳抬升，但沉积层孔隙压实导致的地表沉降量要大于抬升量.若GPS站点建立在沉积层而非基岩层上，则观测不到这种隆升；地下水持续开采扰动了区域构造应力的加载过程，震源深度处应力变化可达到几十kPa.根据历史地震震源参数计算得到华北平原地下水开采造成开采区断层面上的有效正应力减小，库仑应力变化为负，一定程度上延缓区域内强震的发生.开采区外弹性骨架的应力变化占主导作用，使库仑应力变化为正，尤其是开采区西侧边缘，可能对部分微震的发生产生影响.

致谢感谢匿名评审专家对文本完善提出的宝贵修改意见和建议.