斷層演化過程中的周期和非周期地震滑移及非地震滑移
----一維斷層模型分析*

解孟雨 史保平

(中國北京100049中國科學(xué)院大學(xué)地球科學(xué)學(xué)院)

?

斷層演化過程中的周期和非周期地震滑移及非地震滑移
----一維斷層模型分析*

解孟雨*史保平

(中國北京100049中國科學(xué)院大學(xué)地球科學(xué)學(xué)院)

本文首先根據(jù)Dieterich和Ruina提出的含速率和狀態(tài)的摩擦定律(Dieterich-Ruina定律)，基于一維彈簧-滑塊模型推導(dǎo)了地震復(fù)發(fā)周期的解析表達(dá)式，然后將該近似解與數(shù)值模擬結(jié)果以及Barbot等的相關(guān)研究進(jìn)行了對比分析．此外，本文還利用數(shù)值模擬與理論分析研究了斷層周期和非周期演化的力學(xué)成因機(jī)制以及非地震滑移形成的另類力學(xué)機(jī)制，并討論了一維彈簧-滑塊模型的優(yōu)點(diǎn)及其局限性．結(jié)果表明：① 震后滑移和自加速/成核階段的持續(xù)時間在整個演化過程中不能被忽略；② 在修正后的復(fù)發(fā)周期模型中，復(fù)發(fā)周期的長短除了與斷層特征尺度、作用于斷層面上的有效正應(yīng)力和遠(yuǎn)場加載速率相關(guān)外，還受Dieterich-Ruina定律中摩擦參數(shù)的取值以及臨界滑移距離的影響；③ 當(dāng)給定各個物理參數(shù)和幾何參數(shù)時，目前所得到的解析近似解可以很好地估計地震的復(fù)發(fā)周期，其相對誤差可小于5%；④ 在斷層演化過程中，施加剪切應(yīng)力加載會產(chǎn)生非周期的地震滑移，而在自加速/成核階段后期或震后滑移階段早期，施加較大的剪切應(yīng)力加載，則會出現(xiàn)非地震滑移．

含速率和狀態(tài)的摩擦定律 地震循環(huán) 地震復(fù)發(fā)周期 非地震滑移

引言

華北地區(qū)是中國大陸中強(qiáng)震多發(fā)地區(qū)之一，也是歷史和現(xiàn)代儀器記錄地震相對完整的區(qū)域．華北地震活動區(qū)主要包括4個地震帶：郯城—營口地震帶、華北平原地震帶、汾渭地震帶和銀川—河套地震帶．自公元1300年以來，M≥6的地震有49次，其中包括至少4次M≥8地震．京津唐地區(qū)位于自渤海灣至張家口NW--SE向的華北平原邊緣，該地區(qū)曾發(fā)生過至少14次的M≥6破壞性地震，其中1679 年三河-平谷M8地震是該地區(qū)發(fā)生的最大破壞性地震．20世紀(jì)以來，發(fā)生在該地區(qū)的最強(qiáng)地震為1976年7月28日唐山MS7.8地震，該地震對唐山及其周邊地區(qū)造成了約24萬的人員傷亡以及巨大的財產(chǎn)損失．唐山MS7.8主震之后約15個小時，灤縣又發(fā)生MS7.1地震．1976年11月15日，寧河又發(fā)生MS6.9地震，進(jìn)一步加劇了對唐山及其周邊地區(qū)的破壞程度．鑒于華北地區(qū)中強(qiáng)震位置分布的復(fù)雜性(Liuetal，2012)，40年來關(guān)于唐山大地震的成因機(jī)制以及學(xué)術(shù)界和社會公眾所關(guān)注的地震預(yù)測方法一直是學(xué)術(shù)研究的難題．事實上，由于缺乏必要的觀測手段和對發(fā)震斷層物理/力學(xué)過程的深入了解，目前有關(guān)地震預(yù)測的研究仍然處于初步探索階段．

從地質(zhì)構(gòu)造角度來講，地殼內(nèi)部構(gòu)造地震的發(fā)生與斷層內(nèi)部的剪切形變相關(guān)．彈性回跳理論則給出了斷層摩擦失穩(wěn)的簡單圖像(Stein，Wysession，2003)：在長時間尺度的遠(yuǎn)場剪切應(yīng)力加載下，斷層兩側(cè)的剪切形變能不斷累積，剪切應(yīng)力也隨之增加，當(dāng)剪切應(yīng)力超過斷層本身的摩擦強(qiáng)度時，斷層兩側(cè)突發(fā)相對滑動，部分形變能以地震波的形式傳播至地表，從而造成相應(yīng)的地表運(yùn)動．當(dāng)斷層上的滑動停止后，剪切形變能再次逐步累積，形成了所謂的斷層運(yùn)動“黏滑機(jī)制”，該機(jī)制的提出建立在庫侖(Coulomb)簡單摩擦力學(xué)基礎(chǔ)之上(Byerlee，1970，1978)．由此，Shimazaki和Nakata(1980)將地震重復(fù)發(fā)生(周期和非周期)的過程歸結(jié)為3個理想化模型(圖1)：周期模型、時間可預(yù)測模型和滑動可預(yù)測模型．圖1中應(yīng)力隨時間的變化曲線給出了每個模型相對應(yīng)的震前和震后的應(yīng)力變化(靜態(tài)應(yīng)力降)．基于上述黏滑機(jī)制，對地震重復(fù)發(fā)生過程而言，現(xiàn)今廣泛認(rèn)可的控制因素可歸結(jié)為：未來發(fā)震斷層的構(gòu)造剪切應(yīng)力加載速率、斷層本身的摩擦性質(zhì)和周邊地震對未來發(fā)震斷層的影響(Stein，1999)．如果區(qū)域應(yīng)力場出現(xiàn)了擾動(包含靜態(tài)和動態(tài)應(yīng)力變化及地殼內(nèi)部流體壓力變化等)，則斷層上的應(yīng)力加載或卸載會導(dǎo)致斷層演化周期的改變或斷層摩擦強(qiáng)度的變化，從而引起斷層失穩(wěn)時間的提前或推后．

