袁杰，朱守彪，2

1中國地震局地殼應(yīng)力研究所（地殼動力學(xué)重點實驗室），北京 100085

2中國科學(xué)院計算地球動力學(xué)重點實驗室，北京 100049

1 引言

地震斷裂帶通常是由一系列幾何形態(tài)復(fù)雜、空間上不連續(xù)的斷層段所組成.斷裂帶中，若2條斷層走向平行、中間有一定的間距，這樣的走滑斷層系統(tǒng)稱為斷層階區(qū)（stepover）（Rodgers etal.，1980）.斷層階區(qū)對地震破裂有重要影響：一方面，斷層階區(qū)有可能阻止地震破裂由一條斷層向另一條斷層傳播，成為地震破裂的障礙體（Wesnousky，2006）；但在一定條件下，地震破裂也可能會跳躍斷層階區(qū)，從而形成更大震級的地震（Wesnousky，2006）.由于地震震級與破裂總長度直接相關(guān).所以，若斷層階區(qū)對破裂起阻礙作用，導(dǎo)致地震破裂終止，那么破裂總長度就小，相應(yīng)的地震震級就低；反之，若破裂能夠穿越斷層階區(qū)傳播，則破裂總長度就會大，這樣地震震級就高.因此，定量分析地震破裂能否穿越斷層階區(qū)而形成大地震對于地震學(xué)研究及災(zāi)害評估等有重要的科學(xué)意義及應(yīng)用價值.

實際觀測發(fā)現(xiàn)，不少大地震在破裂過程中，出現(xiàn)了破裂從一個斷層跳躍傳播至另一個斷層，即跨越斷層階區(qū)繼續(xù)傳播的現(xiàn)象，如：2001年的昆侖山口西地震（Antolik etal.，2004），2004年的印尼大地震（Bilham，2005），2008年汶川地震（Shen etal.，2009），2010年智利地震（Lay，2011）和日本大地震（Stein and Okal，2011）等.此外，不少研究人員對斷層階區(qū)影響地震破裂傳播的因素進行了定量分析，如：Harris等（1991）使用有限差分的方法來分析拉張階區(qū)和擠壓階區(qū)之間的不同，分析結(jié)果表明第1條斷層破裂結(jié)束到第二條斷層開始發(fā)生破裂，拉張階區(qū)所需要的間隔時間較擠壓階區(qū)更長.Oglesby（2005）利用有限元方法對走滑斷層階區(qū)間夾著一段傾滑斷層的地震破裂過程進行了模擬，發(fā)現(xiàn)由于這一段傾滑斷層的存在會極大地增加破裂跳躍斷層階區(qū)傳播的可能性，從而導(dǎo)致更大地震的發(fā)生.Finzi和Langer（2012）通過遠場加載及改變斷層端部的動摩擦系數(shù)，來改變斷層階區(qū)的庫侖應(yīng)力狀態(tài)，實現(xiàn)模擬地震破裂越過斷層階區(qū)的傳播情況，發(fā)現(xiàn)軟弱斷層帶及斷層界面不同材料都會促使破裂傳播，從而增加破裂跳躍斷層階區(qū)傳播的可能性.Ryan（2012）利用FaultMod有限元程序模擬破裂穿過斷層階區(qū)的情況，并對比了滑移弱化、速率-狀態(tài)相關(guān)摩擦關(guān)系對破裂跳躍階區(qū)傳播的影響.

盡管前人對斷層破裂跳躍階區(qū)傳播的現(xiàn)象做了大量的研究，但很少有人對影響破裂跳躍階區(qū)傳播的具體因素進行深入的分析.本文利用不連續(xù)變形體接觸力學(xué)的動態(tài)有限單元方法，模擬斷層階區(qū)對地震破裂傳播的控制作用，定量分析初始應(yīng)力場、摩擦關(guān)系及斷層階區(qū)間距等對斷層破裂跳躍階區(qū)傳播的影響.

2 計算原理及有限元模型

2.1 物理方程

2.1.1 運動方程

斷層自發(fā)破裂問題的有限元動力學(xué)方程如下：

上式中ü（t）、u（t）和u（t）分別是系統(tǒng)的節(jié)點加速度向量、節(jié)點速度向量和節(jié)點的位移向量.M，C，K和Q（t）分別是系統(tǒng)的質(zhì)量矩陣、阻尼矩陣、剛度矩陣、和節(jié)點荷載向量.式中：

式中Ft（t）和Ff（D（t），D（t））分別表示構(gòu)造力和斷層面上的接觸力，D表示斷層上的相對滑移距離，D 表示斷層上的相對滑移速率.實際上，模型使用的物理方程與袁杰和朱守彪（2014）所使用的相同，此處不再贅述.

另外，斷層面上節(jié)點的運動狀態(tài)滿足庫侖摩擦準(zhǔn)則（Zienkiewicz etal.，2005）：

上式中σn為法向接觸應(yīng)力，τ為切向應(yīng)力，μ為摩擦系數(shù).

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

