華北地區(qū)地下水開采對地殼應力的影響

龐亞瑾， 張懷， 程惠紅， 石耀霖

中國科學院計算地球動力學重點實驗室，中國科學院大學， 北京　100049

華北地區(qū)地下水開采對地殼應力的影響

龐亞瑾， 張懷， 程惠紅， 石耀霖*

中國科學院計算地球動力學重點實驗室，中國科學院大學， 北京100049

摘要近50年來華北地區(qū)遭受持續(xù)大面積過量開采地下水，已形成區(qū)域地下水漏斗、地面沉降、地陷地裂等地質(zhì)災害.然而，地下水的抽取減小了地殼的載荷，造成地殼應力場變化，這一點至今尚未被充分認識.為探索華北地區(qū)地下水超采對地殼應力場的影響，本文建立了二維有限元模型，定量計算地下水超采引起地殼變形和應力場變化.結果表明：華北地區(qū)地下水開采會引起地表抬升達+12.4 cm；漏斗區(qū)上、中地殼的水平拉應力增量分別達到70 kPa和35 kPa；而在地下水開采區(qū)外圍，水平壓應力增量達20 kPa；而華北地區(qū)構造主壓應力積累速率約為0.5 kPa·a-1.通過對比華北地區(qū)1980年前后5級以上地震的分布狀況，本文認為地下水開采對區(qū)域構造應力場的擾動不可忽略, 其卸載過程可能對華北地區(qū)大地震孕震過程存在減緩作用.

關鍵詞華北平原； 地下水開采； 應力場； 地震； 數(shù)值模擬

1引言

地球系統(tǒng)科學是目前地球科學發(fā)展的一個前緣，人們往往更關注水圈、大氣圈及凍土圈等變化速率較明顯的現(xiàn)象，卻較少關注人類活動對地球表層和固體地球的影響.本文將提出華北地下水長期過量開采對地殼應力狀態(tài)的影響問題，希望引起我國地球科學者的關注.

隨著經(jīng)濟、人口的增長，水資源需求不斷增加，地下水過度開采已成為全球性焦點問題.特別是近50年來，華北地區(qū)地下水持續(xù)過量開采引起區(qū)域地下水漏斗廣泛分布.自20世紀70年代以來，華北地區(qū)地下水漏斗加速擴展(Shah et al., 2003; Zhang and Li, 2013)，其中北京、石家莊等城市為主要的地下水沉降中心.地下水位埋深大于10 m的區(qū)域占華北平原總面積的40%，且潛水含水層最大埋深達65 m.目前，華北平原深層地下水水位下降速率高達3 m·a-1，淺層水位下降速率高達0.9 m·a-1(Zheng et al., 2010; 費宇紅等, 2009).

地下水過度開采不但引起水質(zhì)惡化、地面沉降等嚴重環(huán)境地質(zhì)問題(Konikow and Kendy, 2005; Sophocleous, 2002)，同時，對地殼應力場也產(chǎn)生重要的影響.Holzer(1979)指出地下水開采會造成地球表層彈性擴張和地表抬升，長期的地下水開采則可能引起一些構造活動.已有研究表明：冰、雪等水文環(huán)境的變化會對地殼產(chǎn)生加/卸載作用，而這種作用會影響地震活動(Bettinelli et al., 2008; Heki, 2003).González等(2012)通過對Lorca地震斷層滑動的觀測和模擬，得出地下水開采對地殼的卸載作用使斷層淺部庫侖應力增加，控制地震斷層的滑動分布.Amos等(2014)據(jù)GPS觀測數(shù)據(jù)認為San JoaquinValley南部地殼垂向上升速率主要歸因于地下水開采，且卸載作用降低了正應力，促進微震的發(fā)生.

華北平原強震構造發(fā)育，是中國大陸強震多發(fā)區(qū).近100年來華北地區(qū)已發(fā)生多次M>5.0地震，其中邢臺、唐山等M>7.0地震對人民生活造成了嚴重影響，給社會經(jīng)濟帶來了嚴重損失.那么，近50年來華北平原地下水過量開采對地殼的卸載作用有多大？是否可能對華北地震活動性造成影響？本研究建立了華北地區(qū)二維有限元模型，采用黏彈性Maxwell本構關系，定量計算地下水超采引起的地殼應力變化，進一步探討地下水開采對華北平原5級以上地震分布的影響.

