基于邊界積分方程方法的彎折斷層破裂傳播過(guò)程控制因素分析

張麗芬,Bunichiro Shibazaki,廖武林,李井岡,王秋良

1 中國(guó)地震局地震研究所(地震與大地測(cè)量重點(diǎn)實(shí)驗(yàn)室)，武漢　430071　2　中國(guó)地震局地球物理研究所，北京　100081　3 日本國(guó)際地震學(xué)與地震工程研究所，日本茨城　305-0802

?

基于邊界積分方程方法的彎折斷層破裂傳播過(guò)程控制因素分析

張麗芬1,2,Bunichiro Shibazaki3,廖武林1,李井岡1,2,王秋良1

1 中國(guó)地震局地震研究所(地震與大地測(cè)量重點(diǎn)實(shí)驗(yàn)室)，武漢4300712中國(guó)地震局地球物理研究所，北京1000813 日本國(guó)際地震學(xué)與地震工程研究所，日本茨城305-0802

摘要本文利用邊界積分方程方法，以基于三角形網(wǎng)格的全空間格林函數(shù)及離散積分核計(jì)算為基礎(chǔ)，進(jìn)行了最常見(jiàn)的彎折斷層的破裂傳播過(guò)程模擬.為了去除邊界積分方程方法中格林函數(shù)計(jì)算存在的高度奇異性，研究采用分部積分等方法對(duì)動(dòng)力學(xué)方程進(jìn)行了重整化和離散化處理.地震力學(xué)過(guò)程可以被視為斷層由靜摩擦轉(zhuǎn)為動(dòng)摩擦的過(guò)程，對(duì)于震源破裂過(guò)程的動(dòng)力學(xué)模擬，摩擦準(zhǔn)則起著重要作用，本研究采用常用的滑動(dòng)弱化摩擦準(zhǔn)則.計(jì)算引入Courant-Friedrich-Lewy比值來(lái)表達(dá)場(chǎng)點(diǎn)的影響，并控制計(jì)算的收斂性和穩(wěn)定性.通過(guò)與典型算例的比對(duì)，檢驗(yàn)了方法的正確性和有效性.地震破裂能否穿越斷層彎折部位繼續(xù)傳播是震源動(dòng)力學(xué)研究的重要內(nèi)容，基于此，本文建立了多種理論彎折斷層模型，模擬了斷層彎折對(duì)地震破裂傳播的控制作用，并通過(guò)改變斷層周邊初始應(yīng)力場(chǎng)、斷層彎折角度大小以及滑動(dòng)弱化距離大小等來(lái)分析各個(gè)因素對(duì)破裂傳播的影響.模擬結(jié)果表明：斷層面上初始破裂區(qū)域內(nèi)外的應(yīng)力越高，破裂越容易越過(guò)斷層彎折部位繼續(xù)傳播；初始破裂區(qū)域半徑越大，或滑動(dòng)弱化距離越小，破裂也越容易發(fā)生，并越過(guò)彎折部位繼續(xù)傳播.同樣的初始條件，斷層彎折角度越大，斷層彎折作為障礙體，對(duì)破裂傳播的阻礙作用越顯著.小的彎折角，其破裂傳播過(guò)程與平面斷層差別不明顯，基本仍以橢圓方式對(duì)稱向兩側(cè)傳播.

關(guān)鍵詞自發(fā)破裂； 滑動(dòng)弱化摩擦準(zhǔn)則； 三角形網(wǎng)格格林函數(shù)； 彎折斷層； 凹凸體； 障礙體

1引言

作為零階近似，地震斷層可以用簡(jiǎn)單連續(xù)的平面形態(tài)來(lái)模擬.但是，實(shí)際斷層并不是平面構(gòu)造，而是復(fù)雜不連續(xù)的斷裂組合帶，有很多常見(jiàn)的斷層彎折、分叉、階躍、亞平行的分支等(薛霆虓等,2009).理論分析表明，這些斷層幾何結(jié)構(gòu)的非連續(xù)會(huì)引起應(yīng)力場(chǎng)的復(fù)雜變化，對(duì)斷層破裂傳播過(guò)程有重要影響.實(shí)際震例分析發(fā)現(xiàn)，地震破裂的成核、擴(kuò)展以及終止也常常受斷層幾何結(jié)構(gòu)的影響(Sibson,1986;Scholz,1990)，其復(fù)雜性使得變形破壞過(guò)程復(fù)雜化(馬瑾,1987).

Kato等(1999)以及馬瑾等(1996)開(kāi)展了相關(guān)實(shí)驗(yàn)，研究了彎折斷層物理場(chǎng)的演化過(guò)程，并分析了幾種典型斷層與彎折斷層的不同物理場(chǎng)演化特征.Shin和Teng(2001)討論了1999年集集地震斷層幾何形態(tài)造成的斷層附近靜態(tài)位移場(chǎng)變化.Aochi等(2000a,b)研究了二維分支斷層的動(dòng)力破裂過(guò)程，結(jié)果表明，斷層的復(fù)雜幾何形態(tài)造成了斷層面上的應(yīng)力非均勻性，進(jìn)而導(dǎo)致了動(dòng)力學(xué)破裂傳播的多樣化.Aochi和Fukuyama(2002)成功模擬了1992年Landers地震的動(dòng)力學(xué)破裂過(guò)程，由于小的Kickapoo斷層的存在，破裂產(chǎn)生了“跳躍”，使得其傳播方向發(fā)生了改變.Zhang等(2014)利用有限差分方法模擬了由兩條平行斷層Ⅰ和Ⅲ以及一條跨跳斷層II組成的復(fù)雜斷層系的破裂過(guò)程傳播，子斷層I和II的夾角為26.5°.結(jié)果表明，當(dāng)破裂應(yīng)力降較小時(shí)，破裂無(wú)法越過(guò)子斷層Ⅰ和Ⅱ的彎折部位繼續(xù)傳播，而當(dāng)應(yīng)力降較大時(shí)，破裂則可以越過(guò)子斷層II，繼續(xù)向子斷層III傳播.可見(jiàn)，復(fù)雜的非平面斷層幾何形態(tài)確實(shí)對(duì)斷層破裂傳播過(guò)程有重要影響.因此，研究具不同幾何結(jié)構(gòu)的斷層系的變形破壞及演化過(guò)程，對(duì)理解地震機(jī)制有非常重要的意義(馬瑾等,1996).

