基于分裂節(jié)點(diǎn)法的地震同震和震后形變數(shù)值模擬及其在汶川大地震中的應(yīng)用

孫云強(qiáng)，羅 綱，黃祿淵

1.福建農(nóng)林大學(xué)交通與土木工程學(xué)院，福建 福州 350002；2.武漢大學(xué)測繪學(xué)院，湖北 武漢 430079；3.應(yīng)急管理部國家自然災(zāi)害防治研究院，北京 100085

0 引言

自1910年彈性回跳理論提出以來(Reid, 1910)，人們逐漸認(rèn)識到地震與斷層之間的密切相關(guān)性。地震位錯理論得到迅速的發(fā)展。1958年Steketee最早將彈性位錯理論引入地震學(xué)，并推導(dǎo)出彈性半無限空間中點(diǎn)源位錯的地表位移格林函數(shù)(Steketee, 1958)。隨后，該理論被廣泛應(yīng)用于地震同震、震后應(yīng)力應(yīng)變分析中(Maruyama, 1964; Mansinha and Smylie, 1971; 陳運(yùn)泰等，1979)。Okada(1985,1992)在已有研究基礎(chǔ)上，推導(dǎo)出均勻彈性半無限空間中三種位錯源(走滑位錯、傾滑位錯、引張位錯)引起的同震位移場及應(yīng)力場的解析表達(dá)式，由于求解速度快，Okada模型已成為計(jì)算彈性半空間同震效應(yīng)的經(jīng)典方法。但是均勻彈性半無限空間模型不能夠很好地處理地球介質(zhì)的橫向與縱向不均勻，在實(shí)際應(yīng)用中仍具有較大局限性(邵志剛等，2008; 張貝等，2015; 黃祿淵等，2019)。另外，學(xué)者也逐漸認(rèn)識到震后形變，即大地震發(fā)生后，地球介質(zhì)的黏彈性松弛效應(yīng)。為了使模型更符合真實(shí)的地球情況，同時考慮大地震同震及震后效應(yīng)，研究者相繼又開發(fā)了一系列的模型(Pollitz，1992;Sun，1992; Sabadini and Vermeersen, 1997; Wang et al., 2003; 黃祿淵等，2019)。Pollitz(1992)基于自由震蕩簡正模方法，計(jì)算了無重力情況下，黏彈性地球模型的地震位移場與應(yīng)變場。Sun(1992)提出了球?qū)ΨQ分層介質(zhì)的位錯理論。Wang et al.(2003)提出利用矩陣傳播算法計(jì)算地震應(yīng)力場格林函數(shù)，開發(fā)了黏彈性的垂向分層模型下的地震同震和震后形變計(jì)算程序PSGRN/PSCMP等。上述解析/半解析的方法計(jì)算速度快，但是在處理地球介質(zhì)的橫向不均勻性、復(fù)雜的斷層面幾何形狀時仍然遇到了困難。隨著計(jì)算機(jī)技術(shù)的發(fā)展，有限元數(shù)值模擬得到快速發(fā)展。有限元方法能夠根據(jù)已有的觀測資料建立更符合實(shí)際情況的地質(zhì)模型，在地震同震及震后形變數(shù)值模擬研究中有著更廣泛的應(yīng)用(邵志剛等，2008; 李玉江等，2009; Zhu and Zhang, 2013; Xu et al., 2020)。