圖1 地震復(fù)發(fā)模型示意圖(引自Shimazaki，Nakata，1980)(a) 周期模型；(b) 時間可預(yù)測模型；(c) 滑動可預(yù)測模型τ1和τ2分別為地震發(fā)生后和地震發(fā)生前的應(yīng)力，均為常數(shù)(虛線)，Sc為累積同震位移

圖2a給出了受靜態(tài)應(yīng)力擾動后的庫侖模型斷層失穩(wěn)示意圖．可以看出，當(dāng)斷層上累積的剪切應(yīng)力大于或等于斷層摩擦強(qiáng)度時，斷層會突發(fā)失穩(wěn)．因此，當(dāng)給定遠(yuǎn)場剪切應(yīng)力加載速率時，應(yīng)力擾動后斷層失穩(wěn)時間(虛線)的提前或延后量Δt僅與庫侖應(yīng)力變化的絕對大小有關(guān)，而與擾動的加載時刻無關(guān)(Gombergetal，1998; Harris，Simpson，1998; Perfettinietal，2003)．圖2b指出了當(dāng)斷層摩擦強(qiáng)度下降時(虛線)，斷層失穩(wěn)也會出現(xiàn)時間上的提前．考慮到斷層演化的復(fù)雜性，庫侖模型對應(yīng)力擾動下斷層失穩(wěn)時間的描述顯得過于簡單．例如，Harris和Simpson (1998)利用1857年加州M7.9地震和1906年舊金山M7.8地震的位錯數(shù)據(jù)計算了庫侖應(yīng)力的大小，然后估計了應(yīng)力影區(qū)的持續(xù)時間，通過對比兩次地震之后地震發(fā)生的時間和空間信息，結(jié)果表明存在不符合庫侖模型預(yù)測的地震，說明庫侖模型并不能很好地確定地震的延后時間．

圖2 庫侖模型斷層失穩(wěn)示意圖(a) 靜態(tài)應(yīng)力擾動引起的斷層失穩(wěn)(修改自仲秋，史保平，2012)；(b) 斷層摩擦強(qiáng)度降低引起的斷層失穩(wěn)虛線為演化狀態(tài)變化后的斷層演化圖，Δt為斷層演化狀態(tài)變化后第一次失穩(wěn)發(fā)生時間的提前量

對發(fā)震斷層的演化過程，目前地震學(xué)家比較一致的認(rèn)識是斷層滑動至少經(jīng)歷了4個明顯不同的階段(Stein，Wysession，2003)：① 閉鎖階段；② 自加速/成核階段；③ 同震相；④震后松弛/滑移階段．受地震觀測數(shù)據(jù)的限制，我們對每個階段的斷層摩擦動力學(xué)過程及其復(fù)雜性知之甚少．巖石力學(xué)實驗表明，斷層摩擦過程與斷層內(nèi)部的滑移速率和狀態(tài)直接相關(guān)(Marone，1998)．含速率和狀態(tài)的摩擦定律(rate- and state-dependent friction law，簡寫為RSF定律)目前已普遍用于對發(fā)震斷層孕育、成核及發(fā)震深度約束的研究中(Scholz，1998)．在此定律的框架下，有關(guān)斷層演化、地震復(fù)發(fā)周期以及觸發(fā)地震成因機(jī)制的數(shù)值分析與理論研究已在全球范圍內(nèi)得到了廣泛關(guān)注(Perfettini，Avouac，2004a，b; Kanekoetal，2010; Barbotetal，2012; Kameetal，2013a)．而關(guān)于靜態(tài)應(yīng)力擾動后斷層失穩(wěn)時間的提前或推后的論述則最早來自Dieterich (1992，1994)余震觸發(fā)機(jī)制的解析模型．在對帕克菲爾德地震斷層的時空演化特征研究中，基于RSF定律中的Dieterich-Ruina定律(Dieterich，1992，1994；Segall, 2010; Bhattacharya, Rubin, 2014)，Barbot等(2012)給出了地震重復(fù)周期計算的解析表達(dá)式，即

(1)

式中：Δτs為應(yīng)力降；η為描述斷層面形狀的幾何參數(shù)，一般近似等于1(Dieterich，1992)；G為剪切模量；W為斷層面的特征尺度，為斷層長度的一半或斷層面的半徑(圓盤斷層模型而言)；vpl為遠(yuǎn)場加載速率，通常取為常數(shù)；a和b分別為Dieterich-Ruina定律中的參數(shù)；σ為作用于斷層面上的有效正應(yīng)力；vco和vint分別為斷層的同震滑移速率和震間滑移速率．顯然，式(1)可以很好地估算地震重復(fù)周期，但從下文的理論分析與數(shù)值模擬中可以看出，式(1)對于估算斷層的完整演化時間顯然是不完全的，我們將在下文對式(1)進(jìn)行修正．

本文首先利用Dieterich-Ruina定律，從一維彈簧-滑塊模型出發(fā)，結(jié)合不同演化階段所對應(yīng)的近似描述公式，從理論上提出了計算地震復(fù)發(fā)時間的公式；然后利用四階龍格庫塔(4th-order Runga-Kutta)算法進(jìn)行數(shù)值模擬計算，并對比驗證該公式的可行性；另外，本文還將研究發(fā)震斷層的周期和非周期演化的力學(xué)成因機(jī)制以及非地震滑移形成的另類力學(xué)機(jī)制．

1 基本原理

1.1 斷層模型

假設(shè)斷層滿足平面應(yīng)變條件，給定斷層上任意滑動位移分布為δ(x)，那么由此滑動位移所產(chǎn)生的剪切應(yīng)力沿斷層面上的分布τs(x)可表示為(Beeler，2004; Rubin，Ampuero，2005; Bhattacharya，Rubin，2014)

