華北平原地下水開采對區(qū)域地震活動性的影響

竇甜甜, 程惠紅, 周元澤, 石耀霖

中國科學(xué)院大學(xué)地球與行星科學(xué)學(xué)院，計算地球動力學(xué)重點實驗室， 北京 100049

0 引言

研究表明地下水開采、水庫蓄水、頁巖氣開發(fā)、深井注入等對地震發(fā)生存在顯著影響，人類活動誘/觸發(fā)地震已成為地球物理學(xué)研究的一個重要領(lǐng)域(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).對于區(qū)域性地下水開采，大規(guī)模的地殼卸載不僅能直接擾動應(yīng)力和應(yīng)變積累過程(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)，也會引起區(qū)域孔隙壓力變化并以擴散的形式向周圍傳播(Rice and Cleary, 1976; Talwani and Acree, 1984; Kundu et al., 2015; Segall and Lu, 2015; Bhattacharya and Viesca, 2019).Amos等(2014)對地下水開采引起的彈性應(yīng)力變化進行了研究和分析，并認(rèn)為卸載降低了斷層面上的正應(yīng)力.González等(2012)根據(jù)彈性位錯模型計算了Lorca地震的庫侖應(yīng)力變化，認(rèn)為該地震可能是由地下水卸載引起的.Kundu等(2015)通過模擬計算得出2015年Nepal地震以及主喜馬拉雅逆沖斷裂區(qū)發(fā)生的地震均受到印度河—恒河平原地下水超采的影響.Segall和Lu(2015)分別計算了地下水開采引起的彈性應(yīng)力和孔隙壓力變化對地震活動的影響，并強調(diào)了孔彈耦合應(yīng)力的作用.龐亞瑾等(2016)模擬計算并分析了華北平原地下水開采對地殼應(yīng)力的影響，但未考慮孔隙壓力變化.Wetzler等(2019)據(jù)觀測數(shù)據(jù)發(fā)現(xiàn)地下水快速抽取后Kinneret湖區(qū)域發(fā)生兩次異常的淺地震群，通過孔隙彈性模型認(rèn)為與地下水開采有關(guān).根據(jù)孔隙彈性耦合理論，孔隙壓力變化會引起應(yīng)力場改變，同時孔隙壓力也受到平均正應(yīng)力的影響，兩者相互作用共同控制斷層的穩(wěn)定性(Rice and Cleary, 1976; 程惠紅等, 2012; Segall and Lu, 2015; Wetzler et al., 2019; Hemami et al., 2021; Barbour and Beeler, 2021).

目前，地下水仍是人們生活和灌溉的主要來源，特別是我國人口相對密集的華北平原，自20世紀(jì)60年代以來該地區(qū)地下水存在明顯的長期虧損，山前平原區(qū)(包括石家莊、邢臺等地)淺層地下水水位下降達30～50 m(費宇紅等, 2009).2000年后，華北地區(qū)全年水供應(yīng)約74%來源于地下水(Zheng et al., 2010).蘇曉莉等(2012)利用GRACE重力觀測數(shù)據(jù)并與相關(guān)水文模型比較得出華北地下水水位下降速度約0.5 cm·a-1.而隨著“南水北調(diào)”工程的實施，近年來該地區(qū)地下水水位下降速度有所變緩(周志博等, 2020).

華北平原作為中國大陸典型強震多發(fā)區(qū)，近幾十年來持續(xù)的地下水開采造成的地殼卸載和孔隙壓力減小可能會對其產(chǎn)生較大影響.因此，本文基于孔隙彈性耦合理論，定量計算1959—2016年間華北平原地下水開采引起的區(qū)域應(yīng)力和孔隙壓力變化，并根據(jù)歷史地震參數(shù)計算庫侖應(yīng)力變化，分析并探討地下水開采對區(qū)域地震活動性的影響.

1 計算方法及模型建立

1.1 華北地區(qū)三維模型