2模型建立

2.1數(shù)值模型

華北平原地下水沉降漏斗區(qū)域主要集中在北京、保定、石家莊和邯鄲沿線.參考Stephen和Héctor(2004)給出的1960—2000年間華北平原淺層地下水水位下降分布信息，選取北東—南西向剖面AA′(總長500 km)和北西—南東向剖面BB′(總長200 km)來進行分析研究，見圖1.剖面AA′和BB′分別近似平行和垂直于北京—石家莊主要地下水沉降中心沿線，垂向深度為100 km.為降低固定邊界的影響，將AA′和BB′兩剖面兩端分別延伸500 km和300 km，則AA′剖面計算模型總長1500 km，BB′剖面總長800 km.

圖1　華北地區(qū)1960年以來水位下降圖(Stephen and Héctor, 2004)線AA′,BB′為計算模型剖面.Fig.1　Distribution of groundwater decline since 1960 in North China (Stephen and Héctor, 2004)Line AA′, BB′ show the location of model profile.

Zhang等(2011)利用寬角反射剖面對華北地塊研究，初步得出華北平原巖石圈厚度為70 km.Tian等(2014)通過對跨越華北克拉通的東西向地震反射/折射剖面數(shù)據(jù)分析，認為華北平原地殼平均厚度為30 km.Jia和Zhang(2005)據(jù)華北地區(qū)地震測深剖面數(shù)據(jù)，發(fā)現(xiàn)華北平原上地殼厚度約為20 km，下地殼厚度約為10 km.鑒于上述研究，建立了包括上、中、下地殼及地幔和軟流層頂部總厚度為100 km的華北平原分層模型.模型中上、中、下地殼厚度均為10 km，巖石圈地幔40 km和軟流層30 km，如圖2a.采用三角形網(wǎng)格對模型進行剖分，AA′剖面模型分辨率為1.5 km，節(jié)點總數(shù)為70193，單元總數(shù)為138250；BB′剖面模型分辨率為1.2 km，節(jié)點總數(shù)為70686，單元總數(shù)為139736.

由于計算模型水平寬度選取了地下水開采區(qū)寬度3倍的范圍，且垂向計算深度為100 km.考慮到地表卸載對于足夠深度且足夠遠的范圍影響較小，邊界垂向位移可近似為0，所以模型中兩側(cè)邊界法向位移約束為0，切向自由；底邊界垂向位移約束為0，水平方向自由.考慮地下水開采對地殼產(chǎn)生卸載作用，參考AA′和BB′段的地下水位下降信息，對模型上邊界施加垂直向上的法向力.由于本次研究采用Maxwell黏彈性模型，計算時間步長選取為1年，計算總時間為50年，來研究地下水累計開采50年卸載作用對地殼應力的影響.因此，上邊界的初邊值條件設定為：初始法向應力為0，每個時間步法向應力增量為地下水位年平均變化速率、含水層孔隙度和水的容重之積.兩剖面地下水水位年平均下降速率分布見圖2b—2c(Stephen and Héctor, 2004).

2.2動力學模型

由于地下水持續(xù)開采50年對地殼、巖石圈產(chǎn)生卸載作用，地球介質(zhì)在幾十年的黏性松弛效應需要考慮，本文采用黏彈性Maxwell本構方程，計算總時間為50年，時間步長為1年.

黏彈性Maxwell體滿足以下方程形式：

平衡方程

(1)

幾何方程

(2)

Maxwell體應力-應變(本構)關系

(3)

其中：

(4)

(5)

2.3模型參數(shù)選取

