四川大崗山水庫蓄水對地震活動影響的數(shù)值模擬研究

朱家正， 孫玉軍

中國地質(zhì)科學(xué)院地球深部探測中心, 北京 100037

0 引言

川西地區(qū)位于我國南北地震帶的南段,是構(gòu)造形態(tài)復(fù)雜、地殼運(yùn)動劇烈、地震活動頻繁的地區(qū).在深大斷裂的控制下,塊體之間的水平運(yùn)動和垂直差異運(yùn)動都表現(xiàn)得十分明顯,而強(qiáng)烈地震大多數(shù)就發(fā)生在這些深大斷裂帶上(李碧雄等, 2014).同時,我國西南地區(qū)水力資源豐富,近些年在該地區(qū)修建了多個大型水電站,如向家壩、溪洛渡、白鶴灘、烏東德等大型水電站,大型水電工程建設(shè)所面臨的水庫地震問題是重要的科學(xué)問題.一方面我國西南地區(qū)是構(gòu)造活動十分強(qiáng)烈的區(qū)域,構(gòu)造地震活躍；另一方面,修建水庫后,大型水電工程蓄水造成的附加載荷和滲流作用一般都會改變地區(qū)原有的應(yīng)力狀態(tài).水庫地震一般震級較小,但是其造成的附加應(yīng)力改變可能與該區(qū)域初始的構(gòu)造應(yīng)力狀態(tài)疊加,從而造成大的構(gòu)造地震發(fā)生.如2008年發(fā)生在龍門山斷裂帶上的汶川地震,毫無疑問,如此大震級的地震發(fā)生是構(gòu)造運(yùn)動的直接結(jié)果(Chen, 2009),但是后續(xù)的研究表明建立在龍門山斷裂帶上的紫坪鋪水庫對斷層的庫侖應(yīng)力影響為正值(Klose, 2012;周斌等, 2010;孫玉軍等, 2012;雷興林等, 2008;張貝和石耀霖, 2010),即對斷層滑動有一定的促發(fā)作用,可能促使斷層滑動提前發(fā)生(Ge et al., 2009).水庫蓄水可以誘發(fā)地震是一個已經(jīng)被證實(shí)的科學(xué)問題(Gough and Gough, 1970;Gupta, 2002;Gupta and Chadha, 1995;Talwani, 1997;夏其發(fā), 1993;孫玉軍等, 2012).我國最著名的水庫地震是1962年發(fā)生在廣東的新豐江水庫地震,震級達(dá)到了6.1級(丁原章等, 1983;程惠紅等, 2012).因此在構(gòu)造活動強(qiáng)烈的區(qū)域修建大型水庫等水電建設(shè)工程,其對該區(qū)域構(gòu)造活動特別是大型斷裂的影響,是未來值得深入考慮和探究的科學(xué)問題.該問題不僅涉及到大型水電工程的安全穩(wěn)定性問題,同時對建設(shè)區(qū)域的防震減災(zāi)也有重要的參考價值.

四川大崗山水電站修建于2013年,屬于大渡河上游修建的一系列水電站之一.由于龍門山斷裂、鮮水河斷裂、安寧河斷裂等一系列大型斷裂在此區(qū)域交匯,地質(zhì)背景復(fù)雜,構(gòu)造運(yùn)動活躍,屬于地震高烈度區(qū).同時,大崗山水庫屬于高山峽谷型水庫,壩高達(dá)210 m,對于這類特殊性水庫,其蓄水造成的附加災(zāi)害對該區(qū)重要的大型斷裂影響需要重點(diǎn)關(guān)注(程萬正, 2013).

1 研究背景

1.1 構(gòu)造地質(zhì)背景