(2)

式中，ν為泊松比．考慮到當(dāng)給定遠(yuǎn)場剪切應(yīng)力τ∞(x)時，作用于斷層面上的剪切應(yīng)力分布τ(x)可表示為

(3)

對式(3)在時間域求導(dǎo)可得

(4)

1.2 摩擦定律

RSF定律是由巖石力學(xué)實驗結(jié)果所總結(jié)出的經(jīng)驗本構(gòu)關(guān)系，該定律定量地刻畫了斷層內(nèi)部的摩擦規(guī)律，也從物理上闡明了斷層內(nèi)部摩擦的復(fù)雜性，已成為研究斷層演化、地震成核、動態(tài)破裂過程等震源物理現(xiàn)象的基礎(chǔ)，并被廣泛用于描述震源行為的系統(tǒng)性變化(Marone，1998)．RSF定律中描述摩擦系數(shù)的一般性方程為(Guetal，1984; Segall，2010)

(5)

式中：τ為作用在斷層面上的剪切應(yīng)力；θi為狀態(tài)變量(i=1，2，3，…，N)，是任何可以刻畫接觸面狀態(tài)的變量，如顆粒大小、斷層泥的孔隙度等(Segall，2010)；μ為斷層面上的摩擦系數(shù)，與庫侖摩擦過程不同的是，在RSF定律中，μ為斷層滑移速率v和狀態(tài)變量θi的函數(shù)．對于θi的演化，若保持正應(yīng)力σ不變，則滿足以下關(guān)系：

(6)

式中，i=1，2，3，…，N．式(5)和式(6)組成了完整的RSF定律(Guetal，1984; Segall，2010)．

1.3 Dieterich-Ruina定律

為使計算和分析更簡單方便，通常采用單一狀態(tài)變量θ的RSF定律．目前最常用的兩個定律分別是Dieterich-Ruina定律和Ruina定律(Segall，2010; Bhattacharya，Rubin，2014)．在穩(wěn)定狀態(tài)或接近穩(wěn)定狀態(tài)時，兩個定律的公式形式是一致的，但在其它狀態(tài)下兩者是不同的(Bhattacharya，Rubin，2014)．Dieterich-Ruina定律可以很好地解釋滑動—約束—滑動測試中摩擦界面的愈合現(xiàn)象，但不能解釋速度階梯變化實驗中應(yīng)力演化到新的穩(wěn)定狀態(tài)時所需的位移與速率階梯變化大小和正負(fù)符號無關(guān)；與之相反，Ruina定律則可以解釋速度階梯變化實驗中的應(yīng)力變化現(xiàn)象，但不能闡明上述愈合現(xiàn)象(Bhattacharya，Rubin，2014)．Nagata等(2012)針對Dieterich-Ruina定律提出了修正的摩擦定律，可以同時解釋上述滑動—約束—滑動測試和速度階梯變化實驗中出現(xiàn)的現(xiàn)象．為便于討論，考慮到摩擦定律的使用廣泛程度，本文選用Dieterich-Ruina定律進(jìn)行數(shù)值模擬與理論推導(dǎo)，其表達(dá)式為

(7)

(8)

式中：vo為參考速率；μo為穩(wěn)定滑動且滑動速率為vo時的摩擦系數(shù)；θ為凹凸體的平均接觸時間，即凹凸體之間接觸的平均存在時間(Scholz，2002)；Dc為臨界滑動距離，即徹底改變斷層接觸面摩擦狀態(tài)所需的滑動距離(Marone，1998)；a和b分別為實驗常數(shù)，a為“直接影響”系數(shù)，決定速率變化所引起的摩擦強(qiáng)度變化，b為“演化影響”系數(shù)，控制狀態(tài)演化所引起的摩擦系數(shù)變化(Bhattacharya，Rubin，2014)，一般來講，a和b均大于0，且具有相同的數(shù)量級，其取值范圍為0.005—0.01(Dieterich，1994)．在剪切作用下，假定作用于斷層面上的正應(yīng)力始終保持不變，對式(7)在時間域求導(dǎo)，可得

(9)

(10)

顯然，當(dāng)a-b≥0時，斷層的摩擦過程處于速度強(qiáng)化狀態(tài)，此時斷層的滑動是穩(wěn)定的；而當(dāng)a-b<0時，斷層的摩擦過程處于速度弱化狀態(tài)，斷層的滑動可以是條件穩(wěn)定，也可以是不穩(wěn)定的(Scholz，1998)．

結(jié)合式(4)與式(9)，斷層摩擦過程的數(shù)學(xué)描述可表示為

(11)

該式被廣泛應(yīng)用于斷層成核和自發(fā)破裂過程的研究中．對于單一自由度的一維彈簧-滑塊模型(圖3)，滑塊加載過程可由下式給出：

τ(t)=τ(0)+k(vplt-δ)，

(12)

式中：τ(0)為斷層內(nèi)部保留的初始應(yīng)力；t為時間變量；k為彈簧的有效剛度系數(shù)，可表示為k=ηG/W．因此，對于一維彈簧-滑塊模型，式(11)可簡化為

(13)

一維彈簧-滑塊模型看似簡單，但事實上包含了斷層摩擦滑移破裂過程的3個基本要素：① 滑塊的摩擦面，其代表了斷層面；② 彈簧，k的取值代表了地殼的彈性性質(zhì)和幾何特征，彈簧用以積累彈性應(yīng)變能，為斷層自加速階段提供必要的應(yīng)變能；③ 穩(wěn)定的遠(yuǎn)場加載速率，其模擬了穩(wěn)定的板塊運(yùn)動，為斷層應(yīng)變能的積累提供動力來源(Segall，2010)．

圖3 一維彈簧-滑塊模型