前人對中國大陸巖石圈流變性質(zhì)已有大量的研究.例如，Zang等(2005)通過對鄂爾多斯及其鄰近地區(qū)巖石圈地震波速度結構、熱結構和巖石組成等研究，得到鄂爾多斯地塊及其周圍三維巖石圈流變結構；石耀霖和曹建玲(2008)基于巖石圈溫度和應變速率研究成果，利用實驗室流變試驗結果、計算得出中國大陸巖石圈等效黏滯系數(shù)；孫玉軍等(2013)利用地熱和GPS得到的應變速率數(shù)據(jù)采用有限元方法計算得到中國大陸的三維流變結構.基于華北地塊現(xiàn)有的地震反射/折射數(shù)據(jù)得到的巖石圈P波和S波速度結構(Huang et al., 2009; Jia and Zhang, 2005; Tian et al., 2014)，以及Crust2.0模型給出的地殼密度分布(Christensen and Mooney, 1995)，可計算得到巖石圈各分層結構的彈性模量和泊松比等地殼力學參數(shù).表1給出了各分層力學參數(shù).華北平原含水層主要為細砂和粉砂巖，模型選取含水層巖石孔隙度為0.4(弗里澤和徹里，1987；王大純等，2006).

圖2　計算模型及AA′和BB′ 地下水水位變化速率(a) 巖石圈剖面分層模型及邊界條件信息.W為剖面總寬度，AA′剖面W=1500 km，BB′剖面W=800 km; w為上邊界加載寬度，AA′剖面w=500 km，BB′剖面w=200 km; (b)、(c) 加載區(qū)域地下水位年平均下降速率沿兩剖面加載區(qū)分布圖(Stephen and Héctor, 2004).Fig.2　Numerical model and groundwater decline rate on lines AA′ and BB′(a) Layered lithosphere and boundary conditions. W is width of the model, which is 1500 km on line AA′ and 800 km on line BB′, respectively. w is the width of upper boundary which is 500 km on line AA′ and 200 km on line BB′, respectively. (b) and (c) Annual groundwater decline rates on AA′ and BB′, respectively (Stephen and Héctor, 2004).

表1　計算模型采用的力學參數(shù)

3計算結果及討論

3.1華北平原地下水開采引起的區(qū)域位移場變化

圖3給出了AA′和BB′剖面上地下水持續(xù)50年開采后引起的位移場分布.圖3a和3b分別顯示了AA′剖面上垂向和水平向位移分布，可以看出：地下水開采區(qū)地表整體抬升，而外圍區(qū)域的垂向位移較小；開采區(qū)內(nèi)地表垂向位移空間分布變化與水位下降分布相對應，由于保定—邢臺沿線存在大規(guī)模水位下降，該區(qū)地表垂向抬升劇烈，整體大于6 cm.在石家莊附近，地表抬升高達12.4 cm；受地下水開采影響北京地區(qū)地表上升達6 cm.垂向上升量隨深度的增加而減小，在地下20 km深處，保定—邢臺沿線地殼，垂向最大抬升量為7.8 cm.相對垂向位移，地下水開采引起的水平位移較小，地表最大水平位移主要集中在北京附近，約為2.3 cm；而保定—邢臺沿線地表水平位移值均小于0.5 cm.圖3c和3d分別顯示了BB′剖面上垂向和水平向位移分布，同樣地，開采區(qū)地表出現(xiàn)明顯的抬升，抬升最大值約12 cm，集中在石家莊東南附近的地下水最大沉降漏斗區(qū)域；地下20 km深處的最大垂向位移約為9.0 cm；開采區(qū)內(nèi)水平位移大小在0～1.9 cm范圍內(nèi)，地下水超采對開采區(qū)外圍區(qū)域的水平位移影響較小.

圖3　地下水持續(xù)開采50年引起的位移變化(a) AA′剖面垂向位移結果; (b) AA′剖面水平位移結果; (c) BB′剖面垂向位移結果; (d) BB′剖面水平位移結果.Fig.3　Displacement due to groundwater mining in 50 years(a) Vertical displacement on line AA′; (b) Horizontal displacement on line AA′; (c) Vertical displacement on line BB′; (d) Horizontal displacement on line BB′.