斷層上的摩擦關(guān)系非常復(fù)雜（Dieterich，1994；Scholz，1998，2002；Tse and Rice，1986）.前人已嘗試使用多種類型的摩擦關(guān)系來模擬斷層的自發(fā)破裂過程（Lapusta etal.，2000；Shi etal.，2008；Zhang and Chen，2006a，b）.但是，在所有類型的摩擦關(guān)系中，滑移弱化的摩擦本構(gòu)關(guān)系由于其形式簡單且易于實現(xiàn)得到非常廣泛的應(yīng)用.這種摩擦關(guān)系最初是由Ida（1972）、Palmer和Rice（1973）及Andrews（1976a，b）提出.后來，Ohnaka等（1989，1993，1996，1999）也在實驗研究中證實了斷層上的摩擦系數(shù)隨斷層兩側(cè)相對滑移距離增大而減小的真實性.因此本文中的所有模型均采用這個摩擦關(guān)系，其數(shù)學(xué)表達式如下（Andrews，1976a，b；Ida，1972；Palmer and Rice，1973）：

上式中μs為靜摩擦系數(shù)，μd為動摩擦系數(shù)，s為斷層面上兩點之間的相對滑動距離，Dc為特征滑動距離.

2.2 時間積分方法

時間積分方法是非線性有限元計算中的重要方面，為保證計算精度，同時兼顧計算速度，模擬計算中采用中心差分的時間積分方法（顯式有限元算法）.

由于中心差分法是條件穩(wěn)定算法（王勖成，2003），即時間步長Δt必須小于所求解方程性質(zhì)所決定的某個臨界值Δtcr.Δtcr的大小可以通過公式Δtcr＝L／（E／ρ）1／2來近似估算（式中L為模型中尺度最小單元的最小邊長）.經(jīng)過計算，本模型中Δtcr≈0.019s.文中實際計算中所采用的時間步長為1.0×10-3s，遠小于0.019s的臨界值，因而滿足中心差分法的穩(wěn)定性條件.可見，文中的計算過程應(yīng)該是穩(wěn)定、收斂的，能夠滿足計算精度要求.

本文利用ABAQUS／Explicit商業(yè)有限元軟件（有限元顯式計算程序）（Hibbitt etal.，2012）來模擬斷層破裂的動力學(xué)過程，并且使用有限元分析中的接觸單元技術(shù)來描述不連續(xù)斷層的動力學(xué)行為，具體計算中是調(diào)用ABAQUS中的VFRIC用戶子程序進行二次開發(fā)來實現(xiàn)模擬斷層的自發(fā)破裂過程.

2.3 數(shù)值方法驗證

由于斷層的自發(fā)破裂過程是非線性極強的不連續(xù)變形體接觸力學(xué)問題，對于邊界或初始條件稍為復(fù)雜的情況還沒有解析解.因此本文利用非線性有限單元方法來重復(fù)Rojas等（2008）文中的斷層自發(fā)破裂結(jié)果，以驗證本文數(shù)值方法的有效性及可靠性.參照Rojas等的模型，選取驗證模型的空間尺度為60km×60km的正方形，斷層位于模型正中間，其長度為30km.模型中的材料參數(shù)、初始應(yīng)力場、摩擦系數(shù)均與Rojas等模型相同.具體數(shù)值見表1.

表1 模型中的材料參數(shù)、初始應(yīng)力場及摩擦系數(shù)Table 1 Model material parameters，initial stress field and friction coefficients

同Rojas等（2008）的模型一樣，為讓斷層能夠產(chǎn)生自發(fā)破裂行為，模型在斷層中央設(shè)置了成核區(qū).在成核方式上，本文與Rojas等的模型略有不同：本文通過降低成核區(qū)上的摩擦系數(shù)實現(xiàn)成核，發(fā)生自發(fā)破裂行為，而Rojas等的模型通過增加成核區(qū)的初始剪應(yīng)力來實現(xiàn)這一過程.盡管成核方式不同，但實現(xiàn)斷層自發(fā)破裂的效果是等同的.

圖1是模型計算的、斷層上一典型位置（離成核區(qū)中心12.5km處）的相對滑移速率、相對滑移距離以及剪切應(yīng)力隨時間的變化過程.對照圖1給出的相對滑移速率、相對滑移距離、剪切應(yīng)力隨時間變化曲線與Rojas等文中的圖（原文中圖10—12）可以看出，兩者變化格局基本相同，只是某些細節(jié)略有變化.出現(xiàn)兩模型結(jié)果不完全一樣的原因可能是：算法不同（Rojas等（2008）使用的是限差分方法）、成核方式不同以及網(wǎng)格剖分不同等.

總之，通過模型檢驗可以看出，本文所采用的非線性有限元方法完全可以用來模擬斷層的自發(fā)破裂過程，可以定量研究斷層階區(qū)對破裂傳播的影響.

圖1 斷層上一典型點（離成核區(qū)中心12.5km）的相對滑移速率（a）、相對滑移距離（b）、剪切應(yīng)力隨時間的變化（c）Fig.1 Slip velocity（a），sliding distance（b），and shear stress（c）at a typical point on the fault（12.5km from the center of the nucleation patch）vary with time

3 計算結(jié)果

3.1 有限元模型

