四川長寧6.0級地震庫侖應(yīng)力變化及其對余震的觸發(fā)作用

章 陽 ,陳 婷 ,王利兵 ,羅 娜 ,張朋杰

（1.河北省地震局石家莊中心臺,石家莊 050021;2.河北省地震局預(yù)測研究中心,石家莊 050021;3.河北省地震局紅山基準(zhǔn)臺,河北 邢臺 054000）

0 引言

中國地震臺網(wǎng)中心正式測定,2019年6月17日22時55分在四川宜賓市長寧縣（28.34°N,104.90°E）發(fā)生6.0級地震,震源深度8 km,震源機制解為逆沖兼走滑型地震[1]。本次地震序列較為活躍,截至8月20日,共記錄到3級以上余震64次,其中3.0～3.9級地震52次,4.0～4.9級地震8次,5.0～5.9級地震4次,分別為6月17日23時36分5.1級地震,6月18日7時34分5.3級地震,6月22日22時29分5.4級地震,7月4日10時17分5.6級地震。1900年以來,震中50 km內(nèi)未發(fā)生6級以上地震;震中100 km內(nèi)發(fā)生過2次6級以上地震,震中100～200 km內(nèi)發(fā)生6.0～6.9級地震5次,這些6級以上地震均發(fā)生在川滇塊體東邊界的馬邊—大關(guān)地震帶。此次地震位于川東南地區(qū)北東向華鎣山斷裂的東南方向,發(fā)生在四川盆地邊緣的長寧背斜構(gòu)造附近規(guī)模不大的次級斷層上,根據(jù)主震和余震分布方向,推測發(fā)震斷層為長寧—雙河大背斜[1],這與上述6級以上地震以及發(fā)生在龍門山斷裂帶的汶川8.0級地震和蘆山7.0級地震處于不同的構(gòu)造位置。由于本次地震余震比較豐富,對于余震的危險性評估很有必要,因此,本文綜合考慮同震位錯和震后粘滯松弛效應(yīng),采用粘彈性模型計算主震的庫侖應(yīng)力變化,分析主震對余震的觸發(fā)作用,根據(jù)長時間尺度庫侖應(yīng)力動態(tài)變化特征為預(yù)測后續(xù)余震提供地點依據(jù)。長寧地震區(qū)域地質(zhì)構(gòu)造如圖1。

圖1 長寧地震區(qū)域地質(zhì)構(gòu)造圖Fig.1 Geological structure of Changning earthquake region

靜態(tài)庫侖應(yīng)力變化被應(yīng)用在分析大地震對余震分布的影響中,研究認(rèn)為靜態(tài)庫侖應(yīng)力改變量超過0.01MPa（0.1bar）可以影響余震的發(fā)生位置[2]。庫侖應(yīng)力增強觸發(fā)后續(xù)地震,反之則推遲后續(xù)地震發(fā)生[3],因此,庫侖應(yīng)力在余震危險性評估和地震中長期預(yù)報中應(yīng)用廣泛。國內(nèi)眾多學(xué)者開展利用庫侖應(yīng)力動態(tài)演化分析斷裂帶上地震活動時空變化過程,進一步分析區(qū)域地震危險性。萬永革等[4]計算1900年以來青藏高原東北部發(fā)生的7級以上地震庫侖應(yīng)力積累和演化的過程,并且討論大地震發(fā)生之間的影響,發(fā)現(xiàn)庫侖應(yīng)力對大震的觸發(fā)率達85%;邵志剛等[5]計算汶川8.0級地震引起的庫侖應(yīng)力動態(tài)演化,結(jié)合地震學(xué)和地質(zhì)學(xué)資料分析活動斷裂危險性,并指出在分析強震庫侖應(yīng)力時不能忽略震后效應(yīng);徐晶等[6]研究川滇菱形塊體東邊界斷裂帶內(nèi)庫侖應(yīng)力演化及危險性,認(rèn)為強震之間可能存在觸發(fā)作用,并根據(jù)庫侖應(yīng)力的顯著增加,給出斷裂帶上的強震危險性區(qū);靳志同等[7]根據(jù)前續(xù)斷層在后續(xù)斷層的庫侖應(yīng)力變化給出汶川地震的破裂傳播方向;靳志同等[8]基于均勻彈性半空間模型分析九寨溝7.0地震的同震庫侖應(yīng)力變化對余震的觸發(fā);