地震的發(fā)生是斷層上應(yīng)力釋放和調(diào)整的過程，會導(dǎo)致區(qū)域應(yīng)力場發(fā)生變化，進(jìn)而影響區(qū)域的地震活動性。由于地球介質(zhì)的復(fù)雜性和地球深部的不可入性(陳運(yùn)泰, 2009)，獲取地殼的絕對應(yīng)力值仍較為困難。但是基于彈性位錯理論，可以計(jì)算每次地震產(chǎn)生的應(yīng)力張量變化及庫侖應(yīng)力變化。大地震導(dǎo)致的區(qū)域最優(yōu)破裂面及斷層面上的庫侖應(yīng)力變化(同震庫侖應(yīng)力變化及震后的黏彈性松弛過程導(dǎo)致的庫侖應(yīng)力變化)作為一個重要的指標(biāo)，在分析不同地震之間的相互關(guān)系、解釋余震序列分布及區(qū)域地震危險(xiǎn)性評估中都發(fā)揮了重要的作用(King et al., 1994; Stein et al., 1997; Harris and Simpson, 1998;Freed, 2005;萬永革等，2007;胡才博和蔡永恩, 2016; 汪建軍和許才軍, 2017; 朱曉杰和何建坤, 2019; Xu et al., 2020)。例如，King et al.(1994)通過計(jì)算1992年M7.4 Landers地震造成的同震庫侖應(yīng)力變化，認(rèn)為Landers地震觸發(fā)了隨后發(fā)生的M6.5 Big Bear地震，且Landers地震的余震序列大部分發(fā)生在庫侖應(yīng)力變化為正的區(qū)域。Stein et al.(1997)通過計(jì)算庫侖應(yīng)力變化成功解釋了土耳其1939—1992年的地震序列。萬永革等(2007)計(jì)算了青藏高原東北部的庫侖應(yīng)力演化，發(fā)現(xiàn)后發(fā)大地震大部分都發(fā)生在先發(fā)大地震導(dǎo)致的庫侖應(yīng)力為正的區(qū)域，觸發(fā)率達(dá)85%。Freed and Lin(2001)分析了1992年Landers地震的震后庫侖應(yīng)力變化，發(fā)現(xiàn)由于中下地殼上地幔的黏彈性松弛造成的庫侖應(yīng)力變化在1999年M7.1 Hector Mine地震震中附近增加了約0.1～0.2 MPa，從而觸發(fā)了該地震的發(fā)生。朱曉杰和何建坤(2019)通過計(jì)算1970年通海MS7.7地震的同震及震后庫侖應(yīng)力變化，認(rèn)為該地震使得小江斷裂與紅河斷裂上的庫侖應(yīng)力增加，潛在地震活動性增強(qiáng)。

文中首先基于Abaqus的二次開發(fā)平臺，開發(fā)了三維Maxwell(馬克斯威爾體)黏彈性有限元程序，并引入分裂節(jié)點(diǎn)法來模擬地震位錯。其次，通過計(jì)算走滑斷層的地震同震及震后效應(yīng)，并與解析/半解析解對比，驗(yàn)證了程序的正確性與可靠性。最后，將程序應(yīng)用于青藏高原東緣，計(jì)算了2008年MW7.9汶川地震導(dǎo)致的同震應(yīng)力變化及震后中下地殼上地幔的黏彈性應(yīng)力調(diào)整，分析了汶川地震對2013年MW6.6蘆山地震及2017年MW6.5九寨溝地震的影響。

1 計(jì)算方法

1.1 分裂節(jié)點(diǎn)法

地震同震位移可以處理為斷層兩側(cè)大小相等、方向相反的兩部分位移(在斷層兩側(cè)同時施加地震同震位移一半的位移約束，且方向相反)。在有限元數(shù)值模擬中，研究者采用了多種不同方式來處理斷層，包括采用一定寬度的弱化帶(Luo and Liu, 2010)、摩擦接觸(李玉江等，2013;Zhu and Zhang, 2013)、分裂節(jié)點(diǎn)法(Melosh and Raefsky, 1981; 邵志剛等，2008; 黃祿淵等，2018)等。其中分裂節(jié)點(diǎn)法能夠簡單有效地將斷層位移以初始位移的形式引入有限元數(shù)值模擬中，是一種簡單有效且物理意義明確的方法(Melosh and Raefsky, 1981)。因此文中擬使用分裂節(jié)點(diǎn)法來模擬斷層地震。首先以一個包含兩個單元的一維有限元模型為例，介紹分裂節(jié)點(diǎn)法的具體處理過程(圖1)。

箭頭表示地震斷層錯動方向；表示地震斷層的同震位移；數(shù)字表示單元節(jié)點(diǎn)號

假設(shè)模型遵循靜力平衡方程。在有限元分析中，平衡方程可以表示為:

KU=F

(1)

其中，K為單元剛度矩陣；U為位移；F為等效節(jié)點(diǎn)荷載。

用節(jié)點(diǎn)位移表示的單元a的節(jié)點(diǎn)平衡方程為：

(2)