實際地震的斷層幾何是非常復(fù)雜的（如：2001年發(fā)生的昆侖山口西地震（Antolik etal.，2004），2008年汶川地震（Shen etal.，2009）等），一般會發(fā)生彎曲、拐折、分叉等現(xiàn)象，特別是出現(xiàn)斷層階區(qū).為減少計算量，從而更容易研究斷層階區(qū)對地震破裂傳播過程的控制作用，將實際斷層的三維情況抽象為二維模型，并將斷層簡化為直線（見圖2）.圖2為研究所用的模型幾何及初始、邊界條件.由圖可見，模型空間尺度為150km×150km的正方形，直線AB為斷層1所在位置，直線CD為斷層2所在位置，階區(qū)重疊長度a＝10km，斷層間距為b（其值約為1～5km）.同前人的研究一樣（如：Shi etal.，2008；Olsen-Kettle etal.，2008；袁杰和朱守彪，2014），為讓斷層能夠產(chǎn)生自發(fā)破裂行為，本研究通過成核來實現(xiàn)斷層的自發(fā)破裂過程.模型在斷層中央設(shè)置了成核區(qū)L（見圖2），這段區(qū)域的應(yīng)力狀態(tài)在發(fā)生破裂開始前就已達到或超過破裂準(zhǔn)則所要求的臨界狀態(tài)，其長度Lnucl＝2km.為防止地震波通過邊界反射而影響結(jié)果，模型在四周設(shè)置了吸收邊界（圖2中灰色部分）.此外，在整個區(qū)域內(nèi)施加初始應(yīng)力場并在外部邊界Γ上加載統(tǒng)一的正應(yīng)力σ和剪應(yīng)力τ.

模型全部采用3節(jié)點單元來剖分，以此來完全消除沙漏現(xiàn)象對模型計算產(chǎn)生的影響（Hibbitt etal.，2012）.斷層為研究的重點區(qū)域，為了保證精度，在該區(qū)域?qū)W(wǎng)格進行細化，單元邊長為100m.此外，隨著離開斷層距離的增大，單元的尺度也越來越大，模型最外圍部分單元的邊長為500m.這樣整個模型的有限模型的節(jié)點數(shù)為373338，單元數(shù)為742994.模型四周的吸收邊界采用的是無限元單元（Hibbitt etal.，2012）.

圖2 模型的幾何及初始、邊界條件模型幾何為150km×150km的正方形.圖中直線AB表示斷層1位置，直線CD表示斷層2位置，箭頭所指的粗黑線L為成核區(qū)，其長度為Lnucl＝2km，四周灰色部分為吸收邊界.Fig.2 Model geometry，the initial and boundary conditions The model domain is 100km by 100km.Line AB in the figure is the location of fault 1，line CD is the location of fault 2，the bold black line L to which the arrow points is the nucleation zone，its length is Lnucl＝2km，the gray region around the model is absorbing boundary.

由于本文重點是研究初始應(yīng)力狀態(tài)、摩擦關(guān)系以及階區(qū)間距對斷層破裂跳躍階區(qū)傳播過程的影響，為簡單起見，模型中的介質(zhì)選取為各向同性的線彈性材料.楊氏模量取為E＝7.5×1010Pa、泊松比ν＝0.28和密度為ρ＝2700kg·m-3，由材料的P波速度（α）、S波速度（β）、瑞利波速度（CR）與楊氏模量E、泊松比ν和密度ρ的關(guān)系計算可得α＝5959m·s-1，β＝3294m·s-1，CR＝3045m·s-1.

數(shù)值模擬及實驗室研究表明，斷層的自發(fā)破裂是否向外傳播、破裂傳播速度及破裂模式受應(yīng)力場狀況和摩擦本構(gòu)關(guān)系的影響較大（袁杰和朱守彪，2014）.Wesnousky（2006）指出斷層破裂能否跳躍階區(qū)從一個斷層傳播至另一個斷層，與階區(qū)間距有很大的關(guān)系.因此為研究斷層破裂在經(jīng)過階區(qū)時的傳播情況，本文根據(jù)不同的初始應(yīng)力狀態(tài)、不同的摩擦本構(gòu)關(guān)系以及不同的階區(qū)間距，進行了反復(fù)的數(shù)值實驗，以探究不同因素對斷層破裂跳躍階區(qū)傳播的影響.模型參數(shù)的設(shè)置具體如下：

所有模型中的初始正應(yīng)力均為100MPa；初始剪應(yīng)力分別為54MPa、56MPa、58MPa三種應(yīng)力狀態(tài)；靜摩擦系數(shù)有0.6和0.58兩種；動摩擦系數(shù)有0.5、0.52、0.54、0.56四種；斷層階區(qū)間距有1km、3km、5km三種情況.本文利用控制變量法來研究每一種因素對斷層破裂跳躍階區(qū)傳播的影響，因此對上述幾種因素進行了不同的組合，形成多種模型，通過對比這些模型的模擬結(jié)果來探究這些因素的具體控制作用.

3.2 斷層破裂經(jīng)過斷層階區(qū)的傳播過程