大崗山水庫庫區(qū)位于青藏高原東南緣川滇南北向構(gòu)造帶北端,川滇地塊、華南地塊、巴顏喀拉地塊的邊緣交界處(圖1),是地殼運(yùn)動劇烈、構(gòu)造形態(tài)復(fù)雜、地震活動十分頻繁的地區(qū)(張培震等, 2013；趙靜等, 2018).龍門山斷裂帶、鮮水河斷裂帶以及安寧河斷裂帶在此區(qū)域組成Y型構(gòu)造帶.大涼山斷裂、小金河斷裂、玉農(nóng)希斷裂等主要斷裂也在這一區(qū)域附近交匯,其中鮮水河斷裂南段的磨西斷裂以及大渡河斷裂都直接通過水庫蓄水區(qū)域(阮祥等, 2017).龍門山斷裂帶是一條長約500 km、寬約30～50 km、沿NE-SW方向展布的巨大斷裂帶,按照由西向東的順序,龍門山斷裂帶主要包含龍門山后山斷裂(茂縣—汶川斷裂)、中央斷裂(映秀—北川斷裂)和山前斷裂(安縣—灌縣斷裂,亦稱彭縣—灌縣斷裂、江油—灌縣斷裂).這些斷裂(北西段)都以逆沖滑動為主,兼具一定的右旋走滑分量(Xu et al., 2008;鄧起東等, 1994;陳國光等, 2007).從地震剖面上看,逆沖面呈鏟狀幾何形態(tài),上地殼向西陡傾,中地殼向西逐漸傾斜(Xu et al., 2008 ).張培震等的研究結(jié)果顯示,龍門山斷裂西南段存在(4±2.5)mm·a-1的縮短速率,以及(7.5±2)mm·a-1左旋走滑速率(張培震等, 2003).趙靜等采用DEFNODE負(fù)位錯反演程序估算了龍門山斷裂西南端的閉鎖程度與變形狀態(tài),并綜合GPS反演結(jié)果和跨斷層水準(zhǔn)結(jié)果,認(rèn)為目前龍門山斷裂帶西南端在大部分段落處于強(qiáng)閉鎖狀態(tài),且依然有發(fā)生大地震的可能性(趙靜等, 2018).鮮水河斷裂帶北起甘孜東谷附近,向南經(jīng)過爐霍、道孚、康定一線,至石棉縣安順場一帶逐漸減弱消失(郭長寶等, 2015).晚新生代以來,鮮水河斷裂表現(xiàn)出強(qiáng)烈的左旋走滑運(yùn)動,是松潘—甘孜造山帶內(nèi)部一條大型走滑斷裂,橫切了松潘—甘孜造山帶的主體,系造山運(yùn)動后期陸內(nèi)變形的產(chǎn)物,晚新生代以來的位移總規(guī)模在60 km左右(許志琴等, 2007).以乾寧惠遠(yuǎn)寺一帶為界,鮮水河斷裂帶北西段和南東段的活動速率有所差異,北西段活動速率為10～15 mm·a-1(聞學(xué)澤等, 1989),南東段小于10 mm·a-1,一般為 5 mm·a-1左右(郭長寶等, 2015).李鐵明等通過重力、GPS資料計算結(jié)果,得出鮮水河斷裂磨西段先進(jìn)整體呈現(xiàn)左旋走滑,速率為4.41 mm·a-1(李鐵明等, 2019).

圖1 大崗山水庫區(qū)域構(gòu)造背景(a)及庫區(qū)地形圖(b)(a) 中橙色實(shí)線代表二級地塊邊界,藍(lán)色實(shí)線代表一級地塊邊界,黑色實(shí)線代表主要斷裂,黃色方框是研究區(qū)域所在位置,黃色五角星代表水庫大壩位置; (b) 中藍(lán)色填充區(qū)域為水庫主要蓄水區(qū)域,黃色方塊代表大壩所在位置.Fig.1 Regional tectonic background (a) and topographic map (b) of Dagangshan ReservoirThe solid orange lines in the figure (a) represent the boundaries of the secondary blocks. The solid blue lines represent the boundaries of the primary blocks. The solid black lines represent the main faults. The yellow box represents the location of the study area, and the yellow pentacle represents the location of the reservoir dam. The blue filled area in the figure (b) is the main water storage area of the reservoir, and the yellow square represents the location of the dam.

1.2 水庫蓄水前后地震活動背景

大崗山水庫在2013年建成后,大致分兩階段使得蓄水水位達(dá)到預(yù)定水位.于2014年3月開啟第一階段蓄水,水位達(dá)到1000 m；于2014年11月開啟第二階段蓄水,于2015年5月達(dá)到預(yù)定水位,水庫水位由2013年蓄水前的960 m水位上升到1130 m的正常蓄水位.庫區(qū)的地震數(shù)據(jù)表明(圖2),水庫蓄水后研究區(qū)域內(nèi)地震發(fā)生頻率隨著水位的變化發(fā)生了明顯的增加,特別是伴隨著水位的兩次大幅度變化,庫區(qū)周邊地震頻度以及震級均發(fā)生增大.本研究將101.9°E—102.4°E, 29.2°N—29.8°N作為研究范圍,以國家地震臺網(wǎng)的地震數(shù)據(jù)為依據(jù),繪出蓄水前后三年內(nèi)庫區(qū)地震平面分布圖(圖3).可以看到,在水庫蓄水之前,地震主要分布在庫區(qū)西側(cè)鮮水河斷裂,多為ML3級以下微震,在鮮水河斷裂西側(cè)呈條帶狀分布；2015年5月水庫經(jīng)過第二階段蓄水達(dá)到正常水位后,該區(qū)域的地震分布仍然主要分布在庫區(qū)西北的鮮水河斷裂磨西段,但叢集性更強(qiáng),且出現(xiàn)ML3～4級以上地震8次,ML4～5級地震3次,活動性變強(qiáng).但在庫區(qū)東部的龍門山斷裂南段,在蓄水達(dá)到正常水位后幾乎不再出現(xiàn)微震,庫區(qū)南側(cè)的大涼山斷裂附近的地震活動性也有一定的降低.該現(xiàn)象是否與水庫蓄水有直接關(guān)系值得深入研究.

