九寨溝MS 7.0地震對(duì)巴顏喀拉塊體東北緣活動(dòng)斷裂影響的有限元模擬

石富強(qiáng) 朱琳 王瑩 丁曉光 邱玉榮 石軍

陜西省地震局，西安市碑林區(qū)水文巷4號(hào) 710068

0 引言

巴顏喀拉塊體東北緣位于南北地震帶中段，區(qū)域內(nèi)分布有龍門山斷裂、東昆侖斷裂、西秦嶺北緣斷裂等多條大型活動(dòng)斷裂帶，為現(xiàn)今地殼變形顯著、地震活躍的區(qū)域。僅20世紀(jì)以來(lái)，該區(qū)域就發(fā)生了1933年疊溪7.5級(jí)、1976年松潘平武7.2級(jí)、2008年汶川8.0級(jí)、2017年九寨溝7.0級(jí)等多次7級(jí)以上強(qiáng)震。對(duì)該區(qū)域的地震危險(xiǎn)性跟蹤一直是地震科技工作者關(guān)注的焦點(diǎn)，M7專項(xiàng)工作組根據(jù)歷史地震活動(dòng)、大地測(cè)量及地球物理場(chǎng)數(shù)據(jù)對(duì)該區(qū)域給出了甘青川交界危險(xiǎn)區(qū)、西秦嶺中西段危險(xiǎn)區(qū)以及六盤山南-西秦嶺東危險(xiǎn)區(qū)等強(qiáng)震高風(fēng)險(xiǎn)區(qū)域。汶川 8.0級(jí)地震后，一些學(xué)者（Xiong et al，2010；Toda et al，2008；Wan et al，2010；單斌等，2009）基于歷史強(qiáng)震演化的斷層庫(kù)侖應(yīng)力變化的研究結(jié)果表明，東昆侖斷裂東段以及西秦嶺北緣斷裂西段的地震危險(xiǎn)性值得關(guān)注。此外，該區(qū)域還存在歷史地震破裂空段、低b值、小震活動(dòng)圖像、地震矩加速釋放等異常（陳為濤等，2013；M7專項(xiàng)工作組，2012）以及跨斷層水準(zhǔn)等地球物理場(chǎng)觀測(cè)數(shù)據(jù)異常①中國(guó)地震局第二監(jiān)測(cè)中心，2016，2017年度地震趨勢(shì)會(huì)商報(bào)告。。

庫(kù)侖應(yīng)力研究作為研究地震輪回中構(gòu)造應(yīng)力場(chǎng)時(shí)空變化的重要組成部分，已成為研究斷層應(yīng)力狀態(tài)及區(qū)域地震危險(xiǎn)性的重要手段之一。King等（1994）研究指出，1992年美國(guó)Landers 7.3級(jí)地震在Big Bear產(chǎn)生的0.2MPa左右的應(yīng)力加載量加速了Big Bear 6.5級(jí)地震的發(fā)生。1999年Izm it 7.4級(jí)地震發(fā)生前，Stein等（1997）等基于庫(kù)侖應(yīng)力計(jì)算指出該區(qū)域地震危險(xiǎn)性較高。而2004年P(guān)arkfield 6.0級(jí)地震的發(fā)生則受1983年Coalinga 6.5級(jí)地震及Nu?ez多個(gè)中等地震的應(yīng)力卸載的影響，出現(xiàn)了約11年的時(shí)間延遲（Bakun et al，1985）。此次九寨溝地震就發(fā)生在M7專項(xiàng)工作組（2012）給出的甘青川交界危險(xiǎn)區(qū)內(nèi)，基于庫(kù)侖應(yīng)力的強(qiáng)震作用研究表明，區(qū)域內(nèi)1976年松潘平武7.2級(jí)地震、2008年汶川8.0級(jí)地震等都對(duì)此次九寨溝地震發(fā)震斷層有較強(qiáng)的庫(kù)侖應(yīng)力加載作用（徐晶等，2017；單斌等，2017）。