巖石圈深部和軟流圈流變性質(zhì)較強，深部對地下水卸載的響應較明顯，因此，地下深部水平位移整體大于地表位移.地表垂向位移空間分布與地下水水位下降分布曲線一致.水位下降越大，對應地殼垂向位移越明顯.由于地下水賦存于松散沉積層中，地下水開采引起的沉積物顆粒壓實、地面沉降現(xiàn)象更為明顯.地表沉積物的壓實量抵消了地下水卸載對應的地殼抬升作用，所以沉積層覆蓋區(qū)域地表無法觀測到垂向上升.從地殼、巖石圈規(guī)?？紤]，地下水開采對地殼、Moho面垂向變形的影響因素不可忽視.

3.2華北平原地下水開采引起的區(qū)域應力場變化

地下水卸載會引起巖石圈應力變化，對區(qū)域構造應力場產(chǎn)生干擾.圖4為地下水持續(xù)開采50年后引起的地殼水平應力變化結果.圖4a為AA′ 剖面上水平應力變化，由圖可知，地下水卸載使開采區(qū)產(chǎn)生水平向拉張變形，拉應力增量高達89 kPa，位于石家莊北東方向地下水沉降漏斗區(qū).邢臺市周圍最大拉應力增量約為80 kPa；保定西南附近區(qū)域拉應力增量達62 kPa，其東北側(cè)區(qū)域受地下水非均勻分布卸載影響，產(chǎn)生局部擠壓變形，壓應力增量高達20 kPa；北京地區(qū)拉應力增量最大為78 kPa.在開采區(qū)的外圍約50 km范圍內(nèi)，地表產(chǎn)生擠壓變形，最大壓應力增量達38 kPa.開采區(qū)50 km以外，地殼應力受地下水影響甚微.由于華北地區(qū)地震震源深度范圍為5～30 km，本文分別選取震源深度范圍內(nèi)5 km、15 km和25 km深度，分析不同深度上水平應力變化及其對構造應力的擾動程度，見圖4b.開采區(qū)地殼深部整體受拉張作用；北京附近區(qū)域，地下5 km深處最大水平拉應力達41 kPa，15 km深最大水平拉應力增量為21 kPa，地下25 km深處拉應力達15 kPa；石家莊附近水平最大拉應力增量在5 km深度高達49 kPa，15 km深度為25 kPa，25 km深度為18 kPa；邢臺地區(qū)地下5 km、15 km和25 km處最大水平拉應力增量分別為36 kPa、25 kPa、18 kPa；保定地區(qū)相對其他地區(qū)，受卸載影響地殼應力變化較小.圖4c為BB′剖面水平應力分布，該剖面上開采區(qū)受地下水卸載作用整體呈現(xiàn)拉張變形，水平拉應力增量高達118 kPa，分布在石家莊東南方向地下水最大漏斗區(qū).開采區(qū)外圍一定范圍同樣產(chǎn)生壓縮變形，壓應力增量高達38 kPa.圖4d為BB′剖面不同深度水平應力變化，該剖面5 km、15 km和25 km深度對應最大水平拉應力增量分別為72 kPa、55 kPa和31 kPa.受地下水水位下降空間分布差異影響，BB′剖面上開采區(qū)應力變化高于AA′剖面.

圖4　地下水持續(xù)開采50年引起的水平正應力分布(a) AA′剖面； (b) AA′剖面不同深度； (c) BB′剖面； (d) BB′剖面不同深度.Fig.4　Horizontal normal stress changes due to groundwater mining in 50 years (a) On line AA′; (b) At different depths on line AA′; (c) On line BB′; (d) At different depths on line BB′.

圖5為受地下水持續(xù)50年開采影響AA′和BB′剖面上垂向應力變化結果.地下水開采區(qū)巖石圈垂向上整體表現(xiàn)為拉張變形，區(qū)域垂向拉應力最大值集中在地下水漏斗中心區(qū)，地表垂向拉應力增量與地下水水位下降呈正比.最大垂向拉應力增量達227 kPa，位于石家莊東南方向最大沉降漏斗中心區(qū)域，邢臺、石家莊地區(qū)垂向拉應力增量高達170 kPa.北京地區(qū)垂向拉應力增量高達145 kPa.垂向應力變化隨深度逐漸減??；開采區(qū)外圍受地下水開采影響不明顯，垂向應力增量小于5 kPa.

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