為展示破裂能否跳躍斷層階區(qū)、從一個斷層傳播至另一個斷層的具體過程，首先給出2個具體實例.模型1與模型2的初始應(yīng)力場相同，不同的只是滑移弱化摩擦本構(gòu)關(guān)系中的動摩擦系數(shù).

3.2.1 模型1

模型1中，初始應(yīng)力場為：σ＝100MPa，τ＝54MPa；摩擦本構(gòu)關(guān)系中的參數(shù)：μs＝0.6，μd＝0.5，Dc＝0.1m.

在上述條件下，通過有限元計算發(fā)現(xiàn)，斷層1首先在成核區(qū)（斷層1的中心）開始出現(xiàn)破裂，然后逐步向左右兩側(cè)自發(fā)傳播.圖3是斷層1兩側(cè)的相對滑動距離隨時間的變化.圖中可見，斷層上各點在開始破裂后兩邊的相對滑移距離逐漸增大，斷層上的最大滑移距離為3.74m；滑移以成核中心呈對稱分布，斷層破裂由中心傳至右端終止，歷時8.6s.通過計算，得出平均破裂速度為2907m·s-1，低于S波傳播速度，屬于亞剪切破裂.此外，圖4給出了斷層2兩側(cè)的相對滑移隨時間的變化.從圖中可以看出，從斷層1開始破裂后的8.75s時，斷層2開始從離左邊端部約10km處成核，產(chǎn)生自發(fā)破裂，然后向左右兩側(cè)傳播.且圖4還顯示，當(dāng)斷層1破裂完畢后，斷層2不是立即發(fā)生破裂，而是這中間有一個0.15s的時間延遲.本模型中破裂在跳躍斷層階區(qū)出現(xiàn)的時間延遲現(xiàn)象與Harris等（1991，1993）的結(jié)果一致.產(chǎn)生時間延遲現(xiàn)象的原因，主要是斷層1完全破裂后，盡管斷層面上破裂終止傳播，但斷層上位錯的累積導(dǎo)致斷層端部應(yīng)力在一段時間內(nèi)持續(xù)增大，應(yīng)力波不斷向外傳播.當(dāng)應(yīng)力波到達斷層2時，應(yīng)力在斷層2的斷面上積累，隨著時間的增加，應(yīng)力不斷增大；一旦斷層2的應(yīng)力增大到破裂強度時，斷層2發(fā)生破裂.所以導(dǎo)致了斷層1破裂后，斷層2不是馬上立即就破裂的現(xiàn)象.需要指出的是，斷層2的破裂是被斷層1觸發(fā)的，破裂的成核是自動形成的.

3.2.2 模型2

模型2中，初始應(yīng)力場為：σ＝100MPa，τ＝54MPa；摩擦本構(gòu)關(guān)系中的參數(shù)：μs＝0.6，μd＝0.52，Dc＝0.1m.與模型1不同的只是動摩擦系數(shù)不同.

圖3 模型1中斷層1兩側(cè)相對滑移距離隨時間的變化x軸表示沿斷層距離，y軸表示滑移距離，z軸表示時間；斷層上的破裂由中心不斷向兩邊延伸.斷層破裂由中心傳至整個斷層總共歷時8.6s，在8.6s時斷層上的最大滑移距離為3.74m.Fig.3 Snapshots of the slip profiles along the interface 1various times for Model 1 The xaxis represents the position along the strike，yaxis denotes the slip，and zaxis stand for time.Rupture on the fault continuously extends out from the center of the fault along both opposite directions.Rupture propagates from the center to the end of the fault within 8.6s，and the maximum slip on the fault reaches 3.75mwhen the whole process ends.

圖4 模型1中斷層2兩側(cè)相對滑移隨時間的變化x軸表示沿斷層距離，y軸表示滑移距離，z軸表示時間.在8.75s時，斷層從離左邊端部約10km處開始產(chǎn)生自發(fā)破裂，然后逐步向左右兩側(cè)傳播.由于斷層2左側(cè)較短，因此破裂是不對稱的.Fig.4 Snapshots of the slip profiles along the interface 2various times for Model 1 The xaxis represents the position along the strike，yaxis denotes the slip，and z axis stand for time.The spontaneous generation of rupture begins at 10km from the left end of the fault at the time of 8.75s，then the rupture continuously extends out in both opposite directions.The propagation is asymmetry due to the shorter left side in fault 2.

在上述條件下，模擬結(jié)果顯示，斷層1首先在成核區(qū)開始出現(xiàn)滑移，然后逐步向左右兩側(cè)產(chǎn)生自發(fā)破裂.圖5是斷層1兩側(cè)的相對滑移隨時間的變化情況（類似于圖3）.從圖中可以看出，斷層上各點在開始破裂后兩邊的相對滑移逐漸增大，斷層上的最大滑移量為2.2m；此外，滑移以成核位置為中心呈對稱分布.整個斷層破裂完畢歷時9.2s，通過計算得出平均破裂速度為2717m·s-1，仍低于S波傳播速度，也屬于亞剪切破裂.同樣，圖6給出了斷層2兩側(cè)的相對滑移隨時間的變化.從圖中可以看出，斷層面上的滑移為0，自始至終斷層都未發(fā)生破裂.可見，斷層1上的破裂沒有跳躍傳播至斷層2上.