有限元數(shù)值模擬技術(shù)作為探索構(gòu)造變形等的虛擬實(shí)驗(yàn)手段，在地震危險(xiǎn)性分析中得到了廣泛應(yīng)用（李玉江等，2013；Moreno et al，2010；Hergert et al，2010；陳連旺等，2008；祝愛(ài)玉等，2015）。為進(jìn)一步分析九寨溝地震對(duì)巴顏喀拉塊體東北緣主要活動(dòng)斷裂的加卸載效應(yīng)以及該區(qū)域未來(lái)可能的強(qiáng)震風(fēng)險(xiǎn)，本文構(gòu)建了巴顏喀拉塊體東北緣三維粘彈性有限元模型，模擬計(jì)算了九寨溝地震對(duì)周邊活動(dòng)斷裂的加卸載效應(yīng)，分析了震后地震活動(dòng)及其與九寨溝地震間的可能關(guān)系以及未來(lái)可能的強(qiáng)震風(fēng)險(xiǎn)。

1 模型構(gòu)建

1.1 幾何模型及物性參數(shù)

本文在詳細(xì)調(diào)研地震地質(zhì)資料的基礎(chǔ)上，通過(guò)對(duì)所調(diào)研的地質(zhì)資料進(jìn)行假設(shè)構(gòu)建了巴顏喀拉塊體東北緣有限元模型：模型縱向深100km，地殼分層、速度結(jié)構(gòu)以及密度結(jié)構(gòu)根據(jù)Crust1.0速度結(jié)構(gòu)給出（Laske et al，2013）；模型中包含禮縣-羅家堡斷裂、虎牙斷裂、雪山斷裂、東昆侖斷裂、龍門山斷裂、青川斷裂、青川-平武斷裂、臨潭-宕昌斷裂、白龍江斷裂、塔藏?cái)嗔?、龍日壩斷裂、岷江斷裂、西秦嶺北緣斷裂以及基于有限斷層反演給出的此次九寨溝地震發(fā)震斷裂等等；區(qū)域孕震層深度（段星北，1997；張國(guó)民等，2002），除九寨溝地震發(fā)震斷層外，假定其他斷層面下邊界深度為20km；除九寨溝地震發(fā)震斷層施加位錯(cuò)荷載外，其余斷層均為摩擦接觸斷層，斷層空間幾何結(jié)構(gòu)簡(jiǎn)化為由斷層地表出露沿著傾角向下延伸的空間斜面；下地殼及地幔粘滯性結(jié)構(gòu)根據(jù)邵志剛等（2008）、徐晶等（2017）的研究給出。地殼單元尺寸控制在10km內(nèi)，地幔單元尺寸控制在20km內(nèi)，斷層面上單元尺寸控制在3km內(nèi)，模型共包含359676個(gè)單元、67796個(gè)節(jié)點(diǎn)，構(gòu)建的有限元模型及對(duì)應(yīng)的物性參數(shù)分別見(jiàn)圖1和表1。

1.2 邊界約束及模擬方案

地震是構(gòu)造應(yīng)力作用下斷層應(yīng)力狀態(tài)達(dá)到一定程度而斷層破壞形成的自然災(zāi)害。在數(shù)值模擬中，由于我們并不知道斷層當(dāng)前絕對(duì)應(yīng)力狀態(tài)以及斷層自身的承載極限，因此，很難做到對(duì)地震自發(fā)的模擬。為此，本文首先對(duì)模型施加以GPS約束的邊界荷載來(lái)模擬初始構(gòu)造應(yīng)力場(chǎng)；然后，再施加地震同震位錯(cuò)，模擬計(jì)算給出九寨溝地震后區(qū)域的模擬構(gòu)造應(yīng)力場(chǎng)；最后，從震后模擬應(yīng)力場(chǎng)中扣除模擬初始應(yīng)力場(chǎng)，給出九寨溝地震的發(fā)生對(duì)區(qū)域主要活動(dòng)斷裂的加卸載效應(yīng)。其中，九寨溝地震位錯(cuò)模型參考波形反演結(jié)果（http：//www.cea-igp.ac.cn/tpxw275883.html）給出。所有模擬計(jì)算過(guò)程均在有限元開源代碼 Pylith（https：//geodynamics.org/cig/software/pylith/）平臺(tái)下進(jìn)行。