然而,基于彈性變形的模擬[9]只能估計應(yīng)力變化的彈性響應(yīng),能模擬地震后較短時間內(nèi)的應(yīng)力應(yīng)變場變化特征。但對于較長的時間間隔（數(shù)十年）,便不能忽略粘彈性松弛效應(yīng)。在地震應(yīng)力轉(zhuǎn)移的應(yīng)力觸發(fā)研究及地震危險性估計過程中,考慮粘彈性應(yīng)力松弛十分重要。程佳等[10]利用粘彈性地殼模型計算歷史地震對九寨溝地震的同震和震后庫侖應(yīng)力作用,認(rèn)為歷史強震促進九寨溝地震的發(fā)生,并計算九寨溝地震對周邊斷層的地震危險性影響;尹鳳玲等[11]基于分層半無限空間粘彈性模型,計算紅河斷裂帶上25次強震由于同震應(yīng)力階變、震后粘滯松弛和震間構(gòu)造應(yīng)力加載綜合作用的庫侖應(yīng)力變化的演化過程,進而分析地震危險區(qū);劉琦等[12]基于位錯理論和粘彈本構(gòu)關(guān)系,分析汶川地震和蘆山地震的同震位錯和震后粘彈松弛效應(yīng)對周邊區(qū)域的地殼形變場、應(yīng)力場、斷層活動狀態(tài)及地震危險性等的影響;黃祿淵等[13]綜合考慮同震靜態(tài)應(yīng)力和震后粘彈性應(yīng)力分析汶川地震對九寨溝地震的影響。Pollitz等[14]認(rèn)為,用Burgers體模型能更好模擬震后的粘滯松弛效應(yīng),計算震后庫侖應(yīng)力變化。Wang等[15]編寫粘彈性應(yīng)力應(yīng)變程序PSGRN/PSCMP程序,該程序基于Burgers粘彈性模型[16],考慮地球重力對形變場的影響,可計算同震和震后長時間尺度的應(yīng)力變化。

因此,本文綜合考慮同震位錯和震后粘滯松弛效應(yīng),基于Burgers粘彈性模型,利用PSGRN/PSCMP程序計算地震破裂面上產(chǎn)生的庫侖破裂應(yīng)力變化,分析主震對3級以上余震的觸發(fā)作用,進而為后續(xù)的余震預(yù)測提供地點依據(jù)。同時,計算粘彈性模型的長時間尺度庫侖應(yīng)力動態(tài)變化,根據(jù)粘彈性模型的震后長時間尺度庫侖應(yīng)力動態(tài)演化特征進一步分析區(qū)域地震危險性。

1 庫侖應(yīng)力變化

根據(jù)庫侖破裂準(zhǔn)則,當(dāng)巖石發(fā)生破裂時,破裂面上的庫侖應(yīng)力（CFS）可定義為

式（1）中,τ為破裂面上的剪應(yīng)力;μ為內(nèi)摩擦系數(shù);σn為正應(yīng)力,張開為正,壓縮為負(fù);p為地殼內(nèi)部的孔隙產(chǎn)生的張性應(yīng)力,S為剪切強度。當(dāng)剪應(yīng)力τ 越接近于S-μ(σn+p),庫侖破裂應(yīng)力CFS越接近0,材料越容易破裂。倘若μ 和S不隨時的變化而變化,那么,庫侖應(yīng)力變化可以表示為

當(dāng)ΔCFS＞0時斷層面處于加載狀態(tài),觸發(fā)后續(xù)地震,反之則推遲后續(xù)地震發(fā)生。研究中常用的庫侖應(yīng)力變化的近似表達形式

其中,μ’=μ(1-β’)為視摩擦系數(shù),β’是斯肯普頓系數(shù),取值范圍為0.7～1.0[17],研究表明,視摩擦系數(shù)不是材料的特質(zhì),主要依靠介質(zhì)的應(yīng)力變化率,不同取值只影響應(yīng)力變化的大小,不影響庫侖破裂應(yīng)力變化的空間分布[18]。在計算庫侖應(yīng)力變化過程中依據(jù)研究區(qū)域的背景資料確定視摩擦系數(shù)的大小,一般取0.2～0.8。本文采用汪榮江研究員開發(fā)的PSGRN/PSCMP程序計算強震的震后庫侖應(yīng)力變化,模擬同震響應(yīng)和震后松弛效應(yīng)。

2 計算模型

2.1 地殼介質(zhì)模型

地殼介質(zhì)模型采用分層半無限空間粘彈性模型,分為完全彈性的上地殼和中地殼、Burgers粘彈性的下地殼和上地幔。本研究中的P波速度參考趙珠等[19]的研究成果,S波速度根據(jù)求得;殼幔分層結(jié)構(gòu)、地層介質(zhì)密度參考王芃等[20]的研究結(jié)果;粘滯系數(shù)參考石耀霖等[21]、李艷娥等[22]的結(jié)果,下地殼的穩(wěn)態(tài)粘滯系數(shù)取為1.0×1021Pa?s,上地幔的粘滯系數(shù)取為1.0×1022Pa?s,下地殼瞬時與穩(wěn)態(tài)粘滯系數(shù)一般相差1～1.5個數(shù)量級[23,11]。因此,下地殼的瞬時粘滯系數(shù)取為1.0×1020Pa?s,。計算庫侖應(yīng)力變化時的具體地殼介質(zhì)模型見表1。