華北平原位于中國大陸東部，西鄰山西斷陷盆地帶，東鄰郯廬斷裂帶(張培震等, 2013).自晚中生代以來該地區(qū)經(jīng)歷太平洋深俯沖和弧后擴展作用，活動斷裂發(fā)育，地震活動強烈，歷史上曾多次發(fā)生8.0級以上的破壞性地震(Liu et al., 2011).近代華北地區(qū)也多次記錄到7.0級以上的強震，如1966年邢臺7.2級地震、1976年唐山7.8級地震(Liu et al., 2011; 王輝等, 2011) (圖1).研究表明華北地區(qū)地震活動存在明顯的“期幕”特征，時間上呈現(xiàn)活動-平靜期交替出現(xiàn)的周期性，空間上表現(xiàn)為地震帶之間的叢集遷移特征(Liu et al., 2011; 尹曉菲等, 2020).

圖1 華北地區(qū)強震分布及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.

中、新生代時期，西太平洋板塊俯沖導(dǎo)致華北克拉通破裂(朱日祥等, 2012; Zheng et al., 2017; Ma et al., 2022)，發(fā)育了一系列活動斷裂，以北東向的高傾角右旋走滑張扭性斷層為主(鄧起東等, 2003; 張培震等, 2013).區(qū)域內(nèi)主要活動斷裂有夏墊—鳳河營斷裂(XFF)(何付兵等, 2013)、唐山—寧河斷裂(TNF)(李志義和虢順民, 1979; 張素欣等, 2020)、河間斷裂(HJF)(楊家亮等, 2010)、新河斷裂(XHF，又稱百尺口斷裂)(徐杰和方仲景, 1988)和湯西斷裂(TXF)(韓慕康和趙景珍, 1980).遠(yuǎn)震波形數(shù)據(jù)及流動臺站地殼觀測數(shù)據(jù)顯示華北地區(qū)地殼厚度平均約40 km(危自根等, 2015; 王椿鏞等, 2017)，其中上地殼10～12 km，中地殼8～10 km，下地殼10～15 km(王椿鏞等, 2017).

地貌上，華北平原屬典型平原地貌，其上覆蓋第四紀(jì)含水層(Foster et al., 2004; 孟素花等, 2013).地貌構(gòu)成以山前沖、洪積傾斜平原，中部沖、湖積多層疊積平原和東部沖、湖積為主夾數(shù)層海積層的濱海平原為主(孟素花等, 2013).根據(jù)水文地質(zhì)特征，可分為山前沖洪積平原、中部沖積湖積平原、古黃河沖洪積平原和東部沖積海積平原，其中山前平原顆粒粗，具有較強的滲透能力，中東部平原巖性逐漸變細(xì)，多為粉土、粉質(zhì)黏土和黏土.

基于以上地質(zhì)背景資料建立了華北地區(qū)三維有限元模型(圖2).模型采用分層分塊結(jié)構(gòu)，深度上包括沉積層及上、中、下地殼在內(nèi)，共40 km，斷層切割深度15 km，地表地形采用SRTM 90 m數(shù)字高程數(shù)據(jù)(http:∥dds.cr.usgs.gov/srtm/).根據(jù)水文地質(zhì)特征，開采區(qū)沉積蓋層可劃分為4個單元.采用三棱柱單元進行網(wǎng)格剖分，整個模型節(jié)點數(shù)為378690，單元數(shù)為698572，如圖2a.

圖2 華北地區(qū)有限元模型(a) 三維地質(zhì)模型; (b) 邊界條件.Fig.2 The finite element model of NCP(a) Three-dimensional geological model; (b) Boundary conditions.

1.2 計算方法

地下水開采不僅對地殼產(chǎn)生卸載作用，同時引起孔隙壓力變化，孔隙彈性耦合理論可以較好地分析流體與固體骨架之間的耦合作用(Talwani and Acree, 1984).根據(jù)Biot(1941)固結(jié)理論，Rice和Cleary(1976)推導(dǎo)得出了孔隙彈性耦合方程，經(jīng)眾多學(xué)者研究和發(fā)展，目前已經(jīng)廣泛應(yīng)用于分析巖體在流體注入/抽取及水體載荷作用下巖石應(yīng)力場和孔隙壓力場的變化(雷興林等, 2008; 程惠紅等, 2012; 孫玉軍等, 2012; Segall and Lu, 2015)，其本構(gòu)關(guān)系為