受太平洋板塊和菲律賓板塊向歐亞大陸的俯沖聯(lián)合作用，華北平原現(xiàn)代構造應力場的主體特征表現(xiàn)為北東東—南西西方向的擠壓和北西—南東方向的拉張作用(謝富仁等, 2004).主壓應力方位為北東至北東東方向，中間主應力基本是垂直的，最大主應力和最小主應力接近水平，處于走滑應力狀態(tài)(陳連旺等, 2001; 崔效鋒和謝富仁, 2001; 鄧起東等, 1979)，華北地區(qū)構造壓應力積累為0.5 kPa·a-1(柳暢等,2012; 朱守彪等,2010)，而本文計算得到華北地區(qū)地下水開采引起的水平拉應力增量在震源深度范圍高達70 kPa.對比區(qū)域構造0.5 kPa·a-1的壓應力增量，地下水開采對華北地區(qū)區(qū)域應力場的擾動不可忽視.近50年來地下水的累積卸載有效降低了華北平原北東—南西方向的主壓應力構造積累，同時又增加了北西—南東向區(qū)域性的拉應力.

圖5　地下水持續(xù)開采50年引起的垂向正應力變化分布(a) AA′剖面； (b) BB′剖面.Fig.5　Vertical normal stress changes due to groundwater mining in 50 years(a) On line AA′; (b) On line BB′.

地震的發(fā)生受現(xiàn)今區(qū)域構造應力狀態(tài)控制(于慎諤等,2000).改革開放以來隨著經(jīng)濟的發(fā)展，水資源需求迅速增加，地下水水位下降速率明顯增加.即近50年來地下水位下降主要集中在1980年以來.參考研究區(qū)1980年前后35年兩個時間段5.0級以上地震分布，見圖6，研究區(qū)內(nèi)1980年之前地震事件數(shù)目明顯多于1980年之后：河北寧晉—束鹿一代，1966—1968年間發(fā)生了近20個5級以上地震，震源深度大于10 km；1980之后，附近區(qū)域僅發(fā)生一次震源深度為20 km的地震.唐山地區(qū)，于1976年發(fā)生7.8級地震，其后4年伴隨有一系列余震；1980年之后，僅在1991年和1995年發(fā)生兩次5.0級地震.由于構造應力是控制地震發(fā)生的主要驅(qū)動力，構造應力積累造成地震發(fā)生的周期性.華北地區(qū)地震構造應力積累受多種因素干擾，2011日本東北M9.0地震對華北地區(qū)壓應力釋放達1.0 kPa，其應力降抵消了5～10年的構造應力積累(張貝等, 2015).而華北地區(qū)地下水幾十年內(nèi)的累積超采引起震源深度范圍壓應力釋放達70 kPa，有效降低了區(qū)域構造應力積累.由于地震復發(fā)周期具有不確定性，通過20世紀80年代前后地震空間分布的差異，本文不能斷定地下水卸載減少了華北地震的發(fā)生.但考慮到卸載引起的地殼應力變化，本文認為地下水開采是干擾構造應力積累的重要因素之一， 地下水開采對巖石圈的卸載過程可能在一定程度上減緩了地震的孕育過程.

圖6　1980年前后35年時間段內(nèi)5級以上地震分布(a) 1945—1979年前地震； (b) 1980—2015年地震.Fig.6　Earthquake (>MS5.0) epicenter distribution of two periods before and after 1980(a) Earthquakes in 1945—1979； (b) Earthquakes in 1980—2015.

4結論

綜上，華北地區(qū)持續(xù)50年的地下水開采對區(qū)域地殼抬升和應力產(chǎn)生重要的影響，同時在一定程度上影響到華北地區(qū)地震活動性.具體表現(xiàn)在以下幾個方面：