當(dāng)vpl恒定時，給定初始的速度v(0)和狀態(tài)變量θ(0)，通過聯(lián)立方程(8)與(13)，并采用數(shù)值計算方法求解方程組，即可模擬斷層在時間域內(nèi)的演化過程，其中當(dāng)v→∞時，斷層發(fā)生破裂失穩(wěn)．為了避免v→∞這一非物理現(xiàn)象的發(fā)生，通常在數(shù)值求解過程中加入輻射衰減項 (Segall，2010)，如對式(12)在時間域內(nèi)求導(dǎo)并加入衰減項后可得：

(14)

或

(15)

2 數(shù)值解

結(jié)合式(8)與式(15)，采用四階龍格庫塔算法可以求出在給定不同參數(shù)條件下一維斷層的演化特征．表1給出了模擬計算中所使用的參數(shù)值．

表1 模擬計算中使用的參數(shù)值(引自Kame et al, 2013b)

注：狀態(tài)變量初始值θ(0)由v(0)和μ(0)求得.

圖4給出了加入輻射衰減項后模擬得到的包括滑移速率和應(yīng)力變化在內(nèi)的斷層演化過程．可以看出，在模擬計算開始的0—160年(從起始時刻至第一次失穩(wěn)的時間段)內(nèi)，計算結(jié)果與之后時間段的結(jié)果并不相同．在上述時間段內(nèi)斷層的演化受初始參量v(0)，θ(0)和介質(zhì)參數(shù)的控制，所以改變v(0)和θ(0)會使得從起始時刻至第一次失穩(wěn)時間段內(nèi)的演化曲線發(fā)生明顯變化，不過改變v(0)和θ(0)卻不會影響之后時間段的演化情況．實際上第一次失穩(wěn)后的演化僅受介質(zhì)參數(shù)的影響，可以反映介質(zhì)參數(shù)所控制的一維斷層模型本身的演化特征，因此之后分析中不考慮0—160年內(nèi)的演化情況．

圖4 基于一維斷層模型的滑移速率周期演化圖(a)和相平面內(nèi)的閉合軌跡圖(b)Tr為斷層模型的完整演化周期，黑色三角形為各個階段的分界標(biāo)志，下同

圖4a給出了滑移速率的變化特征．可以看出，滑移速率隨時間的變化呈現(xiàn)出明顯的周期性特征，斷層在失穩(wěn)時刻的滑移速率可達(dá)0.14 m/s，圖中曲線所呈現(xiàn)出的特征均與斷層運(yùn)動的黏滑機(jī)制一致．圖4b給出了地震循環(huán)的4個主要階段．模擬結(jié)果顯示：斷層閉鎖階段至少占總的演化周期Tr的70%，而自加速階段、同震相和震后松弛/滑移階段這3個階段只占了30%Tr．

圖5為一個地震循環(huán)周期內(nèi)的一維斷層演化圖．由圖4和圖5可以看出，斷層演化在時間上由4個主要階段組成，每個階段的持續(xù)時間是不同的. 如果以Tint，Tpre，Tco和Tpost分別代表斷層演化過程中閉鎖、自加速/成核、同震相和震后松弛/滑移階段所經(jīng)歷的時間尺度，那么完整的斷層演化周期Tr可表示為

Tr=Tint+Tpre+Tco+Tpost.

(16)

顯然，斷層所處的閉鎖階段的持續(xù)時間遠(yuǎn)遠(yuǎn)大于其它3個階段持續(xù)時間的總和(圖4和圖5), 而成核階段和震后滑移階段的持續(xù)時間基本相當(dāng)，即Tpre≈Tpost．與其它3個階段相比，同震滑移的持續(xù)時間Tco一般為幾秒至幾分鐘．另外，對應(yīng)于每個階段，Ω=vθ/Dc值的大小表現(xiàn)出明顯差異：當(dāng)斷層處于震后滑移和閉鎖階段時，vθ/Dc<1；當(dāng)斷層處于自加

圖5 一個地震循環(huán)周期內(nèi)的一維斷層演化圖(a) 斷層滑移速率演化圖；(b) 斷層模型中Ω=vθ/Dc隨時間的變化

速/成核階段時，vθ/Dc>1．如上所述，同震滑移在幾秒到幾分鐘內(nèi)完成，即Tco相對很小，因此在地震復(fù)發(fā)周期的計算中Tco可以忽略不計．由此，地震演化周期Tr可寫為

Tr≈Tint+Tpre+Tpost.

(17)

3 斷層演化周期的近似計算與相應(yīng)數(shù)值解的對比

3．1 解析近似解

當(dāng)斷層演化處于自加速階段并滿足vθ/Dc?1時，則Dieterich-Ruina定律中的狀態(tài)方程，即式(8)可近似為

(18)

結(jié)合式(13)，則滑移速率的近似解為(Dieterich，1992，1994)

(19)

(20)

(21)

震后滑移開始于同震滑移的結(jié)束，因此在震后滑移起始時刻，滑移速率v近似等于vco，θ近似等于θco(Rubin，Ampuero，2005)．如果設(shè)定震后滑移階段始于t=0，則震后滑移和閉鎖階段的滑移速率隨時間變化的近似解為(Rubin，Ampuero，2005)

(22)

狀態(tài)方程的解為

θ=θco+t.

(23)

(24)

式中，tp為震后滑移的持續(xù)時間尺度(震后滑移階段始于t=0)，也是震后滑移的結(jié)束時刻和閉鎖階段的開始時刻．因此，我們稱tp為斷層演化過程中由震后滑移到閉鎖階段的轉(zhuǎn)折

圖6 滑移速率v(a)和狀態(tài)變量θ(b)的近似解與數(shù)值解對比圖

點(diǎn)．震后滑移所持續(xù)的時間尺度可近似為

(25)

(26)

(27)

對式(27)兩邊取對數(shù)后合并同類項，可得