圖5 模型2中斷層1兩側(cè)相對滑移距離隨時間的變化x軸表示沿斷層距離，y軸表示滑移距離，z軸表示時間；斷層上的破裂由中心不斷向兩邊延伸.斷層破裂由中心傳至整個斷層總共歷時9.2s，在9.2s時斷層上的最大滑移距離為2.2m.Fig.5 Snapshots of the slip profiles along the interface 1various times for Model 2 The xaxis represents the position along the strike，yaxis denotes the slip，and zaxis stand for time.Rupture on the fault continuously extends out from the center of the fault along both opposite directions.Rupture propagates from the center to the end of the fault within 9.2s，and the maximum slip on the fault reaches 2.2mwhen the whole process ends.

圖6 模型2中斷層2兩側(cè)相對滑移距離隨時間的變化x軸表示沿斷層距離，y軸表示滑移距離，z軸表示時間.圖中顯示，斷層上的滑移均為0，可見斷層沒有發(fā)生破裂.Fig.6 Snapshots of the slip profiles along the interface 2various times for Model 2 The xaxis represents the position along the strike，yaxis denotes the slip，and zaxis stand for time.It shows all slips on the fault are 0，so no rupture is generated on the fault.

通過比較以上2個模型，可以看出，在初始應(yīng)力場、靜摩擦系數(shù)和階區(qū)間隔相同時，動摩擦較大者不利于破裂跳躍斷層階區(qū)，此時斷層階區(qū)阻止破裂繼續(xù)傳播，成為地震的終止位置；相反，對于動摩擦系數(shù)較小的斷層，破裂將容易跳躍階區(qū)繼續(xù)傳播，從而有可能導(dǎo)致更大震級的地震發(fā)生.

3.3 決定破裂跳躍斷層階區(qū)的控制因素

3.3.1 影響破裂繼續(xù)傳播的主要因素分析

上面分析了由于動摩擦系數(shù)不同而導(dǎo)致斷層破裂經(jīng)過斷層階區(qū)時的兩種不同情形.為較全面地研究決定破裂跳躍斷層階區(qū)的主要控制因素，下面我們將考慮更多的因素，考察初始應(yīng)力場、摩擦本構(gòu)關(guān)系（動、靜摩擦系數(shù)）、階區(qū)間隔大小對于破裂跳躍傳播過程的影響.計算中模擬了幾十種不同模型進行了詳細分析，表2給出了其中16種典型情況.由于本文主要研究斷層階區(qū)對破裂傳播的影響，因此在這些模型中沒有展示詳細的破裂過程.表2只指示斷層階區(qū)是阻止破裂傳播還是讓其破裂繼續(xù)傳播的情況.

根據(jù)表2的實驗結(jié)果，斷層階區(qū)對破裂過程的影響有以下幾種情況：

表2 不同模型的初始應(yīng)力場、摩擦本構(gòu)關(guān)系及不同階區(qū)間隔時的模擬結(jié)果Table 2 Simulation results in different models in which initial stress fields，friction laws，widths of stepover are different

①比較模型2和模型5，我們可以看出，在其他條件相同時，當(dāng)斷層上的靜摩擦系數(shù)越大，斷層破裂越不易穿越階區(qū).

②比較模型1和模型2、模型7和模型8以及模型14和模型15，可以看出斷層上的動摩擦系數(shù)對破裂能否跳躍階區(qū)的影響一致，其值越大破裂越不易發(fā)生跳躍傳播.

③比較模型6和模型7、模型8和模型9以及模型11和模型12，可以看出，當(dāng)初始剪應(yīng)力越大，即斷層周邊初始剪應(yīng)力與初始正應(yīng)力的比值越大，斷層破裂越容易跳躍階區(qū)傳播.

④比較模型1和模型6、模型3和模型8、模型4和模型15、模型5與模型11、模型7和模型13以及模型10和模型15，可以看出當(dāng)階區(qū)間距越大，斷層破裂越不易跳躍階區(qū)傳播.由此可知，當(dāng)區(qū)域內(nèi)斷層段之間的相隔距離越小，破裂越容易跳躍而發(fā)生大的地震.

由上述分析可見：斷層上的摩擦系數(shù)減小或斷層周邊區(qū)域內(nèi)的初始剪應(yīng)力增大時，都將導(dǎo)致斷層破裂跳躍階區(qū)傳播的可能性.此外，當(dāng)斷層階區(qū)間距越小，斷層破裂也越容易跳躍階區(qū)傳播.

3.3.2 破裂傳播的機制分析

通過以上的對比分析可知，斷層面上的摩擦關(guān)系等模型參數(shù)對破裂通過階區(qū)傳播具有重要的影響.但這是通過模擬實驗總結(jié)出來的一些現(xiàn)象，下面從力學(xué)的角度來定量分析發(fā)生上述現(xiàn)象的主要原因.