彎折斷層是一種非常典型常見(jiàn)的斷層幾何結(jié)構(gòu)，一般認(rèn)為，大角度的彎折斷層應(yīng)力難以傳遞，相互作用比較困難，而小角度的彎折斷層相互作用強(qiáng)烈，大震易于連發(fā).而且斷層彎折部位也常常是構(gòu)造地震發(fā)生的區(qū)域，易于導(dǎo)致不同斷層段發(fā)生錯(cuò)動(dòng)，如安縣附近斷層走向的變化，可能是2008 年汶川8.0級(jí)地震造成北川嚴(yán)重?fù)p失的一個(gè)重要的原因.基于上述，本文欲模擬斷層彎折對(duì)地震破裂傳播的控制作用，定量分析初始應(yīng)力場(chǎng)及斷層彎折角度等對(duì)斷層破裂傳播的影響.

2研究方法

2.1三維斷層動(dòng)力學(xué)破裂方程

研究采用邊界積分方程方法，該方法最早由Das和Aki(1997)引入斷層破裂動(dòng)力學(xué)模擬，是一種半數(shù)值半解析方法.對(duì)于斷層的破裂傳播問(wèn)題，基于邊界積分方程方法求解的出發(fā)點(diǎn)是位移表示定理，最終要建立的彈性動(dòng)力學(xué)方程主要是尋找應(yīng)力和滑動(dòng)量之間的關(guān)系.常用的求解途徑有兩種，分別是位移邊界積分方程方法(直接基于位移的積分表示求解)和牽引力邊界積分方程方法(從位移表示定理出發(fā)，通過(guò)求空間導(dǎo)數(shù)得到應(yīng)力的積分表示，并以此為問(wèn)題求解).在牽引力邊界積分方程方法中，積分方程是基于全空間格林函數(shù)建立的，求解的問(wèn)題不限于平面斷層，對(duì)于任意彎曲斷層都成立.因此，從90年代以來(lái)，基于邊界積分方程方法的斷層破裂問(wèn)題研究基本都是采用牽引力邊界積分方程方法(Cochard and Madariaga,1994;Fukuyama and Madariaga,1995,1998;Chen and Aki,1996).

(1)

其中，Uk(x,t)是在接收點(diǎn)x、t時(shí)刻沿著k方向的位移，Δui(ξ,τ)是源點(diǎn)ξ，τ時(shí)刻沿i方向沿?cái)鄬拥幕瑒?dòng)量.Γ是斷層表面，Cijpq是彈性常數(shù)，nj(ξ)是垂直于斷層表面在源點(diǎn)ξ處的法向單位向量，Gkp(x,t-τ;ξ,0)是格林函數(shù)，它表示的是源點(diǎn)ξ一個(gè)沿著p方向的單位體力，將會(huì)在t-τ時(shí)刻，在接收點(diǎn)x處產(chǎn)生沿著k方向相應(yīng)的位移響應(yīng).

因?yàn)槲覀兛紤]的是動(dòng)力學(xué)破裂問(wèn)題，所以需要建立應(yīng)力和斷層面上滑動(dòng)量的時(shí)空分布關(guān)系.當(dāng)接收點(diǎn)接近斷層表面Γ時(shí)，斷層面上的應(yīng)力T(x,t)可以表示為

(2)

但這個(gè)方程僅僅是非常抽象和概念性的表示，如果將該理論應(yīng)用到實(shí)際的地震學(xué)問(wèn)題中，我們需要用更為實(shí)用的邊界積分方程來(lái)描述這種復(fù)雜的破裂傳播問(wèn)題.對(duì)于三維問(wèn)題而言，斷層面垂直于x3坐標(biāo)軸，滑動(dòng)沿著x1和x2軸，也就是說(shuō)滑動(dòng)量有兩個(gè)自由度Δu1和Δu2：

(3)

該方法最大的優(yōu)勢(shì)在于格林函數(shù)存在形式簡(jiǎn)單的閉合解析解，但代價(jià)是在應(yīng)力的積分表示中，由于全空間格林函數(shù)本身當(dāng)場(chǎng)點(diǎn)趨向于源點(diǎn)時(shí)存在弱奇異性，因此對(duì)它的二階導(dǎo)數(shù)會(huì)出現(xiàn)高度奇異性.在動(dòng)力學(xué)研究領(lǐng)域，通常采用分部積分的方法(FukuyamaandMadariaga,1995;Tadaetal.,2000)，將對(duì)格林函數(shù)的高階導(dǎo)數(shù)降低一階到滑動(dòng)量上，從而達(dá)到去除高度奇異性的目的.盡管面積分仍有強(qiáng)奇異性，但利用全空間格林函數(shù)的解析解中含有Dirac函數(shù)的性質(zhì)，最終回避了這個(gè)問(wèn)題，并可以得到積分核的解析表達(dá)式.

采用積點(diǎn)法和箱形離散化方案對(duì)邊界積分方程進(jìn)行離散化處理，網(wǎng)格數(shù)目和時(shí)間步長(zhǎng)數(shù)用NX和NT表示，第n個(gè)網(wǎng)格所在的斷層面用Γn表示，第q個(gè)時(shí)間步長(zhǎng)的結(jié)束點(diǎn)用tq(tq-1≤t≤tq)表示.經(jīng)過(guò)離散化處理后，最終得到了適用于數(shù)值模擬的彈性動(dòng)力學(xué)方程：