圖1 本文構(gòu)建的巴顏喀拉塊體東北緣有限元模型

表1 本文有限元模型力學(xué)參數(shù)

1.3 庫(kù)侖應(yīng)力計(jì)算

三維有限元模擬軟件并不能直接給出斷層面上的庫(kù)侖應(yīng)力，因此，需要通過(guò)計(jì)算給出的應(yīng)力張量來(lái)計(jì)算斷層面上的正應(yīng)力及剪應(yīng)力，進(jìn)而計(jì)算斷層面庫(kù)侖應(yīng)力變化。圖2為斷層面主方向（）與模型坐標(biāo)系（）相對(duì)關(guān)系，其中，分別為斷層走向（面向走向，右手在上盤）、傾向（逆沖為正）以及張拉方向（張為正）；分別為地理E、N向以及垂直地表向上的方向。α為斷層走向；β為斷層傾角。將斷層面應(yīng)力分量在模型坐標(biāo)系（、）下表示為

其中，T為矩陣轉(zhuǎn)置符號(hào)。

圖2 斷層變形模式示意圖（a）以及斷層面主方向與模型坐標(biāo)間的關(guān)系（b）

由圖2（b）可知，將模型坐標(biāo)系下的應(yīng)力張量繞z軸逆時(shí)針旋轉(zhuǎn)90°－α，然后，再繞主軸1逆時(shí)針旋轉(zhuǎn)β，便可得到主方向下的應(yīng)力張量。則

由式（1）、（2）給出的新的應(yīng)力分量σ33即為斷層面上的正應(yīng)力變化 Δσn，σ23、σ13分別為斷層面上沿傾向和走向的剪應(yīng)力分量。令斷層面優(yōu)勢(shì)破裂方向的滑動(dòng)角為θ，則斷層面優(yōu)勢(shì)破裂滑動(dòng)方向的剪應(yīng)力變化為

根據(jù)斷層面庫(kù)侖應(yīng)力變化的定義（King et al，1994），給出斷層面優(yōu)勢(shì)破裂方向的庫(kù)侖應(yīng)力變化為

其中，μ為斷層面有效摩擦系數(shù)。

2 模擬結(jié)果分析

2.1 模擬初始應(yīng)力場(chǎng)

以1999～2007年GPS觀測(cè)速度為約束，通過(guò)本文構(gòu)建的有限元?jiǎng)恿W(xué)模型計(jì)算了巴顏喀拉塊體東北緣的主應(yīng)力場(chǎng)（圖3）。模擬結(jié)果與基于小震震源機(jī)制反演的區(qū)域構(gòu)造應(yīng)力場(chǎng)（王曉山等，2015）以及基于GPS推算的區(qū)域應(yīng)變率場(chǎng)（陳長(zhǎng)云等，2013）基本一致。區(qū)域主要處于NW-SE拉張、NE-SW西向擠壓的動(dòng)力環(huán)境。九寨溝地震震中南側(cè)至龍門山斷裂，主張應(yīng)力逐漸減小，主壓方向逐漸偏轉(zhuǎn)并呈現(xiàn)出近EW向擠壓態(tài)勢(shì)。九寨溝地震震中附近主應(yīng)力方向與 USGS給出的震源機(jī)制解（https：//earthquake.usgs.gov/earthquakes/eventpage/us2000a5x1＃moment-tensor）的 P軸（方位角 104°）和 T軸（方位角 204°）基本一致。