(1)

(2)

在已知斷層參數(shù)的情況下，庫侖應(yīng)力變化可以由下式計算得到：

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

(3)

其中，Δσn是斷層面上正應(yīng)力變化，Δτ是斷層面上沿滑動方向的剪應(yīng)力變化.由于大多數(shù)情況下缺少研究區(qū)域初始構(gòu)造應(yīng)力場，庫侖應(yīng)力變化成為定量分析地震危險性的常用方法，被廣泛應(yīng)用于地震觸發(fā)問題研究(King et al., 1994; Stein et al., 1997;石耀霖和曹建玲, 2010).由孔隙彈性模型計算得到應(yīng)力張量和孔隙壓力變化量，可以計算特定震源參數(shù)下的庫侖應(yīng)力變化.若結(jié)果為正，表明斷層滑動危險性增加；反之，則表明斷層趨于穩(wěn)定，從而了解分析地震相對更容易發(fā)生于哪些地方，何種斷層.

1.3 計算參數(shù)和邊界條件

根據(jù)華北地區(qū)巖石圈P波和S波速度結(jié)構(gòu)(Pang et al., 2020; 石富強等, 2020)可以給出模型彈性參數(shù).相關(guān)研究表明(Talwani et al., 2007)孔隙彈性介質(zhì)的擴散系數(shù)范圍一般在0.1～10 m2·s-1，本文基于華北平原沉積蓋層的水文地質(zhì)特征給定計算所用的擴散系數(shù).模型具體計算參數(shù)見表1，計算總時間為1959—2016年(58年)，計算時間步長為1年.

表1 模型計算參數(shù)Table 1 Model calculation parameters

由于華北平原地下水開采主要集中于上部含水層，其深度相對模型深度較淺，且在開采時間尺度內(nèi)可忽略地幔松弛效應(yīng)而僅考慮彈性作用.因此，地下水的卸載可近似為對模型上邊界施加垂直向上的法向力.模型彈性位移場邊界條件設(shè)定為：底邊界、側(cè)邊界法向固定，切向自由，上邊界自由，其中開采區(qū)內(nèi)施加垂直向上的拉力，大小為地下水位變化、孔隙度和水的容重之積；孔隙壓力場設(shè)定為：底邊界、側(cè)邊界滲流梯度為0，上邊界由水體卸載量給定，如圖2b.相關(guān)研究表明華北地區(qū)含水層孔隙度范圍在0.03～0.21(Pang et al., 2020)，本文計算中設(shè)模型具有均勻孔隙度，大小為0.18.

2 計算結(jié)果

2.1 華北平原地下水開采引起的位移場變化

地下水開采對地殼產(chǎn)生卸載作用，其引起的位移變化主要受彈性卸載控制，空間上與地下水水位下降分布相對應(yīng).圖3給出了華北平原地下水開采引起的不同深度上的垂直位移變化.由于水體卸載，華北地區(qū)不同深度上都呈現(xiàn)出一定抬升，最大抬升出現(xiàn)在水位下降最明顯的邢臺、石家莊附近.開采區(qū)內(nèi)地表平均抬升量大于5 mm，最大抬升35 mm，位于邢臺以北(圖3a).位移變化結(jié)果隨深度增加而減小，上地殼上表面(地下5 km)和中地殼上表面(地下15 km)最大垂向位移分別為28 mm(圖3b)和18 mm(圖3c).相較于垂直方向上的位移變化，水平方向上位移變化較小，結(jié)果在5 mm以下.開采區(qū)以外的位移場變化不大，基本可忽略.

圖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.