(28)

那么則有

(29)

綜合上述所得到的Tpre，Tpost和Tint，則斷層演化周期的近似解可寫為

(30)

(31)

其具體過程請參見附錄.

比較式(30)與式(31)可以看出，Barbot等(2012)對地震循環(huán)周期的估計忽略了震后滑移和自加速/成核階段所需的演化時間．在Dieterich(1994)的余震觸發(fā)模型中，ta也為主震后余震的持續(xù)時間．Stein和Liu(2009)的研究也表明，在較小的加載速率下，例如大陸內(nèi)部構(gòu)造地震，余震的持續(xù)時間可能高達(dá)幾百年．仲秋和史保平(2012)對1976年7月28日唐山MS7.8地震余震序列的分析表明，如果余震的觸發(fā)機(jī)制可由Dieterich模型表示，那么該地震的余震持續(xù)時間ta大約為80—100年．在Dieterich-Ruina定律中，a和b的取值處于相同的數(shù)量級，如果斷層能夠完成自加速/成核過程，則必須有b>a，由此推斷震后滑移階段的持續(xù)時間Tpost可能會大于自加速/成核階段的持續(xù)時間Tpre，且應(yīng)處于相同的數(shù)量級．因此，式(31)等號右側(cè)的第一、三項在斷層演化周期的估算中是不能省略的．關(guān)于閉鎖階段滑移速率vint的估算，Rubin 和Ampuero(2005)給出了相應(yīng)的近似解．對于一維斷層模型，vint主要受k和kc的影響，一般來講vint比遠(yuǎn)場加載速率vpl小2—3個數(shù)量級(Rubin，Ampuero，2005)．

3.2 近似解與數(shù)值模擬結(jié)果的對比

圖7 不同a/b值情況下斷層滑移速率和剪切應(yīng)力的相圖

表2 不同a/b值情況下近似解與數(shù)值解的對比

Table 2 Comparison of analytic approximate solutions with numerical solutions for differenta/bvalues

a/bTpost/a數(shù)值解近似解Tint/a數(shù)值解近似解Tpre/a數(shù)值解近似解Tco/sTr/a數(shù)值解近似解ta/a2πta/a0.43.9054.10655.09354.0195.0624.7751.2864.06062.9011.64210.3200.54.9144.92753.87654.0197.0336.5501.6165.82465.4962.46415.4800.66.1916.15952.75854.0199.6398.8062.0468.58868.9843.69623.2200.78.2318.21251.47854.01913.39111.8592.6273.09974.0915.74936.1190.812.28112.31849.54454.01919.63816.4073.6081.46382.7449.85561.9190.923.56324.63744.47354.01934.23724.4565.71102.273103.11222.173139.317

D'cDcTpost/a數(shù)值解近似解Tint/a數(shù)值解近似解Tpre/a數(shù)值解近似解Tco/sTr/a數(shù)值解近似解ta/a2πta/a10100.925100.786505.352540.192163.818141.25531.18770.094782.23376.150478.463550.26950.393252.878270.09681.90070.62715.49385.047391.11738.075239.232220.19420.157101.064108.03832.76128.2516.20154.019156.44715.23095.693110.05910.07950.57254.01916.37814.1263.1177.00978.2237.61547.8460.55.0335.03925.27627.0108.1957.0631.5538.50539.1123.80823.9230.22.0122.01610.11410.8043.2762.8250.6215.40215.6451.5239.5690.11.0051.0085.0575.4021.6381.4130.317.7017.8220.7624.785

(32)

(33)

4 剪切應(yīng)力加載作用下的斷層非周期地震滑移

當(dāng)斷層在演化過程中受到應(yīng)力擾動后，會使得斷層運(yùn)動加速或減速，進(jìn)而產(chǎn)生非周期的地震滑移．以成核階段為例(圖8)，由于斷層在該階段滿足vθ/Dc?1，所以可利用式(8)和式(13)推導(dǎo)出當(dāng)t0時刻滑移速率為v0時，斷層的失穩(wěn)時間(未經(jīng)應(yīng)力擾動)為(Dieterich，1992，1994)

圖8 施加和未施加正向靜態(tài)剪切應(yīng)力擾動時的滑移速率演化示意圖實線為未擾動時的演化，虛線為存在擾動時的演化，t0為施加靜態(tài)剪切應(yīng)力擾動的時刻，T和Tv分別為未擾動時和擾動后從t0時刻至斷層模型失穩(wěn)所需的時間

(34)

而在t0時刻受靜態(tài)剪切應(yīng)力擾動后，斷層的失穩(wěn)時間為

(35)

式中，Δτ為斷層受到的靜態(tài)剪切應(yīng)力擾動. 失穩(wěn)時間的提前或推后量為Δt=|T-Tv|．不難看出，只要施加在成核階段的靜態(tài)剪切應(yīng)力擾動非零，那么一定會導(dǎo)致失穩(wěn)時間的提前或推后，即成核階段的持續(xù)時間會發(fā)生變化，進(jìn)而產(chǎn)生非周期的地震滑移．而對于地震循環(huán)的其它階段，若施加靜態(tài)剪切應(yīng)力擾動也會產(chǎn)生同樣的現(xiàn)象，更為詳細(xì)的討論請參閱解孟雨和史保平(2016)文章．

5 非地震滑移

從靜態(tài)剪切應(yīng)力加載產(chǎn)生的非地震滑移圖(圖9)中可以看出，當(dāng)斷層隨時間的演化過程接近于自加速/成核階段的后期或震后滑移階段的早期，且斷層出現(xiàn)絕對值較大的剪切應(yīng)力加載時，斷層可以發(fā)生非地震滑移，從而也能打破斷層演化的周期性．圖9a給出了類似于脈沖的滑移速率擾動，其脈沖幅值遠(yuǎn)遠(yuǎn)小于同震滑移的速度值(0.14 m/s)，基本與遠(yuǎn)場加載速率vpl相當(dāng)．圖9b則分別給出了相應(yīng)的狀態(tài)變量隨時間的變化．在非地震滑移發(fā)生之后，可以再次出現(xiàn)新的同震滑移現(xiàn)象，但所需的時間間隔遠(yuǎn)小于未經(jīng)應(yīng)力擾動時的演化周期(圖9a)．