表1 地殼介質(zhì)模型Table 1 Crustal medium model

圖2 余震序列的震源深度分布圖Fig.2 Focal depth distribution map of aftershock sequences

2.2 長寧6.0級地震的震源參數(shù)及同震位錯模型

長寧6.0級地震斷層參數(shù)引用臺網(wǎng)中心的震源機制解：走向323°,傾角57°,滑動角65°。同震位錯模型參數(shù)中,破裂長度和破裂寬度參考程佳等[24]經(jīng)驗公式,計算結(jié)果約為長28 km和寬9 km。走滑同震位錯參考朱航等[25]經(jīng)驗公式,然后根據(jù)斷層面滑動角估算出傾滑位錯,計算得到走滑位錯為0.638 m,傾滑位錯為1.368 m。前人研究認(rèn)為改變摩擦系數(shù)對庫侖應(yīng)力變化影響不大,本文借鑒經(jīng)驗取最優(yōu)值為0.4[26]。

3 計算結(jié)果

3.1 同震庫侖應(yīng)力變化及對余震的觸發(fā)作用

根據(jù)前文給出的地殼介質(zhì)模型和震源參數(shù)及同震位錯模型計算庫侖應(yīng)力變化,接收斷層參數(shù)以該斷層周邊歷史大震的破裂參數(shù)為依據(jù),選取為走向145°,傾角90°,滑動角0°[24]。計算庫侖應(yīng)力變化時,計算深度是一個很重要的控制參數(shù),余震震源深度分布如圖2。采用基于粘彈性模型的PSGRN/PSCMP程序計算不同典型深度處的庫侖應(yīng)力變化,進而分析6.0級主震對3級以上余震的觸發(fā)作用。所得結(jié)果見圖3。

圖3 不同深度處長寧6.0級地震同震庫侖應(yīng)力變化及對3級以上余震的觸發(fā)作用Fig.3 The co-seismic Coulomb stress change of Changning M6.0 earthquake and its triggering effect on aftershocks with M≥3 at different depths

圖3（a）、圖3（c）和圖3（e）分別為5 km深度處、8 km深度處和10 km深度處長寧地震的同震庫侖應(yīng)力變化,主震破裂面上的庫侖應(yīng)力變化分別為0.24 bar、0.9 bar和2.1 bar,均超過應(yīng)力觸發(fā)閾值0.1 bar,有利于觸發(fā)后續(xù)地震,而且隨著計算深度的增加,觸發(fā)作用有一定的增強。同震庫侖應(yīng)力變化均呈現(xiàn)“蝴蝶狀分布特征”,主震區(qū)的破裂方向為北西向。應(yīng)力觸發(fā)區(qū)主要集中在震區(qū)的北北西、南南東、北東東和南西西方向。圖3（b）、圖3（d）和圖3（f）分別為5 km深度處、8 km深度處和10 km深度處長寧地震對3級以上余震的觸發(fā)作用,不同深度計算結(jié)果均表現(xiàn)為64個3級以上余震中,有62個發(fā)生在庫侖應(yīng)力增加的應(yīng)力觸發(fā)區(qū),2個發(fā)生在應(yīng)力影區(qū),綜合觸發(fā)率為96.9%。根據(jù)主震庫侖應(yīng)力變化特征可以較好地預(yù)測3級以上余震的發(fā)震地點,后續(xù)余震震中位置主要集中在震區(qū)的北北西、南南東、北東東和南西西方向。

3.2 庫侖應(yīng)力動態(tài)演化

研究長寧地震長時間尺度的庫侖應(yīng)力動態(tài)演化特征,計算不同深度處Burgers粘彈性模型同震、震后50年、震后100年以及震后200年的庫侖應(yīng)力變化,接收斷層選取為走向145°,傾角90°,滑動角0°。計算5 km深度、8 km深度和10 km深度主震破裂面上震中位置處的不同時間尺度庫侖應(yīng)力變化與同震庫侖應(yīng)力變化的差值,結(jié)果如表2所示。以10 km計算深度為例,分析不同粘滯系數(shù)的選取對計算結(jié)果的影響,結(jié)果如圖4、圖5。

表2 主震破裂面上震中位置處的不同深度不同時間尺度庫侖應(yīng)力變化與同震庫侖應(yīng)力變化差值Table 2 Difference between Coulomb stress changes and co-seismic Coulomb stress changes at different depths and different time scales on the main shock rupture surface of the epicenter location