根據(jù)模型計算結(jié)果，華北平原地下水開采引起的地殼隆升主要發(fā)生于邢臺、石家莊和北京等地，平均位移變化速度約0.09 mm·a-1，最大0.6 mm·a-1.GPS和GRACE聯(lián)合觀測數(shù)據(jù)(Liu et al., 2014)顯示華北地區(qū)四個基巖站點(BJFS, BJSH, JIXN, TAIN)(圖3a)的抬升速度在0.4～2.0 mm·a-1.計算結(jié)果相較于觀測數(shù)據(jù)較小的原因可能是現(xiàn)有P、S波速度結(jié)構(gòu)模型對近地表分辨率不夠，相關(guān)研究表明(沈偉森等, 2010)首都圈近地表(100～500 m)的平均S波速度可能低至300～800 m·s-1，較低的S波速度意味著更大的地表位移(Chen et al., 2011).實際上，地下水開采導(dǎo)致含水層松散沉積物孔隙壓實，其造成的地表沉降可能會抵消地下水卸載引起的地殼隆起.因此華北地區(qū)大部分建立在沉積層而非基巖上的GPS站點無法觀測到這種垂向抬升(趙斌等, 2014; 呂健和楊超, 2015).

2.2 華北平原地下水開采引起的應(yīng)力場變化

華北地區(qū)歷史地震目錄顯示，其震源深度集中于地下5～20 km，平均約10 km(石富強等, 2020).圖4給出了華北平原地下水開采在地下10 km深度處引起的應(yīng)力變化結(jié)果.從圖中可以看出，地下水開采主要引起垂直方向上的正應(yīng)力Δσzz變化，變化明顯的區(qū)域主要為石家莊、邢臺等地，最大92 kPa(圖4c).水平方向上的正應(yīng)力變化Δσxx和Δσyy結(jié)果相似(其中x指東，y指北)，最大分別為18 kPa和24 kPa(圖4a,b)，且后者略大于前者，這可能與開采區(qū)整體上沿山前平原呈NNE向分布有關(guān).水平方向上的剪應(yīng)力變化Δσxy結(jié)果最小，僅在-3～3 kPa范圍內(nèi)(圖4d).垂直方向上的剪應(yīng)力變化Δσyz和Δσzx結(jié)果分別在-12～10 kPa和-8～16 kPa范圍內(nèi)(圖4e,f).

圖4 華北平原地下水開采引起的地下10 km處的應(yīng)力變化Fig.4 Stress changes at -10 km caused by groundwater mining in NCP

地震發(fā)生是斷層面上積累應(yīng)力的釋放過程，對于大型地震，構(gòu)造應(yīng)力是最主要的驅(qū)動因素.研究區(qū)域位于古老的華北克拉通(朱日祥等, 2012)之上，受太平洋板塊和菲律賓板塊的聯(lián)合俯沖作用，區(qū)域構(gòu)造應(yīng)力場表現(xiàn)為NE-SW向壓縮及NW-SE向拉張，最大和最小主應(yīng)力方向接近水平，主要活動斷層為高傾角走滑斷層，構(gòu)造應(yīng)力積累速率約0.5 kPa·a-1(柳暢等, 2012; 張培震等, 2013).計算結(jié)果顯示：地下水開采導(dǎo)致地下10 km深度處水平正應(yīng)力增加約20 kPa，相當(dāng)于該地區(qū)40年的構(gòu)造應(yīng)力積累.可見華北平原地下水開采對該地區(qū)的應(yīng)力擾動不可忽略，一定程度上干擾了構(gòu)造應(yīng)力加載進程.Pang等(2020)建立黏彈性有限元模型也對華北平原地下水開采引起的應(yīng)力變化進行了計算，其水平方向上的正應(yīng)力及剪應(yīng)力變化結(jié)果與本文相近，但垂直方向上正應(yīng)力變化存在差異，黏彈性模型計算結(jié)果最大約70 kPa，相對本文結(jié)果偏小.考慮地殼黏彈性性質(zhì)時，下地殼較低的黏滯系數(shù)會使得地殼上部呈現(xiàn)拉張，下部呈現(xiàn)擠壓的狀態(tài).由于華北地區(qū)地殼黏滯系數(shù)普遍大于1021Pa·s(Pang et al., 2020)，地下水開采的幾十年時間內(nèi)松弛作用很小，因此本文研究中沒有考慮地殼的黏性性質(zhì).

2.3 華北平原地下水開采引起的孔隙壓力變化