圖9c則以相圖的形式給出了斷層上剪切應(yīng)力與滑移速率之間的演化關(guān)系. 圖9c中出現(xiàn)的現(xiàn)象與Kame等(2013a)采用修正的RSF定律(Nagataetal, 2012)和一維彈簧-滑塊模型得到的數(shù)值模擬結(jié)果一致．Kaneko和Lapusta(2008)利用隨深度變化的摩擦本構(gòu)關(guān)系，也得到了非地震滑移現(xiàn)象．Barbot等(2012)在對帕克菲爾德地震重復(fù)發(fā)生過程的數(shù)值模擬中也出現(xiàn)了重復(fù)的非地震滑移．不過這些非地震滑移現(xiàn)象的出現(xiàn)有違于Dieterich(1992，1994)近似解的預(yù)測．

圖9 正向(左)和負(fù)向(右)靜態(tài)剪切應(yīng)力加載產(chǎn)生的非地震滑移(a) 滑移速率v演化圖；(b) 狀態(tài)變量θ演化圖；(c) 斷層面上剪切應(yīng)力τ和滑移速率v的相圖黑線和紅線分別為未受應(yīng)力擾動和受應(yīng)力擾動時一維彈簧-滑塊模型的演化情況，τss和vss分別為穩(wěn)定狀態(tài)下斷層面上的剪切應(yīng)力和滑移速率

6 討論與結(jié)論

本文在模擬計算中使用的是一維彈簧-滑塊模型，該模型相對簡單，雖然未考慮斷層面的形態(tài)、非均勻性和非彈性的影響，但卻能突出斷層的基本物理特征，而且計算簡單，有助于我們定性理解地震過程，進(jìn)而為更深入的研究提供方向．一維彈簧-滑塊模型等價于空間上滑移和應(yīng)力均勻分布的斷層模型(Rice，1983)，可以用來研究不同情況下斷層模型的一階近似情況 (Dieterich，1981)，能夠突出其基本物理特征，進(jìn)而為研究震源機(jī)制提供線索．例如，通過對一維彈簧-滑塊模型的分析認(rèn)識到，在描述該模型穩(wěn)定性時需引進(jìn)臨界彈簧剛度系數(shù)kc，并且只有當(dāng)彈簧的剛度系數(shù)滿足k

基于Dieterich-Ruina定律和一維彈簧-滑塊模型，本文利用解析模型與數(shù)值模擬相結(jié)合的原則探討了斷層演化過程中地震重復(fù)發(fā)生周期以及非地震滑移形成的可能力學(xué)機(jī)制．初步結(jié)論如下：

1) 數(shù)值模擬結(jié)果顯示斷層的整體演化至少包含4個主要階段，即閉鎖階段、自加速/成核階段、同震相和震后滑移階段，其間如果斷層受到外來應(yīng)力加載的擾動，斷層演化周期會被打斷，從而造成斷層失穩(wěn)時間的提前或推后．

2) 與庫侖模型所給出的時間提前量或推后量不同的是，當(dāng)斷層演化處于自加速/成核階段時，失穩(wěn)時間的提前量或推后量Δt的大小取決于應(yīng)力加載的時刻．當(dāng)正向應(yīng)力加載作用于震后滑移階段的早期或負(fù)向應(yīng)力加載作用于自加速階段的晚期(或同震滑移發(fā)生前)時，斷層演化過程中會出現(xiàn)非地震滑移現(xiàn)象，具體表現(xiàn)為脈沖形的滑移速率，其幅值與遠(yuǎn)場加載速率相當(dāng)．如今有關(guān)非地震滑移現(xiàn)象的討論已成為全球關(guān)于地震成因和地震預(yù)測研究的熱點(diǎn)(Perfettiniet al，2010;Noda，Lapusta，2013;Shirzaeiet al，2014;Avouac，2015;Uchidaet al，2016)，因此，目前的模擬計算結(jié)果可以為非地震滑移的成因提供新的認(rèn)識．

3) 基于對數(shù)值模擬結(jié)果的認(rèn)識，本文給出了估算地震重復(fù)發(fā)生周期(斷層演化循環(huán))的解析近似解，其中，斷層演化周期除了與斷層的特征尺度、遠(yuǎn)場加載速率和作用于斷層面上的正應(yīng)力相關(guān)外，同時也是Dieterich-Ruina定律中摩擦系數(shù)a，b以及臨界滑移距離Dc的函數(shù)．同前人工作的最主要區(qū)別在于，本文工作進(jìn)一步考慮了震后滑移階段和自加速/成核階段對斷層整體演化的貢獻(xiàn)．與數(shù)值模擬結(jié)果的比較也進(jìn)一步表明，目前的近似解對演化各個階段以及整個地震演化持續(xù)時間均給出了很好的近似．

需要強(qiáng)調(diào)的是，目前的工作對于斷層演化過程的認(rèn)識仍然是初步和粗淺的，真實斷層幾何狀態(tài)和斷層內(nèi)部物理/力學(xué)結(jié)構(gòu)的復(fù)雜性尚未考慮，因此，探索地震成因背后的摩擦力學(xué)機(jī)制仍需更為深入細(xì)致的研究工作，包括對理論、實驗、觀測以及方法等的進(jìn)一步認(rèn)識．

附錄 地震演化周期Tr的推導(dǎo)過程

應(yīng)力降Δτ滿足(Barbotetal，2012)

(A1)