將地震引起的位移量移動到方程式的右側(cè)：

(3)

同理，可以得到單元b的節(jié)點(diǎn)平衡方程：

(4)

將單元剛度矩陣組裝到全局，形成總體剛度矩陣，得到整個模型的節(jié)點(diǎn)平衡方程組：

(5)

通過上式，可以看到：在不改變整體剛度矩陣的情況下，分裂節(jié)點(diǎn)法將地震位移以初始位移邊界的形式加載到右端項(xiàng)，可以有效處理地震不連續(xù)面。

1.2 三維黏彈性本構(gòu)關(guān)系

地震同震的持續(xù)時間較短，同震變形可以視為地球介質(zhì)的彈性響應(yīng)。而地震發(fā)生后，形變隨著時間還將繼續(xù)積累。研究指出這是由于不同機(jī)制，如震后余滑、黏彈性松弛等共同作用的影響(邵志剛等，2008)，其中黏彈性松弛現(xiàn)象被廣泛應(yīng)用于地震的震后形變研究(Pollitz, 1992; Wang et al., 2001;Hu et al., 2004; 邵志剛等，2008; Shao et al., 2011)。因此，作為簡化，文中只考慮黏彈性松弛的影響，采用Maxwell黏彈性材料來模擬中下地殼上地幔物質(zhì)。模型采用增量法求解每一時間步的應(yīng)變增量。應(yīng)變增量由黏性應(yīng)變增量和彈性應(yīng)變增量組成：

{dε}={dεv}+{dεe}

(6)

其中{dεv}和{dεe}分別表示黏性和彈性的應(yīng)變增量；{ }為張量形式。

三維黏彈性本構(gòu)關(guān)系描述為：

{dεv}=Q-1{σt}dt=Q-1({σt-dt}+

{dσ})dt{dεe}=D-1{dσ}

(7)

其中，{σt}為t時刻的應(yīng)力張量；dt為時間增量；{dσ}為應(yīng)力張量增量；D為彈性材料矩陣；Q為與黏度相關(guān)的材料矩陣。

(8)

(9)

其中，E為楊氏模量;υ為泊松比;η為黏度。

2 概念性模型

根據(jù)上述理論，文中開發(fā)了模擬地震同震及震后效應(yīng)的三維黏彈性有限元程序。首先以一個概念性模型(包含一條直立的走滑斷層)為例(圖2)，對程序的可靠性進(jìn)行驗(yàn)證。模型尺寸設(shè)置為100 km×100 km×100 km。在垂向方向上模型介質(zhì)分為兩層：上層20 km為彈性層，模擬上地殼物質(zhì)；楊氏模量和泊松比分別設(shè)置為8.75×1010Pa和0.25。往下80 km為黏彈性層，模擬中下地殼上地幔；楊氏模量和泊松比分別為1.1×1011Pa和0.25，黏度設(shè)置為1×1020Pa·s(算例1)。模型的斷層長度為20 km，斷層延伸深度為20 km；通過在斷層兩側(cè)各施加0.5 m的位移來模擬同震位移為1 m的右旋走滑地震(其中斷層左側(cè)施加0.5 m的同震位移，斷層右側(cè)施加-0.5 m的同震位移)。在模型的四個側(cè)面邊界施加滑移邊界條件(法向固定，切線方向可以自由移動)，底部邊界固定。

模型計(jì)算的地表同震位移分布結(jié)果顯示：同震位移分布表現(xiàn)出明顯的關(guān)于斷層成對稱(圖3a)或反對稱(圖3b，3c)分布的花瓣圖，與已有計(jì)算得到的結(jié)果也基本一致(邵志剛等，2008; 黃祿淵等，2018)。為了更進(jìn)一步地驗(yàn)證程序的正確性，便于與解析/半解析解進(jìn)行比較，文中選取跨斷層的一個剖面AA′(圖2)，并畫出剖面節(jié)點(diǎn)上與斷層走向平行的同震位移Uy(圖4a藍(lán)色線段)。同時文中利用Okada模型(Okada, 1985)計(jì)算了平行斷層走向方向的同震位移Uy(圖4a橙色虛線)。從圖中可以看到有限元模型計(jì)算的結(jié)果與解析解得到的結(jié)果是非常吻合的(圖4a)。這也驗(yàn)證了文中所編寫程序的正確性。