圖2 研究區(qū)域地震活動性與水庫水位關(guān)系Fig.2 The relationship between seismic activity and reservoir water level in the study area

圖3 水庫蓄水前后研究區(qū)域地震活動情況(a) 2011年11月—2014年11月庫區(qū)地震分布; (b) 2014年11月—2017年11月庫區(qū)地震分布. 黃色五角星標(biāo)識為水庫所在位置.a, b, c, d分別代表龍門山斷裂、鮮水河斷裂、大涼山斷裂以及錦屏山—小金河斷裂.○代表震中位置,綠色為ML0～2級,黃色為ML2～3級,橙色代表ML3～4級,紅色代表ML4級以上地震.紅色橢圓所圈出區(qū)域為蓄水后地震活動性增強(qiáng)區(qū)域,藍(lán)色橢圓所圈出的區(qū)域為蓄水后地震活動性降低區(qū)域.Fig.3 Seismic activity in the study area before and after the reservoir storage(a) The distribution of earthquakes in the reservoir area from November 2011 to November 2014; (b) The distribution of earthquakes in the reservoir area from November 2014 to November 2017. The yellow star indicates the location of the dam. a, b, c, d represent Longmenshan fault, Xianshuihe fault, Daliangshan fault and Jinpingshan-Xiaojinhe fault respectively. ○ represents the location of the epicenter. Green represents magnitude ML0～2, yellow represents magnitude ML2～3, orange represents magnitude ML3～4, and red represents above magnitude ML4. The red ellipse circles the area of increased seismicity after water storage, and the blue ellipse circles the area of decreased seismicity after water storage.

為了定量分析大崗山水庫蓄水對該區(qū)域斷層的影響,本文依據(jù)構(gòu)造地質(zhì)特征及精細(xì)的DEM數(shù)據(jù)建立了三維孔隙彈性有限元數(shù)值模型(圖4),通過計算大崗山水庫蓄水過程造成庫區(qū)斷層的庫侖應(yīng)力變化,定量分析了大崗山水庫蓄水對地震活動性的影響.

圖4 三維孔隙彈性有限元模型(模型中不同顏色代表不同材料參數(shù))Fig.4 Three-dimensional poroelastic finite element model (Different colors in the model represent different material parameters)

2 研究方法

2.1 孔隙彈性體基本方程

孔隙彈性體的一般本構(gòu)關(guān)系為

(1)

水的滲流擴(kuò)散方程(暫不考慮飽和孔隙水不排水效應(yīng)的影響)為(Bell and Nur,1978;Roeloffs, 1988)

(2)

(3)

其中c為水力學(xué)擴(kuò)散系數(shù),單位m2·s-1,k是介質(zhì)的滲透率,單位m2,η是水的流體動力學(xué)黏滯性系數(shù),單位Pa·s,泊松比、剪切模量、Skempton系數(shù)對結(jié)果的影響主要是通過改變擴(kuò)散系數(shù)，本文將c作為計算參數(shù).c的取值范圍主要參考Talwani等(2007).

在小形變的幾何方程和平衡方程的約束下,施加合理的邊界條件就可以進(jìn)行有限元求解.

幾何方程：

(4)

平衡方程(忽略體力項)：

(5)

本文采用孫玉軍等開發(fā)的三維孔隙彈性有限元模擬程序進(jìn)行計算(孫玉軍等, 2012).流體滲流和孔隙彈性耦合的研究方法不僅可以用于研究水庫地震，同時在研究頁巖氣開采過程中壓裂誘發(fā)地震時也被廣泛應(yīng)用(Hemam et.al., 2021; Wong et.al., 2021; Yang et.al., 2021).

2.2 建立模型

三維模型的建立考慮地層的垂向分層,模型表面層采用實(shí)際的地表地形數(shù)據(jù),本研究采用SRTM30 m數(shù)字高程數(shù)據(jù)(https:∥www.usgs.gov/products/maps/topo-maps).設(shè)定模型的深度為25 km,模型范圍為101.9°E—102.4°E,29.2°N—29.8°N,采用三棱柱單元進(jìn)行網(wǎng)格剖分,并對庫區(qū)部分進(jìn)行網(wǎng)格加密,模型總結(jié)點(diǎn)數(shù)702259,總單元數(shù)1399412.數(shù)值模型可以考慮實(shí)際的分層性,但是由于缺乏實(shí)際滲流參數(shù),本研究只區(qū)分了斷層和非斷層區(qū)域,對其設(shè)置了不同的模型參數(shù).水平層狀網(wǎng)格中,對蓄水區(qū)域進(jìn)行加密,分辨率為60 m,非蓄水區(qū)域網(wǎng)格大小為900 m.垂向上網(wǎng)格分為15層,深度為25 km,底部四層網(wǎng)格無斷層穿過,斷層產(chǎn)狀傾向北西西,傾角83°.

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

在研究區(qū)域范圍內(nèi),對水庫蓄水位以上的其他區(qū)域,給定力的邊界和孔隙壓力邊界均為0,對水位以下地區(qū)給定力的邊界和孔隙壓力邊界條件為