T(x,t)=

(4)

破裂過(guò)程的時(shí)空演化由時(shí)間推進(jìn)算法確定，在每一時(shí)間步長(zhǎng)內(nèi)，我們先確定破裂的區(qū)域，評(píng)估破裂尖端的剪應(yīng)力，通過(guò)比較同一位置的剪應(yīng)力與屈服應(yīng)力的關(guān)系判斷破裂是否能擴(kuò)展.

2.2摩擦本構(gòu)關(guān)系

地震力學(xué)過(guò)程可以被視為斷層由靜摩擦轉(zhuǎn)為動(dòng)摩擦的過(guò)程，對(duì)于震源破裂過(guò)程的動(dòng)力學(xué)模擬，摩擦準(zhǔn)則是求解震源動(dòng)力學(xué)方程的關(guān)鍵.它控制著斷裂的起始、發(fā)展和愈合，合理地理解摩擦準(zhǔn)則對(duì)于破裂過(guò)程的控制作用能為我們更好地理解震源過(guò)程的多樣性和非均勻性提供物理基礎(chǔ)(張麗芬和姚運(yùn)生,2013).

目前，較為常用的摩擦準(zhǔn)則是滑動(dòng)弱化準(zhǔn)則(Ida,1972;Andrews,1976a,b)、速率弱化準(zhǔn)則(Cochard and Madariaga,1994;Aagaard et al.,2001)

和速率-狀態(tài)依賴摩擦準(zhǔn)則(Dieterich,1979;Scholz,1998).本研究只考慮破裂的擴(kuò)展和停止，而忽略破裂成核和愈合的過(guò)程，所以采用滑動(dòng)弱化準(zhǔn)則.該準(zhǔn)則最早由Ida(1972)提出，有很好的物理基礎(chǔ)，Dieterich(1979)巖石破裂的實(shí)驗(yàn)結(jié)果表明，破裂點(diǎn)的應(yīng)力降不是瞬時(shí)產(chǎn)生的，而是要經(jīng)過(guò)一定的滑動(dòng)距離，應(yīng)力才能下降到動(dòng)摩擦力，這與滑動(dòng)弱化模型中應(yīng)力與滑動(dòng)距離成反比的假定一致.該準(zhǔn)則數(shù)學(xué)上形式非常簡(jiǎn)單，便于數(shù)值操作，因此在震源動(dòng)力學(xué)的模擬中被廣泛使用(Andrews,1976a,b;Harris et al,1991,張海明,2004).在這個(gè)模型中，假定內(nèi)聚力可以用位移間斷的函數(shù)表示，內(nèi)聚區(qū)的應(yīng)力與破裂面之間的滑動(dòng)距離有關(guān).巖石摩擦實(shí)驗(yàn)表明，隨著滑動(dòng)量的增大，物質(zhì)性質(zhì)變?nèi)?，最終進(jìn)入穩(wěn)態(tài)滑動(dòng)模式.在滑動(dòng)弱化準(zhǔn)則中，破裂尖端是應(yīng)力集中的部位(Ohnaka,1992)，當(dāng)斷層上一點(diǎn)的剪應(yīng)力超過(guò)屈服應(yīng)力強(qiáng)度τp時(shí)，該點(diǎn)開(kāi)始破裂并發(fā)生滑動(dòng).隨著位錯(cuò)Δu的增加，剪應(yīng)力以線性方式下降，直至位錯(cuò)達(dá)到臨界滑動(dòng)弱化距離Dc后，剪應(yīng)力下降到穩(wěn)定的剩余應(yīng)力水平τr：

(5)

2.3時(shí)間推進(jìn)算法及CFL穩(wěn)定性

Courant-Friedrich-Lewy (CFL)比值是控制計(jì)算穩(wěn)定性的重要參量，顯式時(shí)間推進(jìn)算法必須滿足CFL條件(Fukuyama and Madariaga,1998;Tada and Madariaga,2001)，確保在某一給定單元在時(shí)間t產(chǎn)生的波，在t+Δt之前不影響到其他相鄰單元.為保證計(jì)算精度，同時(shí)兼顧計(jì)算速度，模擬計(jì)算中采用中心差分的時(shí)間積分方法.為了保證數(shù)值模擬的穩(wěn)定性，研究中引入了CFL比值ωα=αΔt/Δx來(lái)表達(dá)場(chǎng)點(diǎn)的影響.對(duì)于序號(hào)為(i,j)的場(chǎng)點(diǎn)，只有Δt時(shí)間間隔內(nèi)發(fā)出的波才能影響它，P波速度最大，因此，只有以場(chǎng)點(diǎn)為圓心，αΔt為半徑的圓所覆蓋的單元才對(duì)該場(chǎng)點(diǎn)有貢獻(xiàn).

2.4數(shù)值算法的驗(yàn)證

(1)格林函數(shù)驗(yàn)證

圖1　矩形網(wǎng)格和三角形網(wǎng)格的CFL比值對(duì)當(dāng)前網(wǎng)格貢獻(xiàn)分布示意圖

圖2　應(yīng)力格林函數(shù)計(jì)算結(jié)果對(duì)比(a) 本研究； (b) Tada(2006))；黑色實(shí)線：三角形ABC的格林函數(shù)Lij/k；黑色虛線：三角形CDA的格林函數(shù)Lij/k；灰色虛線：矩形ABCD的格林函數(shù)Lij/k.

(2)典型算例對(duì)比