白色箭頭表示模型的斷層為右旋走滑；AA′為剖面位置

圖中白色線段表示斷層位置；Ux，Uy，Uz分別表示x方向、y方向和z方向的同震位移分布

另外，為了探討地球介質(zhì)橫向不均勻性對地震同震位移的影響，文中在算例1的基礎(chǔ)上又設(shè)置了兩組不同的算例(表1)：算例2將位于斷層左側(cè)的上地殼楊氏模量設(shè)置為1.75×1011Pa；算例3將位于斷層左側(cè)的上地殼楊氏模量設(shè)置為4.375×1010Pa；其他參數(shù)保持不變。同樣取剖面AA′上平行斷層走向的同震位移的結(jié)果，如圖4b所示。可以看到，斷層兩側(cè)楊氏模量的不同，會改變沿著斷層走向方向的同震位移關(guān)于斷層成反對稱的趨勢。斷層左側(cè)上地殼的楊氏模量增大，位于斷層左側(cè)上地殼的同震位移減小，而斷層右側(cè)上地殼的同震位移增大(圖4b紫色線)；斷層左側(cè)上地殼的楊氏模量減小，位于斷層左側(cè)上地殼的同震位移增大，而斷層右側(cè)上地殼的同震位移減小(圖4b黃色線)。可見地球介質(zhì)的橫向不均勻?qū)Φ卣鹜鹞灰朴酗@著的影響，在實(shí)際的應(yīng)用中不容忽視。

同時，文中分析了中下地殼上地幔的黏彈性松弛對上地殼形變的影響，以算例1為計(jì)算模型，計(jì)算剖面AA′的同震及震后200年累積的形變Uy(圖4c，4d藍(lán)色線)，可以看到地震后上地殼的位移仍繼續(xù)積累(圖4c，4d)。為了進(jìn)一步分析不同中下地殼上地幔黏度對震后效應(yīng)的影響，在算例1的基礎(chǔ)上，文中還分別設(shè)置了中下地殼上地幔黏度為1×1018Pa·s(算例4)、1×1019Pa·s(算例5)、1×1021Pa·s(算例6)的模型，并計(jì)算不同算例的同震及震后形變(圖4c，4d，表1)。通過對比不同的算例結(jié)果，可以看到，中下地殼上地幔的黏度對震后形變有顯著的影響。并且在震后相同的時間段里，中下地殼上地幔的黏度越小，上地殼的震后形變持續(xù)積累越大；而中下地殼上地幔的黏度越大，上地殼的震后形變持續(xù)積累越小。這個結(jié)果(圖4)與相關(guān)學(xué)者用PSGRN/PSCMP程序計(jì)算得到的結(jié)果(徐昊等，2018)也比較一致。

a—地震同震位移分布(實(shí)線表示有限元模型的結(jié)果，虛線為Okada模型計(jì)算得到的結(jié)果)；b—不同算例的地震同震位移；c—震后200年剖面同震位移變化圖；d—同震和震后200年的位移總變化量

表1 單斷層有限元模型的部分材料參數(shù)

在驗(yàn)證了程序可靠性的基礎(chǔ)上，接下來文中將該程序應(yīng)用于青藏高原東緣實(shí)例中；計(jì)算了2008年MW7.9汶川地震的同震及震后庫侖應(yīng)力變化，分析汶川地震對蘆山地震及九寨溝地震的影響。

3 三維黏彈性有限元程序在青藏高原東緣的應(yīng)用

3.1 青藏高原東緣地震活動