p(x,y,z,t)=ρωgh(x,y,z,t),

(6)

ρω為水的密度,g=9.8 m·s-2為重力加速度,h(x,y,z,t)為t時刻水庫水位減去不同地點(diǎn)的地形高程值,即實(shí)際水深,p(x,y,z,t)為對應(yīng)的水庫蓄水造成的庫底水壓力,對于力的邊界,模型側(cè)邊界和底邊界均為滑動邊界條件,對于孔隙壓力,側(cè)面給定為等滲流梯度,底面固定孔隙壓力為0.

2.4 庫侖應(yīng)力計算

庫侖應(yīng)力的變化量定義為(Harris,1998;石耀霖和曹建玲, 2010)

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

(7)

Δτ為所考慮的斷層面上剪應(yīng)力的變化量,μ是摩擦系數(shù),Δσn為斷層面上的正應(yīng)力改變量(拉伸為正),Δp為孔隙壓力改變量.本文以coulomb-s代表庫侖應(yīng)力變化量.

同時,為了比較不同參數(shù)對計算結(jié)果的影響,分別建立不同的模型進(jìn)行比較(表1).

表1 各個模型中計算參數(shù)的選取Table 1 Selection of calculation parameters in each model

鮮水河斷裂南東段的磨西斷裂主要走向為160°(李大虎等, 2015),在靠近庫區(qū)的位置取其走向為166°,根據(jù)程佳、徐錫偉等對鮮水河南端水平和垂直滑動速率的總結(jié),將斷層滑動角參數(shù)設(shè)置為21°(徐錫偉等, 2003;程佳等, 2012).龍門山斷裂南段大渡河斷裂的斷層參數(shù)參考前人的區(qū)域應(yīng)力場及震源機(jī)制解研究來選取(易桂喜等, 2013;李君等, 2019;杜瑤等, 2016;聞學(xué)澤等, 2008).計算庫侖應(yīng)力所采用的接收斷層參數(shù)選取如表2所示.

表2 計算庫侖應(yīng)力采用的各斷層參數(shù)Table 2 The parameters of each fault used to calculate the change of coulomb stress (coulomb-s)

3 計算結(jié)果及分析

3.1 不同模型參數(shù)對結(jié)果影響

為討論不同模型參數(shù)對計算結(jié)果影響,在水庫蓄水后庫區(qū)西北側(cè)選取參考點(diǎn)1坐標(biāo)為(102.12°,29.53°,-5 km),分別對比三組模型在該點(diǎn)處的孔隙壓力以及庫侖應(yīng)力變化,庫侖應(yīng)力的計算采用鮮水河斷裂磨西段的斷層參數(shù)(表2).圖5中coulomb-s代表庫侖應(yīng)力值.

圖5a、5b為模型1計算得到結(jié)果.模型1中通過保持其他參數(shù)不變而只改變了斷層的彈性模量以探究其對結(jié)果影響.從圖中可以看出,在2013年12月水位第一次上漲時水庫蓄水就對該參考點(diǎn)位置造成正的庫侖應(yīng)力影響,但該蓄水過程中孔隙壓力幾乎沒有變化,說明此時水庫蓄水的彈性載荷效應(yīng)在起主要作用,隨著水庫在2014年12月水位快速增長110 m,三種不同模型參數(shù)對應(yīng)的庫侖應(yīng)力變化量也增加明顯,同時參考點(diǎn)的庫侖應(yīng)力隨斷層彈性模量的增大而減小,但由于三個模型并未改變滲流狀態(tài),因此三種模型所對應(yīng)的孔隙壓力基本一樣(圖5b),隨著后期孔隙壓力的增加,水庫蓄水過程造成了正的庫侖應(yīng)力增加.

圖5c、5d為模型2計算結(jié)果.模型2只改變了斷層的擴(kuò)散系數(shù)而保持其他系數(shù)不變.可以從結(jié)果中看出孔隙壓力的結(jié)果在2015年4月出現(xiàn)明顯改變；而庫侖應(yīng)力的三條結(jié)果曲線在2014年12月水位大幅度改變之前基本保持一致,計算結(jié)果顯示此時水的滲流還未達(dá)到該深度.在2014年12月之后庫侖應(yīng)力隨著孔隙壓力的增大而逐漸變大,同時斷層的擴(kuò)散系數(shù)越大,在參考點(diǎn)處造成的庫侖應(yīng)力變化量則越小.

圖5e、5f為保持?jǐn)鄬拥臄U(kuò)散系數(shù)不變,而改變非斷層區(qū)域的擴(kuò)散系數(shù).結(jié)果顯示非斷層區(qū)域的擴(kuò)散系數(shù)對結(jié)果影響較為明顯,在第二階段蓄水之后,當(dāng)擴(kuò)散系數(shù)c=1 m2·s-1時,參考點(diǎn)處的庫侖應(yīng)力和孔隙壓力值都快速增長,而擴(kuò)散系數(shù)低于該值的另外兩個模型顯示,庫侖應(yīng)力和孔隙壓力變化較為平緩.