2.2 九寨溝地震及其對(duì)周邊斷裂的加卸載效應(yīng)

以同震位錯(cuò)為荷載模擬給出的九寨溝地震震中附近的同震位移場(chǎng)如圖4所示。由圖4可見(jiàn)，震中最大同震位移約為9cm，九寨溝地震的發(fā)生對(duì)東昆侖斷裂、塔藏?cái)嗔?、西秦嶺北緣斷裂東段以及龍門山斷裂均為左旋作用，同震位移場(chǎng)的左旋作用有利于左旋活動(dòng)為主的東昆侖斷裂、塔藏?cái)嗔岩约拔髑貛X北緣斷裂東段斷層面的應(yīng)力積累；而對(duì)于具有右旋逆沖活動(dòng)的龍門山斷裂則為卸載作用。同理，對(duì)西秦嶺北緣斷裂西段應(yīng)為卸載影響。為進(jìn)一步定量分析九寨溝地震對(duì)各斷裂的加卸載效應(yīng)，我們基于本文模擬結(jié)果，根據(jù)式（2）、（4）計(jì)算了九寨溝地震在這些斷裂斷層面上引起的庫(kù)侖應(yīng)力變化（圖5）。

圖3 本文模擬給出的巴顏喀拉塊體東北緣初始主應(yīng)力場(chǎng)

圖4 數(shù)值模擬給出的九寨溝地震震中附近同震位移場(chǎng)

由圖5可見(jiàn)，九寨溝地震對(duì)塔藏?cái)嗔盐鞫?、虎牙斷裂加載效應(yīng)顯著，達(dá)7056Pa；東昆侖斷裂東段、龍日壩斷裂、迭部-白龍江斷裂西端和東端、臨潭-宕昌斷裂東段以及西秦嶺北緣斷裂東段也均受到不同程度的加載作用；塔藏?cái)嗔褨|段以及岷江斷裂受九寨溝地震的影響，斷層面庫(kù)侖應(yīng)力變化表現(xiàn)為卸載，最大卸載量達(dá)73984Pa。位于此次九寨溝地震破裂面正前方的雪山斷裂斷層面庫(kù)侖應(yīng)力變化較為復(fù)雜：由圖4模擬同震位移場(chǎng)可知，九寨溝地震對(duì)雪山斷裂西段的影響表現(xiàn)為擠壓，對(duì)東段的影響表現(xiàn)為拉張，根據(jù)庫(kù)侖應(yīng)力計(jì)算公式可知，當(dāng)斷層面應(yīng)力變化主要表現(xiàn)為垂直于斷層面的擠壓/拉張時(shí)，有效摩擦系數(shù)的取值對(duì)庫(kù)侖應(yīng)力變化影響顯著，甚至可以影響其極性變化。Wan等（2010）在分析汶川地震對(duì)周邊斷裂影響時(shí)也給出同樣的結(jié)論。從本文模擬結(jié)果看，當(dāng)有效摩擦系數(shù)較?。é?0.1）時(shí)，九寨溝地震對(duì)雪山斷裂西段表現(xiàn)為庫(kù)侖應(yīng)力加載，對(duì)東段表現(xiàn)為卸載；當(dāng)有效摩擦系數(shù)較大（μ=0.7）時(shí)，對(duì)西段表現(xiàn)為卸載作用，而對(duì)東段表現(xiàn)為加載（圖5）。