青藏高原東緣的川西地區(qū)是中國重要的地震活動帶，區(qū)域構(gòu)造變形強(qiáng)烈，活動斷層展布，控制著一系列歷史強(qiáng)震的發(fā)生(徐錫偉等，2005; 張培震等，2008; 鄭文俊等，2020)。2008年5月12日在這個區(qū)域的龍門山斷裂上發(fā)生了汶川MW7.9大地震，造成了巨大的人員傷亡和經(jīng)濟(jì)損失。隨后在2013年4月20日發(fā)生了MW6.6蘆山地震；2017年8月8日在九寨溝地區(qū)又發(fā)生了MW6.5九寨溝地震等(圖5)。如此頻繁的地震也引起了社會和學(xué)界極大的關(guān)注。研究者在這個區(qū)域開展了大量的研究，包括地震破裂過程、地震復(fù)發(fā)周期等，積累了豐富的資料(王衛(wèi)民等，2008; 張培震等，2008; Wang et al., 2011)?；谶@些觀測的數(shù)據(jù)資料建立模型，并計(jì)算地震導(dǎo)致的庫侖應(yīng)力變化，在分析不同地震之間的相互關(guān)系、解釋余震序列分布及區(qū)域未來地震危險(xiǎn)性等方面發(fā)揮了重要的作用(解朝娣等，2010; 賈科和周仕勇, 2018; Luo and Liu, 2018; 董培育等，2019; 黃祿淵等，2019)。許多學(xué)者基于不同模型計(jì)算了汶川地震、蘆山地震與九寨溝地震導(dǎo)致的應(yīng)力變化及幾次地震之間的相互關(guān)系(Parsonset et al., 2008; Toda et al., 2008; 單斌等，2009; Luo and Liu, 2018; 黃祿淵等，2019; Xu et al., 2020)。Parsons et al.(2008)采用地震波反演得到的均勻彈性有限斷層模型，計(jì)算了汶川地震導(dǎo)致的庫侖應(yīng)力變化，結(jié)果顯示汶川地震引起的同震庫侖應(yīng)力變化在龍門山斷裂帶南段增加了0.1 MPa，而該區(qū)域在2013年發(fā)生了MW6.6蘆山地震。單斌等(2009，2017)利用彈性位錯理論計(jì)算，結(jié)果顯示汶川地震造成蘆山地震、九寨溝地震震中庫侖應(yīng)力變化為正，促進(jìn)了這兩次地震的發(fā)生等。黃祿淵等(2019)基于三維黏彈性有限元模型的計(jì)算結(jié)果認(rèn)為，汶川地震使得蘆山地震與九寨溝地震提前發(fā)生。賈科和周仕勇(2018)基于庫侖應(yīng)力變化及背景地震活動變化的結(jié)果認(rèn)為，汶川地震明顯地觸發(fā)了蘆山地震的發(fā)生，但卻在一定程度上延遲了九寨溝地震的發(fā)生?？梢姴煌瑢W(xué)者得到的研究結(jié)果都顯示了汶川地震使得蘆山地震的庫侖應(yīng)力變化增加，即汶川地震觸發(fā)了蘆山地震的發(fā)生。而關(guān)于汶川地震是否觸發(fā)九寨溝地震，不同研究者得到的結(jié)果仍存在差異。

