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

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

中國科學院大學地球與行星科學學院，計算地球動力學重點實驗室， 北京 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.根據歷史地震震源參數計算得到華北平原地下水開采造成開采區斷層面上的有效正應力減小，庫侖應力變化為負，一定程度上延緩區域內強震的發生.開采區外彈性骨架的應力變化占主導作用，使庫侖應力變化為正，尤其是開采區西側邊緣，可能對部分微震的發生產生影響.

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