圖5g、5h為改變非斷層區(qū)域的排水泊松比所得到的計算結(jié)果.計算結(jié)果顯示，在保持?jǐn)U散系數(shù)不變的情況下，改變非斷層區(qū)域的排水泊松比對孔隙壓力的結(jié)果沒有影響.水庫蓄水初期庫侖應(yīng)力隨著排水泊松比的增大而減小，隨著時間的增加三種計算結(jié)果在一點(diǎn)交會后，兩者關(guān)系轉(zhuǎn)為正相關(guān).該結(jié)果可以通過孔隙彈性的本構(gòu)關(guān)系來理解，也即單純的彈性介質(zhì)下，庫侖應(yīng)力隨排水泊松比增大而減小，但是根據(jù)孔隙彈性的本構(gòu)關(guān)系，孔隙壓力的增加，高排水泊松比更有利于庫侖應(yīng)力的增大.

圖5 不同模型計算得到的庫侖應(yīng)力變化量以及孔隙壓力隨時間變化(a、b) 模型1結(jié)果; (c、d) 模型2結(jié)果; (e、f) 模型3結(jié)果. (g、h) 模型4結(jié)果; (b、h)中紅色、綠色和粉色線重合.Fig.5 The change of Coulomb stress and pore pressure calculated by different models vary with time(a,b) The results of Model 1; (c,d) The results of Model 2; (e,f) The results of Model 3. The red, green and pink lines in (b) coincide.

以上幾種模型中都可以看到在水庫水位變化引起庫侖應(yīng)力變化主要受到兩種效應(yīng)影響,一種是蓄水造成彈性載荷變化引起的庫侖應(yīng)力變化,另一種是水的滲流造成孔隙壓力變化而引起的庫侖應(yīng)力變化.但這兩種效應(yīng)在時間上的表現(xiàn)有所不同,蓄水造成彈性載荷變化所引起的庫侖應(yīng)力變化是瞬時的,如第一階段蓄水之后庫侖應(yīng)力都隨之快速增加,此時斷層彈性模量對計算結(jié)果有較大影響,但此時孔隙壓力并未明顯改變.而第二種效應(yīng)中孔隙壓力改變所造成的庫侖應(yīng)力變化相對于水庫蓄水過程有所滯后,如第二階段蓄水之后的庫侖應(yīng)力增長都隨著孔隙壓力的變化而變化,說明水庫達(dá)到設(shè)計水位后,主要是水的滲流作用影響庫侖應(yīng)力,而對此影響較大的參數(shù)是斷層或非斷層區(qū)域的擴(kuò)散系數(shù),擴(kuò)散系數(shù)越大,孔隙壓力越大,從而對參考點(diǎn)造成的庫侖應(yīng)力改變越大.

從孔隙彈性介質(zhì)材料來看,水庫蓄水過程所產(chǎn)生的這兩種效應(yīng),可能會互相疊加增強(qiáng),也有可能會互相抵消減弱.為了更詳細(xì)的分析兩種效應(yīng)的作用,我們在計算庫侖應(yīng)力的時候可以分別對其進(jìn)行考慮,即考慮水庫滲流造成的庫侖應(yīng)力改變(彈性和滲流的綜合效應(yīng))和不考慮水庫滲流造成的庫侖應(yīng)力變化(彈性效應(yīng)).

3.2 水庫蓄水對不同區(qū)域地震活動性影響

綜合上述三種模型,以及前人的研究成果(Talwani et al., 2007;陳建業(yè)等, 2011)選取模型5為最終的模型計算參數(shù).在庫區(qū)西側(cè)的磨西斷裂選取參考點(diǎn)1(102.12°E,29.53°N,-5 km),大渡河斷裂選取參考點(diǎn)2(102.20°E,29.57°N,-5 km)(圖3a),對水庫蓄水后的地震活動性增強(qiáng)和減弱區(qū)域進(jìn)行庫侖應(yīng)力和孔隙壓力進(jìn)行分析.對點(diǎn)1庫侖應(yīng)力計算選取磨西斷裂的斷層參數(shù),對點(diǎn)2的庫侖應(yīng)力計算采用大渡河斷裂的斷層參數(shù)(表2).圖中coulomb-e為不考慮孔隙壓力影響計算得到的庫侖應(yīng)力變化量.

對比兩圖發(fā)現(xiàn),圖6a中在水庫蓄水后不含孔隙壓力的庫侖應(yīng)力coulomb-e(彈性效應(yīng))為正值,并且隨水位增加而迅速增長,使斷層更加危險；而圖6b的結(jié)果與之相反,不含孔隙壓力的庫侖應(yīng)力coulomb-e(彈性效應(yīng))為負(fù)值,隨水位增長而快速下降,對斷層起到抑制作用.這是由于水庫蓄水區(qū)域跨越了龍門山和鮮水河兩大斷裂的南端,位置分別處于磨西斷裂的下盤,以及大渡河斷裂的上盤.因此庫水對斷層施加的彈性載荷使得磨西斷裂的主要區(qū)域地震活動性增加,使得大渡河斷裂的主要控制區(qū)域地震活動性減弱.但是,隨著水位保持高位,水不斷向斷層中滲流,孔隙壓力隨之增大,使得庫區(qū)兩側(cè)的庫侖應(yīng)力coulomb-s都轉(zhuǎn)變?yōu)檎?斷層趨向于更加危險.

