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

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

?

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

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

1 中國地震局地震研究所(地震與大地測量重點實驗室)，武漢4300712中國地震局地球物理研究所，北京1000813 日本國際地震學(xué)與地震工程研究所，日本茨城305-0802

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

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

1引言

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

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

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

2研究方法

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

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

(1)

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

因為我們考慮的是動力學(xué)破裂問題，所以需要建立應(yīng)力和斷層面上滑動量的時空分布關(guān)系.當(dāng)接收點接近斷層表面Γ時，斷層面上的應(yīng)力T(x,t)可以表示為

(2)

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

(3)

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

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

T(x,t)=

(4)

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

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

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

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

(5)

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

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

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

(1)格林函數(shù)驗證

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

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

(2)典型算例對比

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

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

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

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

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

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

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

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

3.1不同彎折角斷層的破裂傳播過程

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

4討論與結(jié)論

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

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

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é)概論.北京：地震出版社.

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

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

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

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

(本文編輯胡素芳)

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

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

doi:10.6038/cjg20160320 中圖分類號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.基于邊界積分方程方法的彎折斷層破裂傳播過程控制因素分析.地球物理學(xué)報，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.