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

朱家正， 孫玉軍

中國地質科學院地球深部探測中心, 北京 100037

0 引言

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

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

1 研究背景

1.1 構造地質背景

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

圖1 大崗山水庫區域構造背景(a)及庫區地形圖(b)(a) 中橙色實線代表二級地塊邊界,藍色實線代表一級地塊邊界,黑色實線代表主要斷裂,黃色方框是研究區域所在位置,黃色五角星代表水庫大壩位置; (b) 中藍色填充區域為水庫主要蓄水區域,黃色方塊代表大壩所在位置.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年建成后,大致分兩階段使得蓄水水位達到預定水位.于2014年3月開啟第一階段蓄水,水位達到1000 m；于2014年11月開啟第二階段蓄水,于2015年5月達到預定水位,水庫水位由2013年蓄水前的960 m水位上升到1130 m的正常蓄水位.庫區的地震數據表明(圖2),水庫蓄水后研究區域內地震發生頻率隨著水位的變化發生了明顯的增加,特別是伴隨著水位的兩次大幅度變化,庫區周邊地震頻度以及震級均發生增大.本研究將101.9°E—102.4°E, 29.2°N—29.8°N作為研究范圍,以國家地震臺網的地震數據為依據,繪出蓄水前后三年內庫區地震平面分布圖(圖3).可以看到,在水庫蓄水之前,地震主要分布在庫區西側鮮水河斷裂,多為ML3級以下微震,在鮮水河斷裂西側呈條帶狀分布；2015年5月水庫經過第二階段蓄水達到正常水位后,該區域的地震分布仍然主要分布在庫區西北的鮮水河斷裂磨西段,但叢集性更強,且出現ML3～4級以上地震8次,ML4～5級地震3次,活動性變強.但在庫區東部的龍門山斷裂南段,在蓄水達到正常水位后幾乎不再出現微震,庫區南側的大涼山斷裂附近的地震活動性也有一定的降低.該現象是否與水庫蓄水有直接關系值得深入研究.

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

圖3 水庫蓄水前后研究區域地震活動情況(a) 2011年11月—2014年11月庫區地震分布; (b) 2014年11月—2017年11月庫區地震分布. 黃色五角星標識為水庫所在位置.a, b, c, d分別代表龍門山斷裂、鮮水河斷裂、大涼山斷裂以及錦屏山—小金河斷裂.○代表震中位置,綠色為ML0～2級,黃色為ML2～3級,橙色代表ML3～4級,紅色代表ML4級以上地震.紅色橢圓所圈出區域為蓄水后地震活動性增強區域,藍色橢圓所圈出的區域為蓄水后地震活動性降低區域.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.

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

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

2 研究方法

2.1 孔隙彈性體基本方程

孔隙彈性體的一般本構關系為

(1)

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

(2)

(3)

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

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

幾何方程：

(4)

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

(5)

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

2.2 建立模型

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

2.3 邊界條件和計算參數

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

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

(6)

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

2.4 庫侖應力計算

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

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

(7)

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

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

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

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

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

3 計算結果及分析

3.1 不同模型參數對結果影響

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

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

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

圖5e、5f為保持斷層的擴散系數不變,而改變非斷層區域的擴散系數.結果顯示非斷層區域的擴散系數對結果影響較為明顯,在第二階段蓄水之后,當擴散系數c=1 m2·s-1時,參考點處的庫侖應力和孔隙壓力值都快速增長,而擴散系數低于該值的另外兩個模型顯示,庫侖應力和孔隙壓力變化較為平緩.

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

圖5 不同模型計算得到的庫侖應力變化量以及孔隙壓力隨時間變化(a、b) 模型1結果; (c、d) 模型2結果; (e、f) 模型3結果. (g、h) 模型4結果; (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.

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

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

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

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

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

圖6 地震活動性增強區域與減弱區域選定點的庫侖應力和孔隙壓力變化(a) 地震活動性增強區域點1結果; (b) 地震活動性減弱區域點2結果.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 對磨西斷裂影響

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

圖7 2017年1月模型計算得到的深度為5 km (a, b),10 km (c, d),15 km (e, f)的庫侖應力變化量以及孔隙壓力分布(單位：Pa)○ 代表第二階段蓄水之后距離剖面垂直距離范圍2 km內的地震震中位置在平面上的投影,灰色代表震級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.

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

圖8 模型4計算得到的2017年1月不同剖面孔隙壓力(a，c)和庫侖應力變化量(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)

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

3.4 對大渡河斷裂影響

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

圖9 模型4計算得到的2017年1月深度為5 km(a, b)和10 km (c, d)的庫侖應力變化量分布 (a, c)和(b, d)分別為考慮水滲流和不考慮水滲流得到了庫侖應力變化量分布(單位：Pa). ○代表第二階段蓄水之后距離剖面垂直距離范圍2 km內的地震震中位置在平面上的投影,灰色為ML2級以下,黃色為ML2～3級,橙色代表ML3～4級,紅色代表ML4級以上地震.藍色三角形代表石棉地震震中投影.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 討論與結論

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

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

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

本文的主要結論如下：

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

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

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