流體抽取引起孔隙壓力減小，圖5給出了華北平原地下水開采導(dǎo)致的地下10 km處孔隙壓力變化隨時間的擴散分布.整體上，華北地區(qū)孔隙壓力隨時間減小，且空間上與地下水水位變化分布對應(yīng).在開采初期孔隙壓力減小發(fā)生于石家莊、邢臺和北京等水位下降明顯的地區(qū)，最大減小10 kPa(圖5a)；隨著地下水不斷開采，孔隙壓力逐漸向周圍擴散，2000年孔隙壓力最大減小37 kPa(圖5b)；2016年達56 kPa(圖5c)，整個區(qū)域平均減小約8.1 kPa.斷層相對巖體的擴散系數(shù)較大，孔隙壓力擴散較快，大小和擴散速度取決于水位下降幅度及周圍介質(zhì)的擴散系數(shù).計算結(jié)果顯示該深度處新河斷裂(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 華北平原地下水開采引起的庫侖應(yīng)力變化

華北地區(qū)發(fā)育有一系列破壞性活動斷裂，方向多為NNE和NWW(方菲, 2020)，且大部分為高傾角走滑斷層.應(yīng)力和孔隙壓力變化結(jié)果顯示：該地區(qū)地下水卸載對應(yīng)力場產(chǎn)生明顯擾動，而開采區(qū)內(nèi)的活動斷裂，如新河斷裂(XHF)、唐山—寧河斷裂(TNF)和夏墊—鳳河營斷裂(XFF)等均位于應(yīng)力變化明顯的區(qū)域.為了定量分析華北平原地下水開采對區(qū)域地震活動性的影響，這里選取1679年三河—平谷8.0級、1937年菏澤7.0級、1966年邢臺7.2級和1976年唐山7.8級地震的震源機制解作為庫侖應(yīng)力變化的計算參數(shù)，如表2(Pang et al., 2020; 方菲, 2020; 石富強等, 2020).

表2 華北地區(qū)強震震源參數(shù)Table 2 Source parameter of strong earthquakes in North China

圖6給出了四個震源參數(shù)計算得到的地下10 km深度處的庫侖應(yīng)力變化結(jié)果(摩擦系數(shù)μ選取為0.6).

圖6 華北平原地下水開采引起的地下10 km深度處的庫侖應(yīng)力變化(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.

為了分析孔隙壓力的影響，這里分別給出考慮和忽略孔隙壓力時的計算結(jié)果.從圖中可以看出，四種參數(shù)的庫侖應(yīng)力變化結(jié)果在空間分布上基本相同，即華北平原地下水開采對該地區(qū)活動斷層的影響相似.由于地下水開采主要沿NNE向分布，相較于NW向斷層，NE向斷層的結(jié)果稍大.

孔隙壓力變化對庫侖應(yīng)力變化計算結(jié)果影響較大，若忽略孔隙壓力變化(圖6e—h)，整個華北地區(qū)的庫侖應(yīng)力增加，最大約16 kPa.考慮孔隙壓力時(圖6a—d)，地下水開采區(qū)內(nèi)的庫侖應(yīng)力減小，最大減小20 kPa，而“開采漏斗區(qū)”邊緣的庫侖應(yīng)力變化結(jié)果為正，尤其是石家莊西北區(qū)域.也就是說，地下水開采區(qū)內(nèi)孔隙壓力減小起主導(dǎo)作用，結(jié)果為負(fù)，斷層滑動可能性降低；而開采區(qū)外圍的彈性卸載作用要大于孔隙壓力的影響，使庫侖應(yīng)力增加了約3 kPa，可能會促進該區(qū)域地震活動的發(fā)生.Amos等(2014)對San Joaquin Valley地區(qū)的計算結(jié)果顯示地下水開采使河谷西部Coslinga逆沖斷層系的庫侖應(yīng)力增加了1.0～1.7 kPa，與我們考慮孔隙壓力變化時的計算結(jié)果較接近.

華北地區(qū)地震普遍發(fā)生于地下5～20 km，但部分地震仍發(fā)生于中下地殼，孕震深度可達25 km(Dong et al., 2018).通過本次計算分析，我們可以給出一定的解釋：孔隙壓力變化隨深度衰減較快，會存在地殼深處的庫侖應(yīng)力變化相對淺部為正的值更大，意味著更危險.例如，地下10 km處的庫侖應(yīng)力變化范圍在-20～4 kPa，而25 km處范圍約在-7～7 kPa.并且，對于臨近破裂的斷層，kPa量級的應(yīng)力擾動就可以觸發(fā)地震(Johnson et al., 2017)，疊加深部巖石圈流變作用對應(yīng)力積累的影響，華北平原地下水開采也可能對深部斷層破裂產(chǎn)生一定的影響.

3 討論

3.1 華北平原地下水開采對區(qū)域地震活動性的影響

地下水開采對區(qū)域應(yīng)力狀態(tài)的改變是一個復(fù)雜的過程，涉及區(qū)域構(gòu)造應(yīng)力場、巖體巖性及構(gòu)造分布等因素，所有影響聯(lián)合控制了地震在何時、何地發(fā)生.雖然物理計算模型經(jīng)過抽象簡化，但本文工作的意義在于考慮了區(qū)域地質(zhì)背景和水位變化，定量計算了華北平原地下水開采造成的區(qū)域孔隙壓力變化、應(yīng)力變化及庫侖應(yīng)力變化，并探討地下水開采對區(qū)域地震孕育過程可能的影響.

根據(jù)庫侖應(yīng)力變化計算結(jié)果，若僅考慮水體卸載作用引起的應(yīng)力改變，華北平原地下水開采使得區(qū)域內(nèi)庫侖應(yīng)力變化結(jié)果為正(圖6e—f)，增大斷層滑動可能性.然而，地下水開采還會引起區(qū)域孔隙壓力減小并隨時間向周圍擴散，開采區(qū)內(nèi)庫侖應(yīng)力變化結(jié)果為負(fù)，開采區(qū)邊緣呈現(xiàn)3 kPa的增加(圖6a—d).Johnson等(2017)根據(jù)GPS數(shù)據(jù)計算了加州斷層面上的水文載荷及由此產(chǎn)生的應(yīng)力變化，得到走滑型地震的剪應(yīng)力幅度在-1.1～1.6 kPa，其結(jié)果增大與地震數(shù)量增加存在相關(guān)性，表明微震活動的周期性變化受水文載荷影響.同時，他們計算了5.5級以上地震事件的庫侖應(yīng)力變化(-0.6～0.6 kPa)，雖然應(yīng)力擾動幅度較小，但仍認(rèn)為水體載荷變化會促進地震成核并調(diào)節(jié)地震活動性.尤其對臨近破裂的斷層，低幅度的應(yīng)力擾動也有觸發(fā)大地震的可能性(King et al., 1994).歷史地震目錄顯示，華北地區(qū)除1966年邢臺和1976年唐山地震及其余震外，地下水大規(guī)模開采以來該地區(qū)以微、小震為主(Pang et al., 2020)，且主要沿山前平原分布，如圖7所示.兩次強震均發(fā)生于20世紀(jì)70年代，釋放了區(qū)域內(nèi)積累的構(gòu)造應(yīng)力.隨著地下水大量開采，區(qū)域內(nèi)孔隙壓力持續(xù)減小，一定程度上降低了該地區(qū)強震發(fā)生可能性.從該地區(qū)地震數(shù)目來看，近幾十年來微小型地震有所增加，尤其是北京地區(qū).由于缺乏完整的地下水開采前的地震目錄，不能證明地震數(shù)目的增加是由地下水開采引起的，但開采引起的彈性卸載可能對部分微震的發(fā)生存在影響.在此研究中，我們發(fā)現(xiàn)華北平原地下水開采引起的孔隙壓力減小會抑制開采區(qū)內(nèi)斷層滑動，延緩強震發(fā)生.由于孔隙壓力擴散效應(yīng)的滯后，唐山地震發(fā)生時地下水開采引起的庫侖應(yīng)力變化基本為0，因此開采初期的強震活動可能更多地受構(gòu)造應(yīng)力的影響；而開采區(qū)外斷層面上的孔隙壓力減小要小于卸載作用引起的正應(yīng)力增加，可能促進部分微震發(fā)生.

圖7 1970—2017年華北地區(qū)3.0級以上地震活動性分布Fig.7 Seismicity distribution with M≥ 3.0 in North China from 1970 to 2017

地下水開采、水庫和水壓致裂等人類活動都會對區(qū)域構(gòu)造應(yīng)力加載過程造成一定干擾，從而影響區(qū)域地震活動性.地下水開采與水庫蓄/放水引起的地殼加/卸載相比較，后者較快的水位變化使得其影響幅度相對較大，大部分庫區(qū)下方的庫侖應(yīng)力變化可達到50～70 kPa(雷興林等, 2008; 孫玉軍等, 2012; 程惠紅等, 2012)，但庫侖應(yīng)力變化隨深度和水平距離的增大而迅速減小，影響區(qū)域范圍較小，主要集中于近地表的庫區(qū)附近.華北平原地下水開采則不同于水庫蓄水，它對區(qū)域應(yīng)力場影響幅度相對較小，但地下水開采區(qū)面積往往是水庫面積的上百倍，地下水水位變化的時滯效應(yīng)會大一些.

此外，為了緩解華北地區(qū)水資源問題的制約，2014年9月國家完成并正式實施“南水北調(diào)“工程(吳海峰, 2016).該工程大大緩解了華北地區(qū)水資源嚴(yán)重缺乏的情況，華北平原地下水損失量不斷降低，淺層地下水水位下降速度減小，甚至有不同程度的回升(崔亞莉等, 2009).地下水水位下降速度的放緩減輕了水體卸載作用的影響，緩解其對構(gòu)造應(yīng)力積累的干擾.然而，隨著部分地區(qū)地下水回灌，地下水水位恢復(fù)，區(qū)域內(nèi)彈性荷載和孔隙壓力將會增加，在一定程度上將會影響區(qū)域斷層的活動性(Hemami et al.,2021).Pang等(2020)僅考慮地殼的黏彈性性質(zhì)，也發(fā)現(xiàn)若地下水水位恢復(fù)會促進一些地震活動發(fā)生.

3.2 摩擦系數(shù)對庫侖應(yīng)力變化計算的影響

近年來，庫侖應(yīng)力變化被廣泛應(yīng)用于研究地震觸發(fā)問題，特別是在初始構(gòu)造應(yīng)力場不明確的情況下，庫侖應(yīng)力變化是分析地震發(fā)生時間-地點以及主震與余震關(guān)系等應(yīng)力觸發(fā)問題的有效方法(King et al., 1994).由式(3)可以看出，庫侖應(yīng)力變化計算結(jié)果不僅受到斷層幾何形狀和滑動參數(shù)的影響，還受到介質(zhì)摩擦系數(shù)(μ)的影響.為了進一步分析μ的影響，我們以1976年唐山7.8級地震的震源機制解(表2)為例，分別選取μ為0.4和0.6進行計算，以分析其對結(jié)果的影響，如圖8所示.可以看出，μ值雖然沒有改變庫侖應(yīng)力變化的極性，但對計算結(jié)果值影響較大，最大差別約為6 kPa(圖8c).實際上，μ并不是一個定值，地震發(fā)生前后會發(fā)生變化.若忽略這一變化，會出現(xiàn)庫侖應(yīng)力隨摩擦系數(shù)增大而增加的矛盾(朱守彪和繆淼, 2016).目前，構(gòu)造應(yīng)力場的影響常常被忽略，使得結(jié)果的解釋與實際情況出現(xiàn)較大偏差，這也是目前庫侖應(yīng)力變化研究觸發(fā)問題的缺點之一，值得進一步深入研究.

圖8 不同摩擦系數(shù)下的庫侖應(yīng)力變化及差值(a) μ=0.4; (b) μ=0.6; (c) 兩者結(jié)果差值.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 不同斷層類型對地下水開采的響應(yīng)

華北地區(qū)主要斷層均為高傾角走滑型，計算結(jié)果顯示地下水開采引起開采區(qū)內(nèi)庫侖應(yīng)力變化結(jié)果為負(fù)，會使得斷層破裂危險性降低.但是，不同斷層類型會對水體卸載作用產(chǎn)生不同響應(yīng)，為此我們進一步分析地下水開采對典型正斷層和逆斷層的影響，如圖9所示.圖9a—c和圖9d—f分別以純逆沖斷層(走向45.0°，傾角30.0°，滑動角90°)和正斷層(走向45.0°，傾角60.0°，滑動角270°)參數(shù)計算了庫侖應(yīng)力變化、正應(yīng)力和剪應(yīng)力結(jié)果.結(jié)果顯示：1)對于正應(yīng)力Δσn，無論逆沖斷層還是正斷層，地下水開采都會引起斷層面上的Δσn增加(圖9b,e)，變化大小受傾角控制，傾角30.0°和60.0°對應(yīng)的最大值分別為75 kPa和40 kPa.考慮到開采引起孔隙壓力減小，斷層面上的有效正應(yīng)力會隨傾角增大而逐漸減??；2)對于剪應(yīng)力Δτ，地下水開采會使逆沖斷層的Δτ增加(圖9c)，正斷層的Δτ減小(圖9f)，傾角在45°時達到最值；3)地下水開采使低角度逆沖斷層面上的有效正應(yīng)力和剪應(yīng)力都增大，庫侖應(yīng)力變化為正(圖9a)，地震發(fā)生可能性增加.而對于正斷層，開采區(qū)內(nèi)孔隙壓力減小占主導(dǎo)，斷層面上的有效正應(yīng)力減小，庫侖應(yīng)力變化為負(fù)；開采區(qū)以外，尤其是開采區(qū)邊緣，正應(yīng)力變化可能大于孔隙壓力減小使得庫侖應(yīng)力變化出現(xiàn)正值，斷層滑動可能性增加.