由于模型中決定斷層是否破裂是根據(jù)庫侖破裂準(zhǔn)則，因此下面利用庫侖破裂準(zhǔn)則來分析破裂跳躍階區(qū)的傳播機制.

根據(jù)公式（3），兩斷層面開始滑動之前，其斷層上可以承受某一大小的剪應(yīng)力，這種狀態(tài)稱為閉鎖狀態(tài)（粘合）.庫侖摩擦定律定義了一個極限剪應(yīng)力τlim，當(dāng)斷層面上的剪應(yīng)力τ達到τlim時，斷層便開始滑動，這種狀態(tài)稱為滑動狀態(tài)（或破裂）.在本文中，τlim＝μσn，其中，μ表示摩擦系數(shù)，σn表示斷層面上的法向應(yīng)力.

因此根據(jù)庫侖破裂假設(shè)，巖石趨近于破裂狀態(tài)的庫侖破裂應(yīng)力σf為

其中τ為斷層面上的剪應(yīng)力，μs為斷層面上的靜摩擦系數(shù)，σn為斷層面上的法向應(yīng)力.當(dāng)斷層所在區(qū)域內(nèi)的σf＞0時，表明斷層將會發(fā)生破裂，否則斷層處在閉鎖狀態(tài).

圖7 斷層1完全破裂后斷層周邊的庫侖應(yīng)力云圖其中左邊一列圖顯示的是整個模型的庫侖應(yīng)力云圖，右邊一列圖是斷層1右端放大的庫侖應(yīng)力云圖.黑色區(qū)域表示σf＞0，灰色區(qū)域表示σf＜0，最底下一條白線表示斷層1所在位置，其他白線分別表示離斷層1的距離為1km、3km及5km斷層2可能出現(xiàn)的位置.若斷層2出現(xiàn)在黑色區(qū)域里，則破裂可以越過階區(qū)繼續(xù)傳播；否則斷層2不能發(fā)生破裂.Fig.7 The Coulomb stress distribution around the fault 1when the rupture finished completely The left side illustrate the Coulomb of whole model.The right column is the amplified Coulomb stress distribution of the right end of fault 1.The black represent the area whereσf＞0，while the gray denote the area whereσf＜0.White lines on the bottom indicate the location of fault 1，and other yellow lines indicate the possible position of fault 2，which probably is 1km，3km or 5km from fault 1.If fault 2appears in black area，then that the rupture could jump the stepover and extend；if not，the rupture will be arrested.

圖7 續(xù)Fig.7 Continue

圖7 為在斷層1完全破裂后斷層周邊的庫侖應(yīng)力等值線云圖.黑色區(qū)域表示σf＞0，灰色區(qū)域表示σf＜0，根據(jù)庫侖破裂準(zhǔn)則可知當(dāng)斷層1完全破裂后，若斷層2某一部分落入黑色區(qū)域，則根據(jù)庫侖準(zhǔn)則，斷層2將被觸發(fā)產(chǎn)生自發(fā)破裂，即破裂越過斷層階區(qū)繼續(xù)傳播.此外，由圖7可以看出，斷層1破裂后，只有斷層兩端的庫侖應(yīng)力σf＞0，而斷層的其他部分則σf＜0.所以斷層破裂后，應(yīng)力只在斷層端部集中，而其他部位應(yīng)力已經(jīng)釋放.

由圖7c和圖7g可以看出，當(dāng)初始應(yīng)力場和動摩擦系數(shù)相同時，靜摩擦系數(shù)較小的斷層兩端引起的黑色區(qū)域較大，圖7c中黑色區(qū)域未達到離斷層1的距離為1km的地方，而圖7g黑色區(qū)域已達到離斷層1的距離為1km.由此可知當(dāng)初始剪應(yīng)力τ＝54MPa、靜摩擦系數(shù)μs＝0.6、動摩擦系數(shù)μd＝0.54時，斷層1的破裂不能觸發(fā)離斷層1的距離為1km的斷層2的破裂；而當(dāng)靜摩擦系數(shù)μs＝0.58，初始應(yīng)力場、動摩擦系數(shù)相同的情況下，斷層1的破裂能觸發(fā)離斷層1的距離為1km的斷層2的破裂.由以上分析可知在其他因素相同，模型中靜摩擦系數(shù)較小將容易促使斷層破裂跳躍階區(qū)發(fā)生傳播.

由圖7（a、c）、圖7（b、d）及圖7（e、f）可看出，當(dāng)初始應(yīng)力場和靜摩擦系數(shù)相同時，動摩擦系數(shù)較小的斷層兩端產(chǎn)生的黑色區(qū)域較大，即能觸發(fā)斷層2破裂的區(qū)域較大.通過同樣的分析可知，模型中動摩擦系數(shù)較小亦容易促使斷層破裂跳躍階區(qū)發(fā)生傳播.

由公式（4）中可以推斷：當(dāng)摩擦系數(shù)相同時，初始剪應(yīng)力較大將增大σf＞0的區(qū)域范圍.再由圖7（a、b）、圖7（d、e）及圖7（g、h）也可證明這一推斷的正確性.因此可知當(dāng)摩擦系數(shù)相同時，初始剪應(yīng)力較大容易使斷層1的破裂觸發(fā)斷層2的破裂.