圖4 粘彈性模型各個時間點庫侖應(yīng)力變化（η1=1.0x1018Pa?s、η2=1.0x1019Pa?s）Fig.4 The Coulomb stress change of viscoelastic model at various time points (η1=1.0x1018Pa?s、η2=1.0x1019Pa?s)

圖5 粘彈性模型各個時間點庫侖應(yīng)力變化（η1=1.0x1020Pa?s、η2=1.0x1021Pa?s）Fig.5 The Coulomb stress change of viscoelastic model at various time points (η1=1.0x1020Pa?s、η2=1.0x1021Pa?s)

由表2可以發(fā)現(xiàn),同一深度處隨著時間尺度的推移,破裂面上的庫侖應(yīng)力變化不斷增大,但幅度逐漸減小,隨著計算深度的增加,庫侖應(yīng)力變化逐漸增加,而且隨著時間尺度的推移,幅度逐漸增大。圖4和圖5為不同粘滯系數(shù)的Burgers模型計算得到的10 km深度處長寧地震同震響應(yīng)、震后50年、震后100年以及震后200年的庫侖應(yīng)力變化圖,均呈現(xiàn)“蝴蝶狀”分布特征。由圖4可以發(fā)現(xiàn),當(dāng)?shù)貙诱硿禂?shù)較小時,應(yīng)力傳遞速度較快,小尺度時間范圍內(nèi)庫侖應(yīng)力變化幅度較大（圖4（a）和圖（b））,隨著時間的推移,同震庫侖應(yīng)力變化將逐漸弛豫耗散,應(yīng)力觸發(fā)區(qū)和應(yīng)力影區(qū)不斷向外圍傳播,震后50年北西向出現(xiàn)新的應(yīng)力影區(qū)（圖4（b）,圖4（c）,圖4（d））,到達穩(wěn)態(tài)所需要的時間短。圖5中地層粘滯系數(shù)較大,應(yīng)力傳遞速度慢,小尺度時間范圍內(nèi)庫侖應(yīng)力變化幅度較?。▓D5（a）和圖5（b））,震后100年北西向出現(xiàn)新的應(yīng)力影區(qū)（圖5（c）,圖5（d））,到達穩(wěn)態(tài)所需要的時間更長,因而百年尺度庫侖應(yīng)力變化幅度較小。依據(jù)長寧地震發(fā)生至震后200年的庫侖應(yīng)力動態(tài)演化（圖4和圖5）,可以看出地震發(fā)生后應(yīng)力觸發(fā)區(qū)主要集中在震區(qū)的北北西、南南東、北東東和南西西方向,其他區(qū)域為應(yīng)力影區(qū)。將庫侖應(yīng)力動態(tài)演化特征應(yīng)用到本區(qū)域的地震危險性分析,重點關(guān)注上述應(yīng)力觸發(fā)區(qū),可為后續(xù)余震發(fā)生提供地點依據(jù)。

4 結(jié)論和討論

本文綜合考慮同震位錯和震后粘滯松弛效應(yīng),利用PSGRN/PSCMP程序計算長寧6.0地震產(chǎn)生的庫侖應(yīng)力變化,分析長寧6.0級地震的庫侖應(yīng)力演化特征及其對3級以上余震的觸發(fā)作用。得出結(jié)論,3級以上余震發(fā)生在應(yīng)力觸發(fā)區(qū)的概率為96.9%,說明3級以上余震與庫侖應(yīng)力變化有較好應(yīng)力觸發(fā)關(guān)系,可以將庫侖應(yīng)力動態(tài)演化應(yīng)用到本區(qū)地震中長期危險性預(yù)測中。粘彈性模型長時間尺度的庫侖應(yīng)力演化呈現(xiàn)“蝴蝶狀”分布特征,震區(qū)北北西、南南東、北東東和南西西向均位于庫侖應(yīng)力增加的應(yīng)力觸發(fā)區(qū),后期重點關(guān)注這些區(qū)域發(fā)生余震的可能性。

在計算庫侖應(yīng)力變化時,存在不確定因素和基本假設(shè)。首先,我們假定長寧6.0地震之前的累積庫侖破裂應(yīng)力為0;其次,建立地殼層狀模型時,由于缺少本地區(qū)精細(xì)的三維速度模型及粘滯系數(shù)等參數(shù),本文建立的模型參數(shù)均是根據(jù)前人研究成果選取的;同時對于初始應(yīng)力場和長期加載場認(rèn)識有限,對于庫侖應(yīng)力估計比較粗糙;這些都會給計算帶來一定誤差。雖不能準(zhǔn)確地預(yù)測破裂時間,但根據(jù)庫侖應(yīng)力動態(tài)演化仍然能夠為地震中長期預(yù)測提供重要信息。