圖5 青藏高原東緣構(gòu)造背景(震源機(jī)制解數(shù)據(jù)來自于GCMT；https://www.globalcmt.org/)

文中綜合考慮活動構(gòu)造、地球物理等多學(xué)科觀測資料，建立青藏高原東緣的三維黏彈性有限元模型；使用Ji and Hayes(2008)基于遠(yuǎn)場體波資料和有限斷層反演得到的汶川地震同震滑動位移模型，數(shù)值模擬汶川地震造成的同震及震后應(yīng)力場，分析汶川地震對蘆山地震及九寨溝地震的影響。

3.2 青藏高原東緣有限元模型設(shè)置

文中建立了青藏高原東緣的三維黏彈性有限元模型(圖6)。模型東西向?yàn)?00 km；南北向?yàn)?000 km；深度為100 km。根據(jù)活動塊體劃分(鄧起東等，2003; 張培震等，2003)，在橫向上將模型劃分為：青藏高原東北緣、巴顏喀拉塊體、川滇塊體、華南塊體這4個塊體(圖6)。在縱向上，采用半空間分層模型，劃分為4層，包括上地殼、中地殼、下地殼和上地幔層，各分層的深度見表2。模型的上地殼為彈性層；中地殼、下地殼及上地幔為黏彈性層。根據(jù)研究區(qū)的三維速度結(jié)構(gòu)(https://igppweb.ucsd.edu/～gabi/crust2.html；吳建平等，2006; 王椿鏞等，2008; 朱介壽, 2008)計(jì)算得到研究區(qū)各塊體區(qū)域的平均介質(zhì)參數(shù)(楊氏模量、泊松比)：

a—青藏高原東緣的三維有限元模型；b—汶川地震同震破裂模型

表2 青藏高原東緣有限元模型的參數(shù)設(shè)置

其他塊體包括青藏高原東北緣、川滇塊體和巴顏喀拉塊體；上地殼的平均密度設(shè)置為2800 kg/m3，中下地殼上地幔的平均密度設(shè)置為3200 kg/m3

(10)

(11)

其中，E為楊氏模量；υ為泊松比；ρ為平均密度；VS為S波速度；VP為P波速度。各參數(shù)的設(shè)置如表2所示。

中下地殼上地幔的黏度結(jié)構(gòu)對震后形變有著重要的作用。許多學(xué)者基于不同的研究方法對該區(qū)域的中下地殼上地幔的黏度結(jié)構(gòu)做了定量分析(石耀霖和曹建玲, 2008; 張晁軍等，2008; Shao et al., 2011)。如Royden(1996)認(rèn)為青藏高原地區(qū)的中下地殼較為柔軟，并基于解析解，得到其黏度約為1×1020Pa·s。Shao et al.(2011)基于二維黏彈性有限元模型分析了汶川地震的震后形變，認(rèn)為川西地區(qū)的中下地殼黏度約為4×1017Pa·s，四川盆地的中下地殼黏度約為1×1020Pa·s。石耀霖和曹建玲(2008)利用實(shí)驗(yàn)室流變結(jié)果估算了中國大陸巖石圈的等效流變結(jié)構(gòu)，其中青藏高原下地殼等效黏度在1×1019～1×1020Pa·s之間；四川盆地中地殼黏度在1×1021～1×1023Pa·s之間，下地殼黏度1×1021～1×1022Pa·s之間等。雖然不同研究者得到的結(jié)果存在差異，但是相關(guān)研究都顯示了青藏高原的中下地殼較為柔軟，而華南地塊的中下地殼硬度較大。綜合已有研究成果，文中設(shè)置華南地塊的中下地殼上地幔黏度為1×1022Pa·s；其他塊體(包括青藏高原東北緣、巴顏喀拉塊體、川滇地塊)的中下地殼黏度為1×1019Pa·s，上地幔黏度為1×1020Pa·s(表2)。

另外文中假設(shè)地震對遠(yuǎn)場的影響較小，在整個有限元模型的四個側(cè)面邊界施加法向固定，切向自由的位移邊界條件。底部邊界為水平方向可以自由移動，垂向固定。

3.3 汶川地震導(dǎo)致的同震庫侖應(yīng)力變化

汶川地震發(fā)生后，研究者根據(jù)GPS觀測、地震波形記錄等資料反演了地震破裂的滑動位移分布(Ji and Hayes, 2008; 王衛(wèi)民等，2008; Wang et al., 2011)。由于選取的觀測數(shù)據(jù)或反演過程中的參數(shù)選擇等不同，不同研究得到的結(jié)果存在一定的差別。解朝娣等(2010)通過對比5個不同的汶川地震震源模型導(dǎo)致的同震庫侖應(yīng)力變化，發(fā)現(xiàn)其中有兩個模型計(jì)算得到的應(yīng)力變化空間分布與余震的空間分布對應(yīng)較好，Ji and Hayes(2008)研究得到的模型為其中之一。因此，文中使用Ji and Hayes(2008)得到的汶川地震同震破裂模型，在模型的斷層處施加相應(yīng)的同震位移約束，計(jì)算汶川地震導(dǎo)致的同震及震后庫侖應(yīng)力變化。庫侖應(yīng)力變化ΔCSC定義為：

ΔCSC=Δτ+μ′Δσn

(12)

其中，Δτ為斷層面上的剪切應(yīng)力變化，Δσn為斷層面上的法向應(yīng)力變化，μ′為有效摩擦系數(shù)?；诘卣鹞诲e理論可以計(jì)算模型每個節(jié)點(diǎn)的應(yīng)力張量變化；再對應(yīng)力張量進(jìn)行空間坐標(biāo)變換即可得到斷層面的法向應(yīng)力變化和剪切應(yīng)力變化。μ′為有效摩擦系數(shù)，其取值范圍通常在0～0.8之間(King et al., 1994)。ΔCSC大于0，表示巖石更趨近于破裂，斷層的地震危險(xiǎn)性增加；ΔCSC小于0，表示巖石更加遠(yuǎn)離破裂，斷層的地震危險(xiǎn)性減小。

地震導(dǎo)致的應(yīng)力變化是一個張量。計(jì)算庫侖應(yīng)力變化時的應(yīng)力投影方向?qū)Y(jié)果有很大的影響。文中采用將應(yīng)力投影到后續(xù)發(fā)震斷層面上的方式來計(jì)算庫侖應(yīng)力變化。有效摩擦系數(shù)取0.4，計(jì)算深度為10 km。另外，由于文中重點(diǎn)關(guān)注汶川地震對蘆山地震及九寨溝地震的庫侖應(yīng)力影響，因此選取蘆山地震與九寨溝地震的震源附近位置，并繪制其在汶川地震后的庫侖應(yīng)力隨時間的演化(圖7)。圖中顯示汶川地震導(dǎo)致的庫侖應(yīng)力變化在蘆山地震的震源附近為正值。汶川地震導(dǎo)致的庫侖應(yīng)力變化在蘆山地震震源附近增加了約0.013 MPa(其中同震庫侖應(yīng)力變化增加了0.011 MPa；恰在蘆山地震發(fā)生之前的震后庫侖應(yīng)力增加了0.002 MPa)?？梢娿氪ǖ卣饘?dǎo)致的庫侖應(yīng)力變化在蘆山地震震源附近已經(jīng)超過地震觸發(fā)的典型值0.01 MPa(King et al., 1994; Freed, 2005)。即2008年汶川地震造成的同震及震后庫侖應(yīng)力變化觸發(fā)了蘆山地震的發(fā)生。汶川地震導(dǎo)致的庫侖應(yīng)力變化在九寨溝地震的震源附近也為正值(0.009 MPa；圖7)：其中，同震庫侖應(yīng)力變化在九寨溝地震震源附近增加了0.006 MPa；震后庫侖應(yīng)力增加了0.003 MPa。可見汶川地震導(dǎo)致的同震及震后庫侖應(yīng)力變化使得九寨溝地震提前發(fā)生。