另外，Barbot等(2012)還假定在應(yīng)力降發(fā)生前后斷層均為穩(wěn)定狀態(tài)且兩個穩(wěn)定狀態(tài)分別處于同震相和閉鎖階段，則有

(A2)

(A3)

式中，τco和τint分別為同震相和閉鎖階段的應(yīng)力值. 式(A3)減去式(A2)可得

(A4)

進(jìn)而有

(A5)

或

(A6)

解孟雨，史保平. 2016. 數(shù)值模擬靜態(tài)應(yīng)力擾動下的斷層失穩(wěn): 結(jié)果分析兼與Dieterich模型和Coulomb模型的對比[J]. 地球物理學(xué)報，59(2): 593--605.

Xie M Y，Shi B P. 2016. Numerical simulation of fault instability due to an arbitrary static stress perturbation: A comparison with the Dieterich model and Coulomb failure model[J].ChineseJournalofGeophysics，59(2): 593--605 (in Chinese).

仲秋，史保平. 2012. 1976年MS7.8唐山地震余震序列持續(xù)時間及對地震危險性分析的意義[J]. 地震學(xué)報，34(4)：494--508．

Zhong Q, Shi B P. 2012. Aftershock time duration of the 1976MS7.8 Tangshan earthquake and implication for seismic hazard[J].ActaSeismologicaSinica, 34(4): 494--508 (in Chinese).

Avouac J-P. 2015. From geodetic imaging of seismic and aseismic fault slip to dynamic modeling of the seismic cycle[J].AnnuRevEarthPlanetSci，43(1): 233--271.

Barbot S，Lapusta N，Avouac J-P. 2012. Under the hood of the earthquake machine: Toward predictive modeling of the seismic cycle[J].Science，336(6082): 707--710.

Beeler N M，Lockner D A. 2003. Why earthquakes correlate weakly with the solid earth tides: Effects of periodic stress on the rate and probability of earthquake occurrence[J].JGeophysRes，108(B8): 2391.

Beeler N M. 2004. Review of the physical basis of laboratory-derived relations for brittle failure and their implications for earthquake occurrence and earthquake nucleation[J].PureApplGeophys，161(9/10): 1853--1876.

Belardinelli M E，Bizzarri A，Cocco M. 2003. Earthquake triggering by static and dynamic stress changes[J].JGeophysRes，108(B3): 2135.

Bhattacharya P，Rubin A M. 2014. Frictional response to velocity steps and 1-D fault nucleation under a state evolution law with stressing-rate dependence[J].JGeophysRes，119(3): 2272--2304.

Byerlee J. 1970. The mechanics of stick-slip[J].Tectonophysics，9(5): 475--486.

Byerlee J. 1978. Friction of rocks[J].PureApplGeophys，116(4/5): 615--626.

Dieterich J H. 1978. Time-dependent friction and the mechanics of stick-slip[J].PureApplGeophys，116(4/5): 790--806.

Dieterich J H. 1979. Modelling of rock friction: 1. Experimental results and constitutive equations[J].JGeophysRes，84(B5): 2161--2168.

Dieterich J H. 1981. Constitutive properties of faults with simulated gouge[C]∥MechanicalBehaviorofCrustalRocks:TheHandinVolume. Washington DC: American Geophysical Union：103--120.

Dieterich J H. 1992. Earthquake nucleation on faults with rate- and state-dependent strength[J].Tectonophysics，211(1/2/3/4): 115--134.

Dieterich J H. 1994. A constitutive law for rate of earthquake production and its application to earthquake clustering[J].JGeophysRes，99(B2): 2601--2618.

Gomberg J，Blanpied M L，Beeler N M. 1997. Transient triggering of near and distant earthquakes[J].BullSeismolSocAm，87(2): 294--309.

Gomberg J，Beeler N M，Blanpied M L，Bodin P. 1998. Earthquake triggering by transient and static deformation[J].JGeophysRes，103(B10): 24411--24426.

Gu J C，Rice J R，Ruina A L，Tse S T. 1984. Slip motion and stability of a single degree of freedom elastic system with rate and state dependent friction[J].JMechPhysSolids，32(3): 167--196.

Harris R A，Simpson R W. 1998. Suppression of large earthquakes by stress shadows: A comparison of Coulomb and rate-and-state failure[J].JGeophysRes，1032(B10): 24439--24451.

Kame N，F(xiàn)ujita S，Nakatani M，Kusakabe T. 2013a. Effects of a revised rate- and state-dependent friction law on aftershock triggering model[J].Tectonophysics，600: 187--195.

Kame N，F(xiàn)ujita S，Nakatani M，Kusakabe T. 2013b. Earthquake cycle simulation with a revised rate- and state-dependent friction law[J].Tectonophysics，600 : 196--204.

Kaneko Y，Lapusta N. 2008. Variability of earthquake nucleation in continuum models of rate-and-state faults and implications for aftershock rates[J].JGeophysRes，113(B12): B12312.

Kaneko Y，Avouac J-P，Lapusta N. 2010. Towards inferring earthquake patterns from geodetic observations of interseismic coupling[J].NatGeosci，3(5): 363--369.

Liu M，Stein S，Wang H. 2012. 2000 years of migrating earthquakes in North China: How earthquakes in midcontinents differ from those at plate boundaries[J].Lithosphere，22(2): 33--44.

Marone C. 1998. Laboratory-derived friction laws and their application to seismic faulting[J].AnnuRevEarthPlanetSci，26(1): 643--696.

Nagata K，Nakatani M，Yoshida S. 2012. A revised rate- and state-dependent friction law obtained by constraining constitutive and evolution laws separately with laboratory data[J].JGeophysRes，117(B2): B02314.

Noda H，Lapusta N. 2013. Stable creeping fault segments can become destructive as a result of dynamic weakening[J].Nature，493(7433): 518--521.

Parsons T. 2005. A hypothesis for delayed dynamic earthquake triggering[J].GeophysResLett，32(4): L04302.