(1) 華北地區(qū)地下水過量開采對巖石圈的卸載作用會引起地表微小的水平位移；而地表垂向上升可達12.4 cm，垂向位移隨深度的增加逐漸減小，卸載作用對開采區(qū)以外的區(qū)域垂向位移影響甚微.

(2) 地下水超采引起開采區(qū)整體拉張變形：震源深度處，水平拉應力增加高達70 kPa；對局部區(qū)域產(chǎn)生擠壓作用，壓應力增量達17 kPa.開采外圍出現(xiàn)局部小范圍壓縮變形：震源深度范圍內(nèi)，水平壓應力增量達20 kPa.而華北平原構造壓應力積累約為0.5 kPa·a-1，可見，地下水過量開采整體上有效釋放了華北地區(qū)構造主壓應力積累，可能對華北地區(qū)地震孕震過程有減緩作用.

由于本次模擬采用二維模型模擬計算地下水開采對巖石圈變形的影響，不能精確地給出受地下水卸載影響華北地區(qū)三維地殼應力場及斷層庫侖應力的變化，該工作有待于三維數(shù)值模擬研究解決.

References

Amos C B, Audet P, Hammond W C, et al. 2014. Uplift and seismicity driven by groundwater depletion in central California.Nature, 509(7501): 483-486.

Bettinelli P, Avouac J P, Flouzat M, et al. 2008. Seasonal variations of seismicity and geodetic strain in the Himalaya induced by surface hydrology.EarthandPlanetaryScienceLetters, 266(3-4): 332-344.

Chen L W, Lu Y Z, Guo R M, et al. 2001. Evolution of 3D tectonic stress field and fault movement in North China.ActaSeismologicaSinica(in Chinese), 23(4): 349-361.

Christensen N I, Mooney W D. 1995. Seismic velocity structure and composition of the continental crust: A global view.JournalofGeophysicalResearch:SolidEarth, 100(B6): 9761-9788.

Cui X F, Xie F R. 2001. The Space-time variations of present tectonic stress field in North China before and after 1976 Tangshan earthquake.EarthquakeResearchinChina(in Chinese), 17(3): 280-288.

Deng Q D, Zhang Y M, Xu G L, et al. 1979. On the tectonic stress field in China and its relation to plate movement.SeismologyandGeology(in Chinese), 1(1): 11-22.

Freeze R A, Cherry J A. 1977. Groundwater (in Chinese). Wu F J, trans. Beijing: Seismological Press.

Fu Y H, Miao J X, Zhang Z J, et al. 2009. Analysis on evolution of groundwater depression cones and its leading factors in North China Plain.ResourcesScience(in Chinese), 31(3): 394-399.

González P J, Tiampo K F, Palano M, et al. 2012. The 2011 Lorca earthquake slip distribution controlled by groundwater crustal unloading.NatureGeoscience, 5(11): 821-825.

Heki K. 2003. Snow load and seasonal variation of earthquake occurrence in Japan.EarthandPlanetaryScienceLetters, 207(1-4): 159-164.

Holzer T L. 1979. Elastic expansion of the lithosphere caused by groundwater depletion.JournalofGeophysicalResearch:SolidEarth(1978—2012), 84(B9): 4689-4698.

Huang Z X, Li H Y, Zheng Y J, et al. 2009. The lithosphere of North China Craton from surface wave tomography.EarthandPlanetaryScienceLetters, 288(1-2): 164-173.

Jia S X, Zhang X K. 2005. Crustal structure and comparison of different tectonic blocks in North China.ChineseJournalofGeophysics, 48(3): 672-683.

Konikow L F, Kendy E. 2005. Groundwater depletion: A global problem.HydrogeologyJournal, 13(1): 317-320.

Liu C, Shi Y L, Zheng L, et al. 2012. Relation between earthquake spatial distribution and tectonic stress accumulation in the North China Basin based on 3D visco-elastic modelling.ChineseJournalofGeophysics(in Chinese), 55(12): 3942-3957, doi: 10.6038/j.issn.0001-5733.2012.12.007.