圖9 不同斷層類型對地下水開采的響應(yīng)(a—c)以純逆沖斷層(走向45.0°,傾角30.0°,滑動角90°)為計算參數(shù); (d—f)以正斷層(走向45.0°,傾角60.0°,滑動角270°)為計算參數(shù).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.

總之，地下水開采是對地殼的卸載過程，會使斷層面上的壓性正應(yīng)力減小，增大逆斷層滑動可能性，尤其是對低角度逆沖斷層(González et al., 2012; Kundu et al., 2015)；而對于有一定角度的正斷層或走滑斷層，其斷層面上的庫侖應(yīng)力變化結(jié)果受到斷層角度及斷層位置的影響.開采區(qū)內(nèi)，斷層孔隙壓力減小，斷層面上的有效正應(yīng)力取決于于斷層傾角；開采區(qū)外，彈性卸載的應(yīng)力增加大于孔隙壓力減小，斷層可能更易滑動.

4 結(jié)論

自20世紀(jì)60年代以來，華北平原地下水水位經(jīng)歷了一個快速下降過程.在前人已有的研究基礎(chǔ)上，本文基于孔隙彈性耦合理論建立了華北地區(qū)三維孔隙彈性模型，定量計算了華北平原1959—2016年間地下水持續(xù)開采引起的區(qū)域形變、應(yīng)力、孔隙壓力和庫侖應(yīng)力變化，分析其對區(qū)域地震活動性的影響.初步得出：華北平原地下水開采造成的卸載作用會引起區(qū)域內(nèi)地殼抬升，但沉積層孔隙壓實導(dǎo)致的地表沉降量要大于抬升量.若GPS站點建立在沉積層而非基巖層上，則觀測不到這種隆升；地下水持續(xù)開采擾動了區(qū)域構(gòu)造應(yīng)力的加載過程，震源深度處應(yīng)力變化可達到幾十kPa.根據(jù)歷史地震震源參數(shù)計算得到華北平原地下水開采造成開采區(qū)斷層面上的有效正應(yīng)力減小，庫侖應(yīng)力變化為負(fù)，一定程度上延緩區(qū)域內(nèi)強震的發(fā)生.開采區(qū)外彈性骨架的應(yīng)力變化占主導(dǎo)作用，使庫侖應(yīng)力變化為正，尤其是開采區(qū)西側(cè)邊緣，可能對部分微震的發(fā)生產(chǎn)生影響.

致謝感謝匿名評審專家對文本完善提出的寶貴修改意見和建議.