Madariaga等(1998)基于笛卡爾坐標(biāo)下的交錯(cuò)網(wǎng)格有限差分方法，研究了三維斷層自發(fā)破裂問(wèn)題.在他的研究中，設(shè)定斷層初始破裂區(qū)(凹凸體)的半徑為10個(gè)網(wǎng)格空間步長(zhǎng)，初始剪應(yīng)力Tasp為1.6τp，其中τp為剪切破裂強(qiáng)度，取值80；凹凸體外初始剪應(yīng)力Te為0.5τp；臨界滑動(dòng)弱化距離Dc為4.0，CFL系數(shù)取0.35.圖7分別給出了Madariaga等(1998)和本文的研究結(jié)果，通過(guò)比較過(guò)凹凸體中心切平面上點(diǎn)的滑動(dòng)量隨時(shí)間的變化圖像，不難發(fā)現(xiàn)，二者的滑動(dòng)位移分布基本一致.破裂從斷層中心開(kāi)始，之后以微小的速率向兩側(cè)對(duì)稱傳播.滑動(dòng)量圖像中心部位出現(xiàn)一個(gè)鼓包，這個(gè)過(guò)剩的滑動(dòng)量是破裂起始于一個(gè)有限的凹凸體的最顯著特征(圖3).在Madariaga等(1998)的研究結(jié)果中，當(dāng)t=30時(shí)，滑動(dòng)量的分布圖像中出現(xiàn)了數(shù)據(jù)缺失的情況，初步分析是因?yàn)檠芯空呖紤]到隨時(shí)間演化數(shù)值模擬的不穩(wěn)定性而做了數(shù)據(jù)截?cái)嗵幚硪鸬?本研究中數(shù)值計(jì)算結(jié)果相對(duì)穩(wěn)定，沒(méi)有出現(xiàn)如圖中的缺口，與同樣利用邊界積分方程方法的Aochi(1999)和張海明(2004)的研究結(jié)果基本一致.

本研究還比較了給定時(shí)刻(Madariaga等(1998)中t=200；本研究中t=140)，沿切平面方向的應(yīng)力隨位置的空間變化特征(圖4)，兩種方法的整體形態(tài)比較一致.在破裂前鋒之外出現(xiàn)了類似于Madariaga等(1998)中的強(qiáng)應(yīng)力擾動(dòng)，這種現(xiàn)象在Aochi(1999)、Zhang等(2006)的研究中均有體現(xiàn)，分析認(rèn)為這種強(qiáng)應(yīng)力擾動(dòng)是與S波有關(guān)的.通過(guò)上述比較，檢驗(yàn)了所用研究方法的正確性.

圖3　無(wú)限矩形斷層的自發(fā)破裂滑動(dòng)量對(duì)比(a) 本研究； (b) Madariaga et al.(1998).

圖4　給定時(shí)刻矩形斷層應(yīng)力的空間分布對(duì)比(a) Madariaga et al.(1998)； (b) 本研究.

圖5　彎折角分別為10°、30°、45°和60°的彎曲斷層模型

3不同彎折斷層的自發(fā)破裂傳播過(guò)程模擬

對(duì)研究方法進(jìn)行正確性檢驗(yàn)后，我們考慮利用笛卡爾坐標(biāo)系建立多種彎折斷層模型來(lái)討論其動(dòng)力學(xué)破裂傳播問(wèn)題，主斷層模型的左下角點(diǎn)設(shè)為坐標(biāo)系的原點(diǎn)，整個(gè)斷層系統(tǒng)位于無(wú)限均勻彈性介質(zhì)中.

為研究斷層破裂在經(jīng)過(guò)彎折部位時(shí)的傳播情況，本文根據(jù)不同的初始應(yīng)力狀態(tài)、不同的滑動(dòng)弱化距離以及不同的彎折角度，進(jìn)行了反復(fù)的數(shù)值實(shí)驗(yàn)，利用控制變量法來(lái)研究每一種因素對(duì)斷層破裂越過(guò)彎折部位傳播的影響，并對(duì)上述幾種因素進(jìn)行了不同的組合，形成多種模型，通過(guò)對(duì)比這些模型的模擬結(jié)果來(lái)探究這些因素的具體控制作用.

3.1不同彎折角斷層的破裂傳播過(guò)程

為檢驗(yàn)斷層彎折角度的影響，我們分別建立了平面斷層模型和不同彎折角度的非平面斷層模型進(jìn)行比較，除了斷層模型不同外，初始應(yīng)力條件、滑動(dòng)弱化距離等條件均相同.

1)模型1——彎折角為0°的平面斷層

建立長(zhǎng)度為20 km，寬為5 km、彎折角度為0°的三維垂直走滑平面斷層模型，未考慮斷層厚度，設(shè)置滑動(dòng)沿走向方向.對(duì)斷層模型進(jìn)行離散化處理，網(wǎng)格大小為0.5 km×0.5 km，共劃分400個(gè)三角形網(wǎng)格，時(shí)間步長(zhǎng)取100.設(shè)定S波速度為3.0 km·s-1，縱橫波速比為1.732，剛度模量為30 GPa.對(duì)于自發(fā)破裂傳播問(wèn)題，采用人為指定一個(gè)初始應(yīng)力大于破裂強(qiáng)度的區(qū)域(即凹凸體)，保證破裂首先在該區(qū)域發(fā)生，破裂在摩擦準(zhǔn)則控制下自發(fā)擴(kuò)展.本研究中給定破裂強(qiáng)度τp為1.0 MPa，凹凸體內(nèi)的應(yīng)力為1.2τp，凹凸體外的初始應(yīng)力為0.8τp，剩余應(yīng)力τr設(shè)定為0.