Shah T, Roy A D, Qureshi A S, et al. 2003. Sustaining Asia′s groundwater boom: An overview of issues and evidence.NaturalResourcesForum, 27(2): 130-141.

Shi Y L, Cao J L. 2008. Effective viscosity of China continental lithosphere.EarthScienceFrontiers(in Chinese), 15(3): 82-95.

Sophocleous M. 2002. Interactions between groundwater and surface water: The state of the science.HydrogeologyJournal, 10(1): 52-67.

Stephen F, Héctor G. 2004. China: Towards Sustainable Groundwater Resource Use for Irrigated Agriculture on the North China Plain.∥Sustainable Groundwater Management, Lessons from Practice, Case Profile Collection. Washington, DC: World Bank.

Sun Y J, Dong S W, Fan T Y, et al. 2013. 3D rheological structure of the continental lithosphere beneath China and adjacent regions.ChineseJournalofGeophysics(in Chinese), 56(9): 2936-2946, doi: 10.3321/j.issn:0001-5733.2004.04.016.

Tian X F, Zelt C A, Wang F Y, et al. 2014. Crust structure of the North China Craton from a long-range seismic wide-angle-reflection/refraction data.Tectonophysics, 634: 237-245.

Wang D C, Zhang R Q, Shi Y H, et al. 2006. General Hydrogeology (in Chinese). Beijing: Geological Publishing House.

Xie F R, Cui X F, Zhao J T, et al. 2004. Regional division of the recent tectonic stress field in China and adjacent areas.ChineseJournalofGeophysics(in Chinese), 47(4): 654-662, doi: 10.3321/j.issn:0001-5733.2004.04.016.

Yu S E, Gong F H, Lu Z X, et al. 2000. The tectonic environment of strong earthquakes in North China Plain and its adjacent areas and their seismic tendency analyses.∥Crustal Structure and Stress Essays (in Chinese). Beijing: Institute of Crustal Stress, China Seismological Bureau, 16-30.

Zang S X, Wei R Q, Liu Y G. 2005. Three-dimensional rheological structure of the lithosphere in the Ordos block and its adjacent area.GeophysicalJournalInternational, 163(1): 339-356.

Zhang B, Zhang H, Shi Y L. 2015. Equivalent-body force approach on modeling elastic dislocation problem using finite element method.ChineseJournalofGeophysics(in Chinese), 58(5): 1666-1674, doi: 10.6038/cjg20150518.

Zhang Y, Li G M. 2013. Long-term evolution of cones of depression in shallow aquifers in the North China Plain.Water, 5(2): 677-697.

Zhang Z J, Chen Q F, Bai Z M, et al. 2011. Crustal structure and extensional deformation of thinned lithosphere in Northern China.Tectonophysics, 508(1-4): 62-72.

Zheng C M, Liu J, Cao G L, et al. 2010. Can China cope with its water crisis? -Perspectives from the North China Plain.GroundWater, 48(3): 350-354.

Zhu S B, Zhang P Z, Shi Y L. 2010. A study on the mechanisms of strong earthquake occurrence in the North China Basin.ChineseJournalofGeophysics(in Chinese), 53(6): 1409-1417, doi: 10.3969/j.issn.0001-5733.2010.06.019.

附中文參考文獻

陳連旺, 陸遠忠, 郭若眉等. 2001. 華北地區(qū)斷層運動與三維構造應力場的演化. 地震學報, 23(4): 349-361.

崔效鋒, 謝富仁. 2001. 1976年唐山地震前后華北地區(qū)現(xiàn)代構造應力場的時空變化特征. 中國地震, 17(3): 280-288.

鄧起東, 張裕明, 許桂林等. 1979. 中國構造應力場特征及其與板塊運動的關系. 地震地質(zhì), 1(1): 11-22.

費宇紅, 苗晉祥, 張兆吉等. 2009. 華北平原地下水降落漏斗演變及主導因素分析. 資源科學, 31(3): 394-399.

弗里澤R A, 徹里J A. 1987. 地下水. 吳靜方譯. 北京: 地震出版社.