圖6 地震活動性增強(qiáng)區(qū)域與減弱區(qū)域選定點(diǎn)的庫侖應(yīng)力和孔隙壓力變化(a) 地震活動性增強(qiáng)區(qū)域點(diǎn)1結(jié)果; (b) 地震活動性減弱區(qū)域點(diǎn)2結(jié)果.Fig.6 Changes in Coulomb stress and pore pressure at selected points in the areas of enhanced and weakened seismic activity(a) The result of point 1 in the seismic activity enhancement area; (b) The result of point 2 in the seismic activity weakened area.

3.3 對磨西斷裂影響

將磨西斷裂參數(shù)作為計算庫侖應(yīng)力的參數(shù),得到2017年1月的庫侖應(yīng)力變化量與孔隙壓力的平面分布(圖7).孔隙壓力增大的區(qū)域?qū)?yīng)水庫主蓄水區(qū)域,且隨著深度的增加,孔隙壓力值逐漸減小,如在深度5 km孔隙壓力最大值為63 kPa,而在10 km和15 km孔隙壓力最大值分別為22 kPa和10 kPa；同時由于斷層的擴(kuò)散系數(shù)一般較大,孔隙壓力在斷層區(qū)域一般較高.庫侖應(yīng)力變化量正值區(qū)域除了在水庫庫區(qū)范圍外,主要分布在庫區(qū)的北西和南東兩個方向,同時隨深度增加逐漸減小,在5 km深度處庫侖應(yīng)力變化量正極大值可達(dá)45 kPa,在10 km和15 km庫侖應(yīng)力變化量正極大值分別為17 kPa和7 kPa.庫侖應(yīng)力變化量的負(fù)值區(qū)域主要分布在庫區(qū)東北側(cè)和西南側(cè),在5 km深度負(fù)的極大值可達(dá)-6 kPa,在10 km和15 km負(fù)的極大值分別為-4 kPa和-2 kPa.深度為10 km以上時,水的滲流造成的孔隙壓力較小,庫侖應(yīng)力變化量的區(qū)域性分布更為明顯,庫區(qū)西南方向出現(xiàn)了庫侖應(yīng)力變化量的負(fù)值區(qū)域,在水庫蓄水后的地震增加區(qū)域(庫區(qū)北西側(cè))看到明顯的2 kPa左右的正值區(qū)域,說明水庫蓄水對該區(qū)域地震的發(fā)生具有促進(jìn)作用.

圖7 2017年1月模型計算得到的深度為5 km (a, b),10 km (c, d),15 km (e, f)的庫侖應(yīng)力變化量以及孔隙壓力分布(單位：Pa)○ 代表第二階段蓄水之后距離剖面垂直距離范圍2 km內(nèi)的地震震中位置在平面上的投影,灰色代表震級ML2以下微震,黃色代表ML2～3級地震,橙色代表ML3～4級地震,紅色代表ML4級以上地震.Fig.7 The coulomb stress and pore pressure distribution at the depth of 5 km (a, b),10 km (c, d), and 15 km (e, f) in January, 2017 (unit: Pa)○ Represents the projection of the earthquake epicenter position on the plane within 2 kilometers of the vertical distance from the profile after the second stage of water storage. Gray represents seismic magnitude below ML2, yellow represents magnitude ML2～3, orange represents magnitude ML3～4, and red represents magnitude above ML4.

以參考點(diǎn)(102.12°,29.53°)為中心位置對計算結(jié)果做縱向剖面(圖8),孔隙壓力主要分布在水庫區(qū)域且有隨著斷層在地下的延展更快速地擴(kuò)散(圖8c),庫侖應(yīng)力變化量的分布與圖7對應(yīng),在庫區(qū)南北兩側(cè)形成負(fù)值區(qū)域,中心位置為正值；庫區(qū)東側(cè)以及南側(cè)的5～10 km深處存在庫侖應(yīng)力變化量的負(fù)值區(qū)域(圖8b、8d),西側(cè)沿斷層的延展分布庫侖應(yīng)力變化量正值區(qū)域.

圖8 模型4計算得到的2017年1月不同剖面孔隙壓力(a，c)和庫侖應(yīng)力變化量(b, d)分布圖剖面A-A′為102.15°E的縱切面,剖面B-B′為29.53°N的橫切面. 剖面位置已在圖3a中顯示.Fig.8 The pore pressure (a, c) and Coulomb stress (b, d) distribution of different profiles calculated by Model 4 in January 2017 (The positions of the two profiles are shown in Fig.3a)