以上分析從力學(xué)角度解釋了斷層破裂跳躍階區(qū)的傳播機制，也進一步驗證了表2的結(jié)果.

4 討論與結(jié)論

文中利用有限單元方法對斷層破裂經(jīng)過斷層階區(qū)的傳播情況進行了分析，給出了決定破裂能否跨越斷層階區(qū)繼續(xù)傳播的主要控制因素.但文中的初始應(yīng)力場是均勻的、材料的物性是線彈性，實際情況比文中的模型要復(fù)雜得多.因此，對于非均勻初始場，特別是考慮重力作用及構(gòu)造加載獲得的應(yīng)力場（朱守彪和張培震，2009；Zhu and Zhang，2013），今后應(yīng)作為研究的重點，從而將地震孕育與同震破裂過程進行有機的統(tǒng)一結(jié)合，使得獲得的破裂模型更符合實際.此外，由于應(yīng)力在斷層端部集中，端部介質(zhì)的物性會發(fā)生變化（Duan，2008，2010；Duan and Day，2008；Shibazaki，2005）；斷層兩盤由于長期的地質(zhì)作用，介質(zhì)的物理性質(zhì)會發(fā)生變化（Xia etal.，2005），因此在后續(xù)的研究中應(yīng)考慮更為復(fù)雜的介質(zhì)情況.

根據(jù)上述分析及模擬結(jié)果，得出以下初步結(jié)論：

斷層面上的摩擦系數(shù)減小、斷層周邊區(qū)域內(nèi)的初始剪應(yīng)力增大及斷層階區(qū)間距越小，都會增加斷層破裂跳躍階區(qū)傳播的可能性，斷層破裂也越容易跳躍階區(qū)繼續(xù)傳播，更容易出現(xiàn)更大的地震.相反，斷層上的摩擦系數(shù)較大、初始剪應(yīng)力較小、斷層階區(qū)間隔大，那么斷層階區(qū)可能會阻止破裂繼續(xù)前進，斷層所在位置將是斷層破裂的終止區(qū)域.本研究對于認識破裂跳躍階區(qū)傳播的動力學(xué)過程及進行地震危險性分析有一定的參考價值.

致謝 十分感謝兩位匿名審稿專家所提出的問題以及對文章修改所給予的建設(shè)性意見.

Andrews D J.1976a.Rupture propagation with finite stress in antiplane strain.Journal of Geophysical Research，81（20）：3575-3582.

Andrews D J.1976b.Rupture velocity of plane strain shear cracks.Journal of Geophysical Research，81（32）：5679-5687.

Antolik M，Abercrombie R E，Ekstr?m G.2004.The 14November 2001Kokoxili（Kunlunshan），Tibet，earthquake：Rupture transfer through a large extensional step-over.Bulletin of the Seismological Society of America，94（4）：1173-1194.

Bilham R.2005.A flying start，then a slow slip.Science，308（5725）：1126-1127.

Dieterich J.1994.A constitutive law for rate of earthquake production and its application to earthquake clustering.Journal of Geophysical Research：Solid Earth（1978—2012），99（B2）：2601-2618.

Duan B，Day S M.2008.Inelastic strain distribution and seismic radiation from rupture of a fault kink.Journal of Geophysical Research：Solid Earth（1978—2012），113（B12311），doi：10.1029／2008JB005847.

Duan B C.2008.Effects of low-velocity fault zones on dynamic ruptures with nonelastic off-fault response.Geophysical Research Letters，35（4），doi：10.1029／2008GL033171.

Duan B.2010.Inelastic response of compliant fault zones to nearby earthquakes.Geophysical Research Letters，37（16），doi：10.1029／2010GL044150.

Finzi Y，Langer S.2012.Predicting rupture arrests，rupture jumps and cascading earthquakes.Journal of Geophysical Research：Solid Earth（1978—2012），117（B12303），doi：10.1029／2012JB009544.

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.Geophysical Research Letters，18（5）：893-896.

Harris R A，Day S M.1993.Dynamics of fault interaction：Parallel strike-slip faults.Journal of Geophysical Research：Solid Earth（1978—2012），98（B3）：4461-4472.

Hibbitt H，Karlsson B，Sorensen P.2012.ABAQUS Theory Manual，Version 6.12.Pawtucket，Rhode Island，USA.

Ida Y.1972.Cohesive force across the tip of a longitudinal-shear crack and Griffith′s specific surface energy.Journal of Geophysical Research，77（20）：3796-3805.

Lapusta N，Rice J R，Ben-Zion Y，Zheng G.2000.Elastodynamic analysis for slow tectonic loading with spontaneous rupture episodes on faults with rate-and state-dependent friction.Journal of Geophysical Research，105（B10）：23765-23789.

Lay T.2011.Earthquakes：a Chilean surprise.Nature，471（7337）：174-175.

Oglesby D D.2005.The dynamics of strike-slip step-overs with linking dip-slip faults.Bulletin of the Seismological Society of America，95（5）：1604-1622.