柳暢, 石耀霖, 鄭亮等. 2012. 三維黏彈性數(shù)值模擬華北盆地地震空間分布與構造應力積累關系. 地球物理學報, 55(12): 3942-3957, doi: 10.6038/j.issn.0001-5733.2012.12.007.

石耀霖, 曹建玲. 2008. 中國大陸巖石圈等效黏滯系數(shù)的計算和討論. 地學前緣, 15(3): 82-95.

孫玉軍, 董樹文, 范桃園等. 2013. 中國大陸及鄰區(qū)巖石圈三維流變結構. 地球物理學報, 56(9): 2936-2946, doi: 10.3321/j.issn:0001-5733.2004.04.016.

王大純, 張人權, 史毅虹等. 2006. 水文地質(zhì)學基礎. 北京: 地質(zhì)出版社.

謝富仁, 崔效鋒, 趙建濤等. 2004. 中國大陸及鄰區(qū)現(xiàn)代構造應力場分區(qū). 地球物理學報, 47(4): 654-662, doi: 10.3321/j.issn:0001-5733.2004.04.016.

于慎諤, 龔復華, 盧造勛等. 2000. 華北平原及其鄰近地區(qū)強震構造環(huán)境與地震活動趨勢分析.∥地殼構造與地殼應力文集. 北京: 中國地震局地殼應力研究所, 16-30.

張貝, 張懷, 石耀霖. 2015. 有限元模擬彈性位錯的等效體力方法. 地球物理學報, 58(5): 1666-1674, doi: 10.6038/cjg20150518.

朱守彪, 張培震, 石耀霖. 2010. 華北盆地強震孕育的動力學機制研究. 地球物理學報, 53(6): 1409-1417, doi: 10.3969/j.issn.0001-5733.2010.06.019.

(本文編輯何燕)

Changes of crustal stress induced by groundwater over-pumping in North China Plain

PANG Ya-Jin, ZHANG Huai, CHENG Hui-Hong, SHI YAO-Lin*

KeyLaboratoryofComputationalGeodynamics,UniversityofChineseAcademyofSciences,Beijing100049,China

AbstractIt is well known that excessive mining of groundwater has induced surface subsidence in North China, especially in the past 50 years. While how groundwater mining unloads the lithosphere and causes stress changes is not fully understand. In this study, we set up a 2-D finite element model across a groundwater depression zone to quantitatively analyze deformation and crustal stress changes of North China Plain. The numerical results show that the maximum surface vertical uplift is 12.4 cm, extensional stress changes are up to 70 kPa and 35 kPa in upper and lower crust in the groundwater exploration area, respectively; while compressional deformation occurs at the edge of exploration zone with values up to 20 kPa. The annual tectonic stress change is 0.5 kPa. By comparing historical earthquakes (>M5.0) in North China Plain between two periods 1945—1979 and after 1980—2015, it is preliminary concluded the influence of groundwater unloading is significant, which may to some extent reduce the risk of great earthquakes in North China Plain.

KeywordsNorth China Plain; Groundwater mining; Stress field; Earthquakes; Numerical simulation

基金項目國家自然科學基金重大項目(41590865，41590864)，中國科學院國際合作創(chuàng)新團隊(KZZD-EW-TZ-19)資助.

作者簡介龐亞瑾，女,1990年生,博士,從事地球動力學數(shù)值模擬.E-mail：pang.yajin@163.com *通訊作者石耀霖,男,1944年生,教授,中國科學院院士,第三世界科學院院士,從事地球動力學研究. E-mail:shiyl@ucas.ac.cn

doi:10.6038/cjg20160420 中圖分類號P313

收稿日期2015-05-26，2015-12-24收修定稿

龐亞瑾， 張懷， 程惠紅等. 2016. 華北地區(qū)地下水開采對地殼應力的影響.地球物理學報，59(4):1394-1402,doi:10.6038/cjg20160420.

Pang Y J, Zhang H, Cheng H H, et al. 2016. Changes of crustal stress induced by groundwater over-pumping in North China Plain.ChineseJ.Geophys. (in Chinese),59(4):1394-1402,doi:10.6038/cjg20160420.