此外，此次九寨溝地震還對(duì)汶川地震破裂區(qū)東北端的青川斷裂以及青川-平武斷裂的西南段有一定的庫(kù)侖應(yīng)力加載，最大加載量可達(dá)1139Pa。邵志剛等（2010）研究汶川地震的震后影響發(fā)現(xiàn)，青川斷裂庫(kù)侖應(yīng)力加載非常顯著，實(shí)際上青川地區(qū)也是汶川地震余震活動(dòng)非?；钴S的區(qū)域（石軍，2008），但從2013年開始，該區(qū)域地震活動(dòng)水平明顯減弱。此次九寨溝地震后不到2個(gè)月，9月30日青川（32.27°N，105°E）再次發(fā)生5.4級(jí)地震。分析認(rèn)為，此次青川地震可能是由九寨溝地震的觸發(fā)作用所致。

圖5 有限元模擬給出的不同有效摩擦系數(shù)下斷層面庫(kù)侖應(yīng)力變化

3 結(jié)論與討論

基于構(gòu)建的三維有限元模型，以同震位錯(cuò)為荷載模擬了九寨溝地震的發(fā)生及其對(duì)巴顏喀拉塊體東北緣主要活動(dòng)斷裂的加卸載效應(yīng)。數(shù)值模擬結(jié)果顯示，九寨溝地震的發(fā)生對(duì)處于甘青川交界危險(xiǎn)區(qū)內(nèi)的東昆侖斷裂東段至塔藏?cái)嗔盐鞫螏?kù)侖應(yīng)力加載效應(yīng)顯著，最大可達(dá)7056Pa；對(duì)六盤山南-西秦嶺東危險(xiǎn)區(qū)的西秦嶺北緣斷裂東段也有一定程度的加載效應(yīng)，約為346Pa。此外，計(jì)算結(jié)果還顯示，九寨溝地震對(duì)龍日壩斷裂、虎牙斷裂、青川-平武斷裂西段、迭部-白龍江斷裂西段和東段、臨潭-宕昌斷裂東段等也均有不同程度的庫(kù)侖應(yīng)力加載效應(yīng)；對(duì)岷江斷裂、塔藏?cái)嗔褨|段庫(kù)侖應(yīng)力卸載效應(yīng)顯著，最大卸載量可達(dá)73984Pa。

地震是區(qū)域構(gòu)造應(yīng)力積累到一定程度且超越斷層面承載極限而引起斷層剪切破壞的動(dòng)力過(guò)程。但實(shí)際上我們很難獲取斷層面的絕對(duì)應(yīng)力狀態(tài)和斷層面的承載極限。庫(kù)侖應(yīng)力作為由地震發(fā)生而引起的構(gòu)造應(yīng)力場(chǎng)變化在斷層面優(yōu)勢(shì)破裂方向的體現(xiàn)，在地震危險(xiǎn)性定性分析判斷中有著廣泛的應(yīng)用。汶川地震后，許多學(xué)者通過(guò)庫(kù)侖應(yīng)力計(jì)算分析了汶川地震對(duì)區(qū)域活動(dòng)斷層的加卸載效應(yīng)及地震危險(xiǎn)性（Toda et al，2008；單斌等，2009；邵志剛等，2010；李玉江等，2013）。隨后的2013年蘆山7.0級(jí)地震、岷縣漳縣6.6級(jí)地震以及此次九寨溝7.0級(jí)地震均發(fā)生于文獻(xiàn)報(bào)道的庫(kù)侖應(yīng)力加載區(qū)域。這些強(qiáng)地震的時(shí)空叢集活動(dòng)也表明，巴顏喀拉地塊東段在未來(lái)一段時(shí)間內(nèi)可能繼續(xù)保持這種活躍態(tài)勢(shì)。從庫(kù)侖應(yīng)力變化的角度來(lái)看，此次九寨溝地震對(duì)甘青川交界危險(xiǎn)區(qū)內(nèi)的主要活動(dòng)斷裂（東昆侖斷裂東段和塔藏?cái)嗔眩┮约拔挥诹P山南-西秦嶺東危險(xiǎn)區(qū)的主要活動(dòng)斷裂（西秦嶺北緣斷裂東段），均具有不同程度的庫(kù)侖應(yīng)力加載效應(yīng)。