對(duì)于全空間平面斷層而言(圖6)，由于假定斷層面上各處破裂強(qiáng)度是均勻的，因此，破裂在沒(méi)有其他障礙體的情況下，破裂永不停止.實(shí)際模擬結(jié)果可以看到，有限矩形斷層凹凸體內(nèi)的初始應(yīng)力高于破裂強(qiáng)度，使得破裂起始成核.隨著破裂的傳播，裂紋尖端的應(yīng)力集中進(jìn)一步觸發(fā)周圍各點(diǎn)，使得破裂沿走向方向以凹凸體為中心呈橢圓形向兩側(cè)對(duì)稱擴(kuò)展.在此模型設(shè)定的初始應(yīng)力條件下，破裂速度前鋒以亞剪切速度傳播.t=50dt開(kāi)始，切平面方向內(nèi)的破裂前鋒接觸到斷層面的上下邊界(可以認(rèn)為是障礙體)產(chǎn)生停止震相，停止震相繼續(xù)向斷層中心傳播，t=80dt時(shí)，從斷層上下邊界反射回來(lái)的停止震相到達(dá)斷層中心并相遇.隨著應(yīng)力的調(diào)整，滑動(dòng)速率逐漸減小至零，在t=100dt時(shí)，破裂前鋒到達(dá)斷層的左右邊界.

2)模型2——不同彎折角度的彎曲斷層

建立長(zhǎng)度為40 km，寬為10 km、彎折角度分別為10°、30°、45°和60°的三維非平面斷層模型，未考慮斷層厚度，設(shè)置滑動(dòng)沿走向方向.對(duì)斷層模型進(jìn)行離散化處理，網(wǎng)格大小為1 km×1 km，時(shí)間步長(zhǎng)取200.設(shè)定S波速度為3.0 km·s-1，縱橫波速比為1.732，剛度模量為30 GPa.對(duì)于自發(fā)破裂傳播問(wèn)題，仍采用上述人為指定一個(gè)初始應(yīng)力大于破裂強(qiáng)度的區(qū)域，保證破裂首先在該區(qū)域發(fā)生的方法，初始應(yīng)力條件設(shè)置同上述平面斷層.

非平面斷層的數(shù)值模擬結(jié)果顯示，對(duì)于彎折角度為10°的情況(圖7)，彎折角較小，破裂的傳播方 式與平面斷層的傳播方式差別不大，沿?cái)鄬幼呦蛞詸E圓形式近對(duì)稱向兩側(cè)擴(kuò)展.成核后初始破裂階段，滑動(dòng)速率較??；隨著破裂的不斷擴(kuò)展破裂前鋒的滑動(dòng)速率逐漸變大，在t=41 dt，破裂前鋒就到達(dá)了20 km處的斷層彎折部位，但由于彎折角度較小，而且成核區(qū)距離彎折部位很近，所以破裂越過(guò)彎折部位后繼續(xù)傳播.滑動(dòng)速率逐漸增大，到t=141 dt時(shí)最大，之后破裂波前到達(dá)斷層的右邊界，由于反射回來(lái)的停止震相使得滑動(dòng)速率逐漸減小為0，破裂停止.可見(jiàn)，即使在均勻介質(zhì)中，凹凸體外的初始應(yīng)力、破裂強(qiáng)度及Dc都相同的情況下，破裂傳播過(guò)程中滑動(dòng)速率也不是一個(gè)常量.

圖6　平面走滑斷層自發(fā)破裂滑動(dòng)量(左圖，色標(biāo)單位m)和滑動(dòng)速率(右圖，色標(biāo)單位m·s-1)瞬像

圖7　彎折角為10°的彎曲斷層自發(fā)破裂滑動(dòng)量(左圖，色標(biāo)單位m)和滑動(dòng)速率(右圖，色標(biāo)單位m·s-1)瞬像

對(duì)于彎折角度為30°的情況(圖8)，破裂的傳播方式沿?cái)鄬幼呦蛞詸E圓形式近對(duì)稱向兩側(cè)擴(kuò)展.與彎折角為10°的相比，前者的破裂呈似對(duì)稱傳播.隨著彎折角的增大，當(dāng)彎折角為45°時(shí)(圖9)，破裂越過(guò)彎折部位后，破裂的傳播方式與平面斷層的傳播方式出現(xiàn)了明顯差別，在彎折部位不再以近對(duì)稱向兩側(cè)擴(kuò)展，而是出現(xiàn)了明顯的阻礙，破裂在傳播至35 km左右停止.當(dāng)彎折角為60°時(shí)(圖10)，破裂傳播的不對(duì)稱特征更為清晰，大角度彎折對(duì)破裂傳播的阻礙影響更為顯著，破裂在傳播至25 km左右即停止.

通過(guò)上述不同彎折角非平面斷層系統(tǒng)破裂行為的研究發(fā)現(xiàn)，破裂即使在均勻介質(zhì)中沿著斷層彎折之后也可能自動(dòng)停止，強(qiáng)度或應(yīng)力分布的不均勻性對(duì)于破裂的停止不是必須的，顯示了局部的構(gòu)造和斷層的非平面結(jié)構(gòu)在地震產(chǎn)生和破裂過(guò)程中的作用非常重要.對(duì)于彎折斷層而言，較大的彎折對(duì)于破裂的傳播行為有較大影響，可以視為障礙體.彎折之前破裂速度及斷層面上的滑動(dòng)分布均發(fā)生改變，這些主要是由于初始應(yīng)力增加和應(yīng)力降的非均勻性造成的，同時(shí)這又與斷層幾何形態(tài)引起的靜態(tài)應(yīng)力變化相關(guān).薛霆虓等(2009)運(yùn)用有限元法建立大尺度幾何彎曲斷層的二維模型，并利用接觸單元技術(shù)模擬斷層間的作用，探討了具有一定幾何形態(tài)斷層對(duì)斷層系統(tǒng)活動(dòng)的影響.結(jié)果表明，幾何彎曲的斷層導(dǎo)致了應(yīng)力的集中，而且在斷層的地震事件中起到了抑制作用，但是也為孕育大震提供了條件.本研究還發(fā)現(xiàn)，隨著彎折角的增大，在斷層的彎折部位，滑動(dòng)量表現(xiàn)出減小的特征，并且越過(guò)彎折之后，滑動(dòng)量明顯地分成了兩個(gè)部分，這與Aochi(1999)對(duì)彎折斷層的動(dòng)力學(xué)破裂過(guò)程模擬結(jié)果也比較一致，當(dāng)彎折角度較小時(shí)，滑動(dòng)量的分布與平面斷層差別不大，但隨著彎折角的增大，障礙體對(duì)滑動(dòng)量分布的影響更為顯著.