a—汶川地震發(fā)生后蘆山地震震源附近的應(yīng)力隨時間的變化；(應(yīng)力投影方向?yàn)樘J山地震震源機(jī)制，走向212°、傾角42°、滑動角100°)；b—汶川地震發(fā)生后蘆山地震震源附近及九寨溝地震震源附近的應(yīng)力隨時間的變化(應(yīng)力投影方向?yàn)榫耪瘻系卣鹫鹪礄C(jī)制，走向150°、傾角78°、滑動角-13°)

另外，文中也計(jì)算了不同的有效摩擦系數(shù)對結(jié)果的影響。分別取有效摩擦系數(shù)為0.0，0.4和0.8，并計(jì)算相應(yīng)的庫侖應(yīng)力變化(圖7)。有效摩擦系數(shù)增大，也就是正應(yīng)力在庫侖應(yīng)力變化中的權(quán)重增大，庫侖應(yīng)力增加的量值減小。反之，有效摩擦系數(shù)減小，正應(yīng)力在庫侖應(yīng)力變化中的權(quán)重減小，庫侖應(yīng)力增加的量值增大(圖7)。但是圖中的結(jié)果顯示，庫侖應(yīng)力變化的趨勢及正負(fù)并沒有隨著有效摩擦系數(shù)的不同而發(fā)生改變。

4 討論

隨著地震位錯理論的進(jìn)一步發(fā)展，現(xiàn)代地震學(xué)的很多研究都基于精細(xì)的地震位錯理論(Sun, 1992; Wang et al., 2003; 王啟欣等，2015)。實(shí)際的地球模型是非常復(fù)雜的；地球介質(zhì)除了有垂向上的分層特征，橫向上也有較大的不均勻性。有限元數(shù)值模擬可以建立更接近實(shí)際情況的地質(zhì)模型，在處理地球介質(zhì)的橫向及縱向不均勻性都具有天然的優(yōu)勢(邵志剛等，2008; 黃祿淵等，2018)。在有限元數(shù)值模擬中，對斷層地震位錯的處理，Melosh and Raefsky(1981)提出了分裂節(jié)點(diǎn)法。該方法將地震引起的位移看作是外力中的一部分；在有限元方程中其他部分不做任何改變的情況下，將斷層引起的地震引入到數(shù)值計(jì)算中，是一種簡單有效且物理意義明確的方法(Melosh and Raefsky, 1981)。