Perfettini H，Schmittbuhl J，Cochard A. 2003. Shear and normal load perturbations on a two-dimensional continuous fault: 1. Static triggering[J].JGeophysRes，108(B9): 2408.

Perfettini H，Avouac J-P. 2004a. Postseismic relaxation driven by brittle creep: A possible mechanism to reconcile geodetic measurements and the decay rate of aftershocks，application to the Chi-Chi earthquake，Taiwan[J].JGeophysRes，109(B2): B02304.

Perfettini H，Avouac J-P. 2004b. Stress transfer and strain rate variations during the seismic cycle[J].JGeophysRes，109(6): 117--132.

Perfettini H，Avouac J-P，Tavera H，Kositsky A，Nocquet J-M，Bondoux F，Chlieh M，Sladen A，Audin L，F(xiàn)arber D L，Soler P. 2010. Seismic and aseismic slip on the Central Peru megathrust[J].Nature，465(7294): 78--81.

Ranjith K，Rice J R. 1999. Stability of quasi-static slip in a single degree of freedom elastic system with rate and state dependent friction[J].JMechPhysSol，47(6): 1207--1218.

Rice J R. 1983. Constitutive relations for fault slip and earthquake instabilities[J].PureApplGeophys，121(3): 443--475.

Rubin A M，Ampuero J-P. 2005. Earthquake nucleation on (aging) rate and state faults[J].JGeophysRes，110(B11): B11312.

Ruina A. 1983. Slip instability and state variable friction laws[J].JGeophysRes，88(B12): 10359--10370.

Ryder I. 2006.ElasticandViscoelasticModellingofPostseismicMotionandFaultStructures[D]. Oxford: Oxford University: 12.

Scholz C H. 1998. Earthquakes and friction laws[J].Nature，391(6662): 37--42.

Scholz C H. 2002.TheMechanicsofEarthquakesandFaulting[M]. New York: Cambridge University Press: 85，266--267.

Segall P. 2010.EarthquakeandVolcanoDeformation[M]. Princeton: Princeton University Press: 332，338--339，349--350.

Shimazaki K，Nakata T. 1980. Time-predictable recurrence model for large earthquake[J].GeophysResLett，7(4): 279--282.

Shirzaei M，Burgmann R，Uchida N，Pollitz F，Matsuzawa T. 2014. Seismicversusaseismic slip: Probing mechanical properties of the northeast Japan subduction zone[J].EarthPlanetSciLett，406: 7--13.

Stein R S. 1999. The role of stress transfer in earthquake occurrence[J].Nature，402(6762): 605--609.

Stein S，Wysession M. 2003.AnIntroductiontoSeismology，Earthquakes，andEarthStructure[M]. Malden: Blackwell Publishing: 215--217.

Stein S，Liu M. 2009. Long aftershock sequences within continents and implications for earthquake hazard assessment[J].Nature，462(7269): 87--89.

Uchida N，Iinuma T，Nadeau R M，Burgmann R，Hino R. 2016. Periodic slow slip triggers megathrust zone earthquakes in northeastern Japan[J].Science，351(6272): 488--492.

The periodic and aperiodic slip during earthquake faulting and aseismic faulting slip: 1D fault model analysis

Xie Mengyu*Shi Baoping

(CollegeofEarthScience，UniversityofChineseAcademyofSciences，Beijing100049，China)

With a simple 1D spring-slide block model governed by Dieterich-Ruina law，this study derived an approximate analytical solution used to calculate the earthquake recurrence time and compared its approximate solution with previous results given by Barbotetal. The comparison between approximate solutions and the results by numerical simulation also demonstrated that current solution gives us a better approximation compared with previous result. Furthermore，it should be emphasized that，because only 40% to 80% of earthquake recurrence time belongs to interseismic stage，therefore，the duration of postseismic and nucleation/preseismic stage could not be ignored in the whole faulting process. Specifically，the time scale of earthquake recurrence in current earthquake recurrence model is not only related to the characteristic dimension of fault，the effective normal stress and remote loading rate，but also strongly depends on the values of the frictional parameters in the Dieterich-Ruina law and the critical distance obtained from laboratory. It is necessary to point out that，if all parameters we used represent a real situation in the earth，the approximate solution proposed in the present paper can give us an excellent estimation of earthquake recurrence time for a given tectonic region，of which relative error is less than 5%. Moreover，numerical simulation and theory are also used to study the mechanical mechanisms which cause the periodic and aperiodic evolution of the earthquake faulting and aseismic faulting slip. In addition，this paper also discussed the limitations and advantages of using 1D spring-slide block model to describe the elasticity in modeling evolution of fault. And it is found that applying shear stress loading during evolution of earthquake faulting can cause aperiodic evolution. Specially，applying shear stress loading at late nucleation/preseismic stage or early postseismic stage will lead to aseismic faulting slip.

rate- and state-dependent friction law; earthquake cycle; earthquake recurrence time; aseismic slip

10.11939/jass.2016.04.006.

國家國際科技合作專項(2015DFA21260)、北京市高校共建項目、國家基金委(41574040)和中國科學(xué)院創(chuàng)新團(tuán)隊項目(KZZD-EW-TZ-19)共同資助.

2016-03-15收到初稿，2016-05-27決定采用修改稿.

10.11939/jass.2016.04.006

P315.3

A

解孟雨, 史保平. 2016. 斷層演化過程中的周期和非周期地震滑移及非地震滑移----一維斷層模型分析. 地震學(xué)報, 38(4): 590--608.

Xie M Y, Shi B P. 2016. The periodic and aperiodic slip during earthquake faulting and aseismic faulting slip: 1D fault model analysis.ActaSeismologicaSinica, 38(4): 590--608.doi:10.11939/jass.2016.04.006.

*通訊作者 e-mail: jiemengyu13@mails.ucas.ac.cn