圖8　彎折角為30°的彎曲斷層自發(fā)破裂滑動(dòng)量(左圖，色標(biāo)單位m)和滑動(dòng)速率(右圖，色標(biāo)單位m·s-1)瞬像

圖9　彎折角為45°的彎曲斷層自發(fā)破裂滑動(dòng)量(左圖，色標(biāo)單位m)和滑動(dòng)速率(右圖，色標(biāo)單位m·s-1)瞬像

圖10　彎折角為60°的彎曲斷層自發(fā)破裂滑動(dòng)量(左圖，色標(biāo)單位m)和滑動(dòng)速率(右圖，色標(biāo)單位m·s-1)瞬像

3.2破裂能否越過(guò)斷層彎折部位的其他影響因素

前述分析了由于彎折角度不同而導(dǎo)致斷層破裂經(jīng)過(guò)斷層彎折部位時(shí)的不同情形，為了較全面地研究決定破裂越過(guò)斷層彎折部位的主要控制因素，我們考察了初始應(yīng)力場(chǎng)、摩擦本構(gòu)關(guān)系、滑動(dòng)弱化距離等對(duì)破裂傳播過(guò)程的影響，模擬了十幾種不同的模型，并進(jìn)行了詳細(xì)分析.本文主要研究斷層彎折對(duì)破裂傳播的影響，因此在這些模型中沒(méi)有展示詳細(xì)的破裂過(guò)程，表1給出的10種典型情況，只顯示斷層彎折對(duì)破裂傳播的阻止與否.

① 比較模型1、2、3和模型4，我們可以看出，在其他條件相同時(shí)，滑動(dòng)弱化距離Dc越小，斷層破裂越容易穿越彎折部位，且滑動(dòng)速率越快.對(duì)于模型 1，破裂波前在t=20 dt時(shí)到達(dá)斷層彎折部位，并越過(guò)彎折繼續(xù)傳播，t=80 dt時(shí)破裂前鋒到達(dá)斷層的左右邊界，破裂停止.模型2，破裂波前在t=30 dt到達(dá)斷層彎折部位后繼續(xù)傳播，t=130 dt破裂前鋒到達(dá)斷層的左右邊界，破裂停止.模型3，破裂波前到達(dá)斷層彎折部位的時(shí)間為t=50 dt，在t=200 dt破裂停止.模型4則因?yàn)榛瑒?dòng)弱化距離較大，破裂在成核后不久(t=10 dt)就停止了.

表1　不同初始應(yīng)力場(chǎng)、凹凸體半徑及滑動(dòng)弱化距離對(duì)30°彎折斷層破裂傳播過(guò)程的影響

② 比較模型2、5和模型6，可以看出,凹凸體內(nèi)初始應(yīng)力越高，破裂越容易越過(guò)彎折部位繼續(xù)傳播，且傳播速度越快.

③ 比較模型2、7和模型8，可以看出,凹凸體外初始應(yīng)力越高，破裂越容易越過(guò)彎折部位繼續(xù)傳播，且傳播速度越快.

④ 比較模型2、9和模型10，可以看出,當(dāng)凹凸體半徑越大，破裂越容易越過(guò)斷層彎折部位繼續(xù)傳播，反之破裂不傳播或在很短時(shí)間內(nèi)破裂即停止.

通過(guò)上述模型的對(duì)比，可以對(duì)凹凸體內(nèi)、外的初始應(yīng)力Tasp、Te，凹凸體半徑Rasp和臨界滑動(dòng)弱化距離Dc對(duì)破裂的影響給出一些定性的認(rèn)識(shí).Tasp、Te、Rasp越大，或Dc越小，破裂越容易發(fā)生，也容易越過(guò)斷層的彎折部位繼續(xù)傳播；反之，破裂較難越過(guò)彎折部位，甚至僅在成核后不久破裂即停止.可見(jiàn)，大地震破裂傳播時(shí)間的長(zhǎng)短與上述各因素均存在重要聯(lián)系.當(dāng)凹凸體半徑越大，即震源尺度越大，破裂傳播的時(shí)間相對(duì)較長(zhǎng)；凹凸體(震源體)內(nèi)的初始應(yīng)力超過(guò)屈服應(yīng)力強(qiáng)度越高，破裂越容易發(fā)生，并向外傳播.

4討論與結(jié)論

本文利用邊界積分方程方法對(duì)彎折斷層的破裂傳播情況進(jìn)行了分析，作為三維非平面斷層數(shù)值模擬的例子，我們研究了均勻彈性介質(zhì)中彎折斷層的簡(jiǎn)單情況.本研究的模擬結(jié)果表明，斷層彎折角對(duì)于破裂傳播過(guò)程或是減速或是停止作用，而且由非平面斷層幾何形態(tài)伴隨的動(dòng)態(tài)應(yīng)力改變引起的擾動(dòng)不足以使得破裂傳播停止.文中給出了決定破裂是否能越過(guò)彎折部位繼續(xù)傳播的控制因素，凹凸體內(nèi)、外的初始應(yīng)力Tasp、Te及凹凸體半徑Rasp越大或臨界滑動(dòng)弱化距離Dc越小，破裂越容易發(fā)生，越容易越過(guò)斷層的彎折部位；反之破裂難以發(fā)生，或在成核后不久就停止.實(shí)際的斷層動(dòng)力學(xué)破裂傳播過(guò)程是非常復(fù)雜的，不單單受復(fù)雜斷層幾何形態(tài)的影響，復(fù)雜的介質(zhì)、非均勻應(yīng)力場(chǎng)分布及自由表面的存在等都是重要影響因素.本研究使用的邊界積分方程方法有其局限性，僅能考慮均勻介質(zhì).且在本研究中我們也僅研究了全空間的情況，沒(méi)有將自由表面的影響考慮在內(nèi)，所以后續(xù)的研究會(huì)嘗試其他數(shù)值模擬方法，探討半空間非均勻介質(zhì)中復(fù)雜斷層模型的破裂傳播問(wèn)題，希望將地震孕育與震源破裂過(guò)程進(jìn)行結(jié)合，使獲得的破裂模型更符合實(shí)際.