從應(yīng)力觸發(fā)區(qū)域和地震活動的對應(yīng)關(guān)系上,可以看到水庫蓄水對10 km和15 km深度范圍的地震在空間位置上具有明顯的對應(yīng)關(guān)系：在10 km深度處的應(yīng)力影響區(qū)邊緣出現(xiàn)了較為集中的ML2～3級地震分布,15 km深度范圍內(nèi)的ML3～5級的地震也集中分布在西北側(cè)應(yīng)力觸發(fā)區(qū)域的邊緣,說明水庫蓄水對該區(qū)域地震發(fā)生具有促進(jìn)作用,而不同深部的應(yīng)力影區(qū)域的地震活動均受到抑制.

3.4 對大渡河斷裂影響

以大渡河斷裂參數(shù)(表2)作為計算依據(jù)得出2017年1月庫侖應(yīng)力變化量(圖9a, 9c)及不考慮水滲流的庫侖應(yīng)力變化量(圖9b, 9d)平面分布圖.不考慮水滲流的庫侖應(yīng)力變化量在庫區(qū)東側(cè)存在一個明顯的負(fù)值區(qū)域,在5 km深度靠近水庫大壩位置的庫侖應(yīng)力變化量負(fù)極大值達(dá)到-12 kPa,考慮水滲流的情況下該值增加至-8 kPa左右,且與采用磨西斷裂參數(shù)時一樣，在蓄水區(qū)域出現(xiàn)正值區(qū),5 km深度處庫侖應(yīng)力變化量正值極大值達(dá)到50 kPa,10 km深度處的庫侖應(yīng)力變化量分布與采用磨西斷裂得到的結(jié)果相比,庫侖應(yīng)力變化量(考慮水滲流)負(fù)的極值增大1 kPa左右,達(dá)到-5.6 kPa.對比采用磨西斷裂參數(shù)計算得到的結(jié)果可以看出,盡管采用的兩條斷層參數(shù)有所改變,但水庫蓄水所造成的庫侖應(yīng)力變化量正、負(fù)區(qū)域(即應(yīng)力觸發(fā)和應(yīng)力影區(qū)域)并未明顯改變,除了造成庫區(qū)的庫侖應(yīng)力變化量為正值外,在庫區(qū)北西和南東為庫侖應(yīng)力變化量正值區(qū)(應(yīng)力觸發(fā)區(qū)域),庫區(qū)北西側(cè)相對較大；而在庫區(qū)北東和南西則為庫侖應(yīng)力變化量負(fù)值區(qū)(應(yīng)力影區(qū)域),庫區(qū)北東負(fù)值更明顯.該影響與水庫蓄水后地震活動性比較一致.

圖9 模型4計算得到的2017年1月深度為5 km(a, b)和10 km (c, d)的庫侖應(yīng)力變化量分布 (a, c)和(b, d)分別為考慮水滲流和不考慮水滲流得到了庫侖應(yīng)力變化量分布(單位：Pa). ○代表第二階段蓄水之后距離剖面垂直距離范圍2 km內(nèi)的地震震中位置在平面上的投影,灰色為ML2級以下,黃色為ML2～3級,橙色代表ML3～4級,紅色代表ML4級以上地震.藍(lán)色三角形代表石棉地震震中投影.Fig.9 The distribution of the change of coulomb stress at the depth of 5 km (a, b) and 10 km (c, d) in January, 2017 for the model 4 (unit: Pa)(a, c) and (b, d) are the change of Coulomb stress with and without considering pore-pressure, respectively. ○ Represents the projection of the earthquake epicenter position on the plane within 2 kilometers of the vertical distance from the profile after the second stage of water storage. Gray represents microseisms below magnitude ML2, yellow represents magnitude ML2～3 earthquakes, orange represents magnitude ML3～4 earthquakes, and red represents earthquakes above magnitude ML4. The blue triangle represents the asbestos earthquake epicenter projection.

4 討論與結(jié)論

本文采用三維孔隙彈性模型,研究大崗山水庫蓄水前后水位變化對庫區(qū)地震及斷層的影響.根據(jù)地震臺網(wǎng)數(shù)據(jù)，已知經(jīng)過水位變化的第二階段蓄水之后,水庫附近的微震活動變化明顯,首先,庫區(qū)西北的鮮水河斷裂南段附近地震震級增大,出現(xiàn)了ML3級以上的地震,且數(shù)量增多,叢集性更明顯,由于蓄水位置處于磨西斷裂下盤,促進(jìn)了斷層活動性,根據(jù)計算結(jié)果,在地震叢集增強(qiáng)的區(qū)域靠近庫區(qū)的磨西斷裂一側(cè)為庫侖應(yīng)力變化量正值區(qū)(應(yīng)力觸發(fā)區(qū))；而在地震活動性降低的庫區(qū)北東區(qū)域,庫侖應(yīng)力變化量為負(fù)值(應(yīng)力影區(qū)).但是隨著滲流的持續(xù),庫侖應(yīng)力變化量負(fù)值區(qū)域也會逐漸變?yōu)檎?地震活動性增強(qiáng).根據(jù)前人研究成果,鮮水河斷裂帶的最大剪切應(yīng)變速率為40～60×10-9/a(Wang and Shen, 2020).如果取最大剪切應(yīng)變速率為50×10-9/a,假設(shè)彈性模量為80 GPa的情況下,應(yīng)力累積速率為4 kPa/a,而蓄水造成的庫侖應(yīng)力變化量增加在淺層5 km深度最大可達(dá)到45 kPa(圖10),這相當(dāng)于11年的構(gòu)造運(yùn)動產(chǎn)生的應(yīng)力累積.