文中分析了走滑斷層實(shí)例的同震及震后效應(yīng)；研究結(jié)果顯示，地球介質(zhì)的橫向不均勻性對地震同震位移有顯著的影響，而中下地殼上地幔的黏彈性松弛對震后效應(yīng)起著主要控制作用。在實(shí)際的震例計(jì)算中，為了分析大地震引起的同震及震后效應(yīng)，應(yīng)該同時考慮地球介質(zhì)的橫向不均勻性及中下地殼上地幔的黏彈性松弛效應(yīng)。最后，根據(jù)上述理論，文中將所開發(fā)的程序應(yīng)用到青藏高原東緣，計(jì)算了2008年MW7.9汶川大地震導(dǎo)致的同震應(yīng)力變化及震后應(yīng)力調(diào)整。研究結(jié)果顯示了汶川地震導(dǎo)致的庫侖應(yīng)力變化在蘆山地震及九寨溝地震的震源附近都為正值(圖7)，說明了汶川地震可能使得蘆山地震與九寨溝地震都提前發(fā)生。文中計(jì)算的結(jié)果與大部分研究得到的結(jié)果是比較一致的(汪建軍和許才軍, 2017; 黃祿淵等，2019)。而也有部分研究指出汶川地震延遲了九寨溝地震的發(fā)生(賈科和周仕勇, 2018)，這可能與計(jì)算地震庫侖應(yīng)力變化時選取的模型參數(shù)、應(yīng)力投影方向等有關(guān)。但是相關(guān)研究都顯示了地震之間的這種相互作用的存在。地震之間的相互作用是地震時空遷移與叢集的重要影響因素(Luo and Liu, 2018; Xu et al., 2020;趙根模等，2020)，定量地研究這種相互作用有助于更進(jìn)一步地分析區(qū)域的地震危險(xiǎn)性。

文中重點(diǎn)討論了地震引起的靜態(tài)的同震效應(yīng)及中下地殼上地幔的黏彈性松弛導(dǎo)致的震后效應(yīng)，在計(jì)算中沒有考慮地震波的動態(tài)傳播、震后余滑等因素的影響，更進(jìn)一步的分析還需要結(jié)合多個方面共同作用的影響。另外，斷層上的構(gòu)造加載是地震孕育發(fā)生的主要因素。在下一步的工作中擬結(jié)合區(qū)域的構(gòu)造加載來模擬歷史地震序列的時空演化，為區(qū)域地震危險(xiǎn)性分析提供參考。

5 結(jié)論

文中開發(fā)了模擬地震同震及震后效應(yīng)的三維黏彈性有限元程序，并分析了地殼橫向不均勻性、中下地殼上地幔的黏度對地震同震及震后效應(yīng)的影響；通過建立青藏高原東緣的三維有限元模型，計(jì)算了汶川地震導(dǎo)致的同震及震后庫侖應(yīng)力變化，探討了汶川地震對蘆山地震和九寨溝地震的影響，得到以下結(jié)論。

(1)地球介質(zhì)的橫向不均勻性對地震同震位移有顯著的影響。如斷層兩側(cè)的楊氏模量不一致，地震同震形變關(guān)于斷層反對稱的性質(zhì)會發(fā)生改變；且楊氏模量大的一側(cè)同震位移較??；楊氏模量小的一側(cè)同震位移較大。

(2)中下地殼上地幔的黏度對大地震震后形變起著主要控制作用。黏度越小，震后形變積累量越大；而黏度越大，震后形變積累量則越小。

(3)汶川地震引起的庫侖應(yīng)力變化在蘆山地震震源附近為0.013 MPa，在九寨溝地震震源附近為0.009 MPa。即汶川地震的發(fā)生有助于蘆山地震和九寨溝地震提前發(fā)生。