致謝十分感謝Hideo Aochi老師和Taku Tada老師的耐心指導(dǎo)，并感謝審稿老師所提出的問(wèn)題以及對(duì)文章修改所給予的建設(shè)性意見(jiàn).

References

Aagaard B T,Heaton T H,Hall J F.2001.Dynamic earthquake ruptures in the presence of lithostatic normal stresses: implications for friction models and heat production.Bull.Seism.Soc.Am.,91(6): 1765-1796.

Aki K,Richards P G.1980.Quantitative Seismology Theory and Methods,Vol.I,W.H.San Francisco,California: Freeman and Company.

Andrews D J.1976a.Rupture propagation with finite stress in antiplane strain.J.Geophys.Res.,81(20): 3575-3582.

Andrews D J.1976b.Rupture velocity of plane strain shear cracks.J.Geophys.Res.,81(32): 5679-5687.

Aochi H.1999.Theoretical studies on dynamic rupture propagation along a 3D non-planar fault system [Ph.D.Thesis].Tokyo: University of Tokyo.

Aochi H,Fukuyama E,Matsu′ura M.2000a.Spontaneous rupture propagation on a non-planar fault in 3-D elastic medium.Pure Appl.Geophys.,157(11-12): 2003-2027.

Aochi H,Fukuyama E,Matsu′ura M.2000b.Selectivity of spontaneous rupture propagation on a branched fault.Geophys.Res.Lett.,27(22): 3635-3638.Aochi H,Fukuyama E.2002.Three-dimensional nonplanar simulation of the 1992 Landers earthquake.J.Geophys.Res.,107(B2): ESE 4-1-ESE 4-12.Chen X F,Aki K.1996.An effective approach to determine the dynamic source parameters.Pure Appl.Geophys.,146(3): 689-696.

Cochard A,Madariga R.1994.Dynamic faulting under rate-dependent friction.Pure Appl.Geophys.,142(3-4): 419-445.

Das S,Aki K.1977.A numerical study of two-dimensional spontaneous rupture propagation.Geophys.J.Int.,50(3): 643-668.

Dieterich J H.1979.Modeling of rock friction: 1.Experimental results and constitutive equations.J.Geophys.Res.,84(B5): 2161-2168.

Fukuyama E,Madariaga R.1995.Integral equation method for plane crack with arbitrary shape in 3D elastic medium.Bull.Seism.Soc.Am.,85(2): 614-628.

Fukuyama E,Madariaga R.1998.Rupture dynamics of a planar fault in a 3D elastic medium: rate- and slip-weakening friction.Bull.Seism.Soc.Am.,88(1): 1-17.

Harris R A,Archuleta R J,Day S M.1991.Fault steps and the dynamic rupture process: 2-D numerical simulations of a spontaneously propagating shear fracture.Geophys.Res.Lett.,18(5): 893-896.

Ida Y.1972.Cohesive force across the tip of a longitudinal-shear crack and Griffith′s specific surface energy.J.Geophys.Res.,77(20): 3796-3805.

Kato N,Satoh T,Lei X L,et al.1999.Effect of fault bend on the rupture propagation process of stick-slip.Tectonophysics,310(1-4): 81-99.

Ma J.1987.Introduction of Tectonophysics.Beijing: Seismological Press

Ma J,Ma S L,Liu L Q,et al.1996.Physical field evolution and instability properties of fault geometry.Acta Seismologica Sinica (in Chinese),18(2): 200-207.

Madariaga R,Olsen K,Archuleta R.1998.Modeling dynamic rupture in a 3D earthquake fault model.Bull.Seism.Soc.Am.,88(5): 1182-1197.

Ohnaka M.1992.Earthquake source nucleation: a physical model for short-term precursors.Tectonophysics,211(1-4): 149-178.

Scholz C H.1990.The Mechanics of Earthquakes and Faulting.New York: Cambridge University Press,439.Scholz C H.1998.Earthquake and friction laws.Nature,391:37-42.Shin T C,Teng T L.2001.An overview of the 1999 Chi-Chi,Taiwan,earthquake.Bull.Seism.Soc.Am.,91(5): 895-913.

Sibson R H.1986.Rupture interaction with fault jogs.∥ Earthquake Source Mechanics,Geophys.Monogr.Ser.37.Washington D C: AGU,157-167.

Tada T,Fukuyama E,Madariaga R.2000.Non-hypersingular boundary integral equations for 3-D non-planar crack dynamics.Computational Mechanics,25(6): 613-626.

Tada T,Madariaga R.2001.Dynamic modelling of the flat 2-D crack by a semi-analytic BIEM scheme.Int.J.Numer.Meth.Engng.,50: 227-251.

Tada T.2006.Stress Green′s functions for a constant slip rate on a triangular fault.Geophysical Journal International,164(3): 653-669.

Tada T.2009.Boundary integral equation method for earthquake rupture dynamics.∥ Fault-zone Properties and Earthquake Rupture Dynamics.International Geophysics,94: 217-267.

Xue T X,Fu R S,Lin F.2009.The simulation of activities of bend fault.Chinese J.Geophys.(in Chinese),52(10): 2509-2518,doi: 10.3969/j.issn.0001-5733.2009.10.009.