地震增強(qiáng)區(qū)域和減弱區(qū)域的平面方向與臨近的鮮水河斷裂和龍門山斷裂存在一定的夾角,這是否對應(yīng)了兩條斷裂帶地下的展布狀態(tài),需要進(jìn)一步討論.由于水庫位于磨西斷裂的下盤且主要蓄水區(qū)域都被其直接穿過,導(dǎo)致無論水的滲流還是庫水的彈性載荷都促使其更加活躍；而大渡河斷裂的上盤由于受到水庫的垂向重力載荷,導(dǎo)致活動性減弱,地震減少.但是隨著水的不斷滲流,以及水位的穩(wěn)定,水庫彈性載荷效應(yīng)造成的對地震的抑制作用會不斷被削弱,庫區(qū)的庫侖應(yīng)力可能都處于增加狀態(tài),即整個區(qū)域的地震活動性更強(qiáng).根據(jù)地震數(shù)據(jù),該區(qū)域在蓄水之后三年的時間中都沒有再發(fā)生地震,說明現(xiàn)實(shí)中庫水的彈性載荷對該區(qū)域地震的抑制效應(yīng)比預(yù)計的更強(qiáng)；另外一種可能是水庫區(qū)域地殼結(jié)構(gòu)的黏彈性效應(yīng)，黏彈性效應(yīng)往往也造成庫侖應(yīng)力分布的不均勻，也可能對該區(qū)域造成負(fù)的庫侖應(yīng)力，從而造成滲流效應(yīng)在該區(qū)域發(fā)生作用需要更長的時間.因此，在未來的研究中對于地殼結(jié)構(gòu)的黏彈性效應(yīng)也應(yīng)該有所考慮.水庫北東側(cè)的斷裂活動受到抑制，產(chǎn)生閉鎖效應(yīng)，造成龍門山斷裂西南段的應(yīng)變累積速率加快，水位的快速下降有可能造成此前積累的應(yīng)變能快速釋放，增加大震發(fā)生的可能性.

由于模型和實(shí)際斷層詳細(xì)數(shù)據(jù)的限制,沒有給出鮮水河和龍門山斷裂的鏟狀結(jié)構(gòu),實(shí)際情況中斷層的傾角可能會隨著深度增加而減小,因此水的滲流影響會在深部隨著斷層擴(kuò)散至更遠(yuǎn)的區(qū)域,造成西北部區(qū)域的孔隙壓力增加值比模型計算得到的值更大,該區(qū)域的斷層可能更加活躍.另外，在鮮水河斷裂的西北區(qū)域也可能存在與斷裂平行的小斷裂，如果這些小斷裂與鮮水河主斷裂參數(shù)一致，那么目前計算得到的庫侖應(yīng)力也適用于這些小斷裂，從而促進(jìn)斷裂活動.因此未來如果獲取了該區(qū)域更為詳細(xì)的活動斷裂探測數(shù)據(jù)，則可以利用該模型開展更為詳細(xì)的計算分析.

本文的主要結(jié)論如下：

(1)從孔隙彈性材料的角度看,水庫蓄水所造成的庫侖應(yīng)力改變主要有兩種效應(yīng),一種是蓄水造成彈性載荷變化引起的庫侖應(yīng)力變化,另一種是水的滲流造成孔隙壓力變化而引起的庫侖應(yīng)力變化.這兩種效應(yīng)在時間上的表現(xiàn)有所不同,蓄水造成彈性載荷變化所引起的庫侖應(yīng)力變化是瞬時的；而孔隙壓力改變所造成的庫侖應(yīng)力變化相對于水庫蓄水過程有所滯后.

(2)依據(jù)斷層實(shí)際滑動參數(shù)計算得到的庫侖應(yīng)力結(jié)果來看,大崗山水庫蓄水過程所造成的庫侖應(yīng)力改變使得鮮水河斷裂磨西段庫侖應(yīng)力增加,而使得大渡河斷裂庫侖應(yīng)力降低.這與水庫蓄水后地震臺網(wǎng)觀測得到的磨西斷裂西側(cè)地震活動性明顯增強(qiáng),而大渡河斷裂東側(cè)地震活動性減弱的事實(shí)一致.

致謝感謝兩位匿名審稿人和編輯提出的建設(shè)性意見.