Ohnaka M，Yamashita T.1989.A cohesive zone model for dynamic shear faulting based on experimentally inferred constitutive relation and strong motion source parameters.Journal of Geophysical Research：Solid Earth（1978—2012），94（B4）：4089-4104.

Ohnaka M.1993.Critical size of the nucleation zone of earthquake rupture inferred from immediate foreshock activity.Journal of Physics of the Earth，41（1）：45-56.

Ohnaka M.1996.Nonuniformity of the constitutive law parameters for shear rupture and quasistatic nucleation to dynamic rupture：aphysical model of earthquake generation processes.Proceedings of the National Academy of Sciences of the United States of America，93（9）：3795-3802.

Ohnaka M，Shen L F.1999.Scaling of the shear rupture process from nucleation to dynamic propagation：Implications of geometric irregularity of the rupturing surfaces.Journal of Geophysical Research：Solid Earth（1978—2012），104（B1）：817-844.

Olsen-Kettle L，Weatherley D，Saez E，etal.2008.Analysis of slip-weakening frictional laws with static restrengthening and their implications on the scaling，asymmetry，and mode of dynamic rupture on homogeneous and bimaterial interfaces.Journal of Geophysical Research：Solid Earth（1978—2012），113（B8）.

Palmer A C，Rice J R.1973.The growth of slip surfaces in the progressive failure of over-consolidated clay.Proceedings of the Royal Society of London.A.Mathematical and Physical Sciences，332（1591）：527-548.

Rodgers D A，Balance P E，Reading H G.1980.Analysis of pullapart basin development produced by en echelon strike-slip faults.Sedimentation in Oblique-Slip Mobile Zones，Spec.Publ，4：27-41.

Rojas O，Day S，Castillo J，etal.2008.Modelling of rupture propagation using high-order mimetic finite differences.Geophysical Journal International，172（2）：631-650.

Ryan K.2012.Dynamic models of earthquake rupture on fault stepovers and dip-slip faults using various friction formulations［Ph.D.］.University of Californica.

Scholz C H.1998.Earthquakes and friction laws.Nature，391（6662）：37-42.

Scholz C H.2002.The Mechanics of Earthquakes and Faulting.Cambridge：Cambridge University Press.

Shen Z K，Sun J B，Zhang P Z，etal.2009.Slip maxima at fault junctions and rupturing of barriers during the 2008Wenchuan earthquake.Nature Geoscience，2（10）：718-724.

Shi Z Q，Ben-Zion Y，Needleman A.2008.Properties of dynamic rupture and energy partition in a solid with a frictional interface.Journal of the Mechanics and Physics of Solids，56（1）：5-24.

Shibazaki B.2005.Nucleation process with dilatant hardening on a fluid-infiltrated strike-slip fault model using a rate-and statedependent friction law.Journal of Geophysical Research：Solid Earth（1978—2012），110（B11308），doi：10.1029／2005JB003741.

Stein S，Okal E A.2011.The size of the 2011Tohoku earthquake need not have been a surprise.Eos，Transactions American Geophysical Union，92（27）：227-227.

Tse S T，Rice J R.1986.Crustal earthquake instability in relation to the depth variation of frictional slip properties.Journal of Geophysical Research：Solid Earth（1978—2012），91（B9）：9452-9472.

Wang X C.2003.Finite Element Method（in Chinese）.Beijing：Tsinghua University Press.

Wesnousky S G.2006.Predicting the endpoints of earthquake ruptures.Nature，444（7117）：358-360.

Xia K，Rosakis A J，Kanamori H，etal.2005.Laboratory earthquakes along inhomogeneous faults：Directionality and supershear.Science，308（5722）：681-684.

Yuan J，Zhu S B.2014.FEM simulation of the dynamic processes of fault spontaneous rupture.Chinese J.Geophys.（in Chinese），57（1）：138-156.

Zhang H M，Chen X F.2006a.Dynamic rupture on a planar fault in three-dimensional half space—I.Theory.Geophysical Journal International，164（3）：633-652.

Zhang H M，Chen X F.2006b.Dynamic rupture on a planar fault in three-dimensional half-space—II.Validations and numerical experiments.Geophysical Journal International，167（2）：917-932.

Zhu S B，Zhang P Z.2009.A study on the dynamical mechanisms of the Wenchuan Ms8.0earthquake，2008.Chinese J.Geophys.（in Chinese），52（2）：418-427.

Zhu S B，Zhang P Z.2013.FEM simulation of interseismic and coseismic deformation associated with the 2008Wenchuan Earthquake.Tectonophysics，584：64-80

Zienkiewicz O C，Taylor R L，Zhu J Z.2005.The Finite Element Method：Its Basis and Fundamentals：Its Basis and Fundamentals.Butterworth-Heinemann.

附中文參考文獻

王勖成.2003.有限單元法.北京：清華大學(xué)出版社.

袁杰，朱守彪.斷層自發(fā)破裂動力過程的有限單元法模擬研究.地球物理學(xué)報，2014，57（1）：138-156.

朱守彪，張培震.2009.2008年汶川Ms8.0地震發(fā)生過程的動力學(xué)機制研究.地球物理學(xué)報，52（2）：418-427.