Zhang H M.2004.Theoretical study on the spontaneous propagation of 3-D seismic rupture on a planar fault in half space[Ph.D.thesis] (in Chinese).Beijing: Peking University.

Zhang L F,Yao Y S.2013.Review on numerical simulation of dynamic rupture process of earthquake source.Acta Seismologica Sinica (in Chinese),35(4): 604-615.

Zhang W B,Iwata T,Irikura K.2006.Dynamic simulation of a dipping fault using a three-dimensional finite difference method with nonuniform grid spacing.J.Geophys.Res.,111(B5): B05301.

Zhang Z G,Zhang W,Chen X F.2014.Three-dimensional curved grid finite-difference modelling for non-planar rupture dynamics.Geophys.J.Int.,199(2): 860-879.

附中文參考文獻(xiàn)

馬瑾.1987.構(gòu)造物理學(xué)概論.北京：地震出版社.

馬瑾,馬勝利,劉力強(qiáng)等.1996.斷層幾何結(jié)構(gòu)與物理場(chǎng)的演化及失穩(wěn)特征.地震學(xué)報(bào),18(2): 200-206.

薛霆虓,傅容珊,林峰.2009.幾何彎曲斷層活動(dòng)性的模擬.地球物理學(xué)報(bào),52(10): 2509-2518,doi: 10.3969/j.issn.0001-5733.2009.10.009.

張海明.2004.半無(wú)限空間中平面斷層的三維自發(fā)破裂傳播的理論研究[博士論文].北京: 北京大學(xué).

張麗芬,姚運(yùn)生.2013.震源動(dòng)力學(xué)破裂過(guò)程數(shù)值模擬研究.地震學(xué)報(bào),35(4): 604-615.

(本文編輯胡素芳)

基金項(xiàng)目國(guó)家自然科學(xué)基金項(xiàng)目(41304046,41404016)和中國(guó)地震局地震研究所所長(zhǎng)基金(IS201506212)共同資助.

作者簡(jiǎn)介張麗芬，女，1981年生，副研究員，主要從事震源動(dòng)力學(xué)及數(shù)字地震研究.E-mail:zhanglf112@163.com

doi:10.6038/cjg20160320 中圖分類號(hào)P315,P631

收稿日期2015-07-14，2015-12-02收修定稿

Controlling factors analysis of dynamic rupture propagation simulation of curved fault based on Boundary integral equation method

ZHANG Li-Fen1,2,Bunichiro Shibazaki3,LIAO Wu-Lin1,LI Jing-Gang1,2,WANG Qiu-Liang1

1KeyLaboratoryofEarthquakeGeodesy,InstituteofSeismology,ChineseEarthquakeAdministration,Wuhan430071,China2InstituteofGeophysics,ChineseEarthquakeAdministration,Beijing100081,China3Internationalinstituteofseismologyandearthquakeengineering,BuildingResearchInstitute,1Tatehara,Tsukuba,Ibaraki 305-0802,Japan

AbstractEarthquakes seldom rupture along single planar faults.Instead,there exist geometric complexities,including fault bends,branches and step overs,which affect the rupture process,nucleation and arrest.In order to understand the influences of nonplanar fault geometry on the earthquake rupture,dynamic numerical simulation provides a new insight.Boundary integral equation method (BIEM) is an appropriate method to model the dynamic rupture process of complex fault geometries,which simplifies the problem and requires small resources in computation by only discretizing the fault surfaces.In addition,it is easy to consider the singularity at the tip of crack using BIEM.When spatially discretizing the fault models,the rectangle meshes are commonly used,however,it is more detailed to describe the nonplanar fault geometries with triangle meshes.It has been recognized that an exponentially growing numerical oscillation resulting from spatiotemporal discretization is well known in the BIEM community.Therefore,an appropriate and optimum combination of space grid and time intervals is very important,which can suppress the oscillation and unstability to some extent.Here,the Courant-Fridrichs-Lewy condition is utilized to achieve the target.In this study,the relationship determined by CFL condition is Vp.

In this paper,the stress Green′s functions for a constant slip rate on a triangular fault are calculated.Theorectially,the mechanics of earthquake rupture process can be regarded as a transformation from static friction to dynamic one.For the dynamic rupture modeling,friction criterion plays a very important role.Because we only focus on the rupture propagation and neglect the nucleation and cessation,slip weakening friction law is applied in our study.Nonplanar fault geometries include many types as mentioned above,and here we just take the curved faults with different bending angles as examples to study the influences of fault geometries on the rupture propagation.After modeling,it is found whether rupture can continue to propagate after curventure part is controlled by many factors,such as bending angles,initial stress in and out of the asperity,the radius of the asperity,slip weakening distance and so forth.Simulation results show that the higher the initial stress in and out of the asperity is,the easier the rupture propagates.And the larger the radius of the asperity or the smaller the slip weakening distance is,the easier the rupture propagates beyond the bend part.Given the same initial conditions,when the inclination angle is bigger,it has more obvious inhibition effect and can be taken as a barrier.However,the curved fault with small inclination angle has similar rupture propagation characteristics,and does not show obvious inhibition effect.

KeywordsSpontaneous rupture propagation; Slip-weakening friction law; Triangular Green′s function; Curved fault; Asperity; Barrier

張麗芬,Bunichiro Shibazaki,廖武林等.2016.基于邊界積分方程方法的彎折斷層破裂傳播過(guò)程控制因素分析.地球物理學(xué)報(bào)，59(3):981-991,doi:10.6038/cjg20160320.

Zhang L F,Bunichiro Shibazaki,Liao W L,et al.2016.Controlling factors analysis of dynamic rupture propagation simulation of curved fault based on Boundary integral equation method.Chinese J.Geophys.(in Chinese),59(3):981-991,doi:10.6038/cjg20160320.