基于速率-狀態(tài)摩擦定律研究剪切應(yīng)力擾動對地震斷層演化的影響

高雅琪， 史保平

1 中國科學(xué)院大學(xué)地球與行星科學(xué)學(xué)院， 北京 100049 2 中國地震局第一監(jiān)測中心， 天津 300180

0 引言

地震成因機制研究是現(xiàn)代地震科學(xué)所面臨的一個重要挑戰(zhàn)(National Research Council，2008；Lay et al.，2009).近幾十年來對于觸發(fā)地震和誘發(fā)地震的研究也已成為地球物理學(xué)的一個重要領(lǐng)域(Dieterich，1994；萬永革等，2002；Ellsworth，2013；Brodsky and van der Elst，2014；Shapiro，2015; 解朝娣等，2021).因此，深入了解觸發(fā)地震的摩擦力學(xué)過程對于現(xiàn)代地震危險性評估和地震預(yù)測有著重要的應(yīng)用價值和科學(xué)意義.對于發(fā)震斷層的應(yīng)力擾動形式以及所對應(yīng)的斷層摩擦滑移響應(yīng)的理解也成為我們深入探討斷層時空演化機制一個有意義的課題，速率和狀態(tài)相關(guān)的摩擦定律(Dieterich，1992, 1994；Marone，1998)則成為當前定量了解地震觸發(fā)機制一個重要組成部分.

由于地震活動的復(fù)雜性，當?shù)貧?nèi)部出現(xiàn)不同時空尺度的應(yīng)力變化時，應(yīng)力躍變和依賴于時間的應(yīng)力對發(fā)震斷層的加載，處于不同演化時期的斷層滑移和滑移速率將會受到不同程度的擾動，進而可改變原來的斷層演化進程.影響區(qū)域地震活動性出現(xiàn)變化的應(yīng)力擾動主要可分為兩種：1)靜態(tài)應(yīng)力擾動；2)動態(tài)應(yīng)力擾動.針對靜態(tài)擾動模型，前人已建立了比較完善的物理模型并用于對區(qū)域地震活動性變化的闡述(Gross and Kisslinger，1997；Gross and Bürgmann，1998；Toda et al.，1998，2012；Dieterich et al.，2000；Perfettini et al.，2003a；孫楠等，2014).針對動態(tài)擾動模型，主要的研究則側(cè)重于周期/非周期應(yīng)力擾動對斷層失穩(wěn)的探討(Heki，2003；Perfettini et al.，2003b；Cochran et al.，2004；Bollinger et al.，2007；Christiansen et al.，2007；Bettinelli et al.，2008；Ader and Avouac，2013；van der Elst et al.，2013；孫楠等，2014; van der Elst and Savage，2015；Yoshida，2018；Yoshida et al.，2020).近期的研究則擴展到由地表及近地表應(yīng)力加載變化同區(qū)域地震活動相關(guān)性等方面的探討，例如，由水文循環(huán)(Heki，2003；Christiansen et al.，2005，2007；Bollinger et al.，2007；Bettinelli et al.，2008；Ellsworth，2013；Kundu et al.，2015)引起的季節(jié)性應(yīng)力變化同地震活動的關(guān)聯(lián)性.喜馬拉雅山地區(qū)季風(fēng)與溫度的變化會所引起當?shù)厮牡淖兓绾黾雍蜏p少相當于地表應(yīng)力的加載和卸載作用.Ader和Avouac(2013)則進一步分析了這樣的準周期性載荷同區(qū)域地震活動的關(guān)聯(lián)性.此外，人類活動如油頁巖開采過程中地下注水引發(fā)的地下孔隙壓力的變化，從而形成近地表應(yīng)力加載的準周期性變化，亦可誘發(fā)地震活動(Shapiro，2015).同樣，潮汐、固體潮等周期性載荷對地震觸發(fā)的可能影響也為大量學(xué)者所研究(Wilcock，2001；Tanaka et al.，2002；Beeler and Locker，2003；Cochran et al.，2004；馮向東和魏東平，2007；吳小平等，2009).

目前，應(yīng)用最廣泛的狀態(tài)演化定律分別是由Dieterich和Ruina提出的(Dieterich，1979a，b；Ruina，1983)，可以用于解釋巖石塊體接觸面摩擦過程的復(fù)雜行為.由于在實際應(yīng)用時所具有的優(yōu)缺點(見下文原理部分的闡述)，許多學(xué)者都對這兩個方程進行了必要的修正.但是，一致的狀態(tài)變量演化的數(shù)學(xué)形式仍然是一個尚未解決的問題.實際上，現(xiàn)有的狀態(tài)演化方程沒有一個能與實驗數(shù)據(jù)完全吻合(Beeler et al.，1994；Marone et al.，1995；Bayart et al.，2006；He and Wong，2014).Nagata等(2012)在狀態(tài)演化方程中引入了應(yīng)力速率項，給出了修正后的Dieterich演化定律.近期的研究結(jié)果表明，Nagata定律綜合了Dieterich和Ruina定律的優(yōu)點，使斷層模型在對正、負向擾動的響應(yīng)具有對稱性的同時，又能表達斷層愈合過程(Bhattacharya and Rubin，2014).Nagata演化定律已應(yīng)用于對斷層成核過程、斷層演化及余震發(fā)生率模型等多個領(lǐng)域的數(shù)值分析研究(Kame et al.，2013a,b，2015； Bhattacharya and Rubin，2014)，Yoshida等(Yoshida，2018；Yoshida et al., 2020)則采用該定律探討了高頻地震波作用下的地震動態(tài)觸發(fā)機制.本文采用1D彈簧-滑塊模型，從Nagata演化定律出發(fā)，進一步探討靜態(tài)應(yīng)力擾動和周期性應(yīng)力擾動，尤其極高頻應(yīng)力擾動對斷層演化過程的影響.我們將介紹當前研究的基本原理，展示靜態(tài)應(yīng)力擾動和周期性應(yīng)力擾動下的斷層演化過程的區(qū)別，分析應(yīng)力擾動量大小及擾動頻率等因素對斷層失穩(wěn)時刻、地震重復(fù)發(fā)生周期的影響.

1 原理

本文采用1D彈簧-滑塊模型來類比真實斷層，該模型由彈性的彈簧和剛性的滑塊組成.根據(jù)彈性理論給出的應(yīng)力-應(yīng)變關(guān)系，斷層內(nèi)部的剪切應(yīng)力可以表示為：

τd=τin+τload-kδ-V,

(1)

公式(1)中，τin為斷層內(nèi)部的初始應(yīng)力或剩余應(yīng)力；τload為外部加載于斷層面的剪切應(yīng)力，可為一依賴于時間的函數(shù)；δ和V分別是滑塊的滑動位移和速率(V=dδ/dt)；k=ηG/l為剛度系數(shù)，其中G為斷層面的剪切模量，l為斷層的特征尺度，η為描述斷層或裂紋幾何狀態(tài)的參數(shù)，與真實斷層面的幾何形態(tài)有關(guān)(Ruina，1983；Dieterich，1992，1994；Scholz，2002；Segall，2010)，其取值一般近于1(Dieterich，1992；Marone，1998)；=G/(2CS)為輻射衰減系數(shù)(Rice，1993)，其中CS表示剪切波速度.對于地殼內(nèi)部發(fā)震構(gòu)造而言，最簡單的設(shè)定是板塊驅(qū)動所形成的恒定剪切應(yīng)力速率加載，即外部剪切應(yīng)力加載僅由遠場板塊運動引起：

(2)

(3)

公式(3)中，Δτ0表征靜態(tài)擾動的方向(正向或負向)和大小，H(t)為階梯函數(shù).與靜態(tài)擾動相對應(yīng)的是動態(tài)擾動.進而，公式(3)可改寫為：

(4)

公式(4)中，f(t)表示一般應(yīng)力隨時間變化的函數(shù).在目前的工作中，我們僅考慮f(t)=Δτ0sin(ωt)下的周期性擾動的情形.其中，Δτ0表征簡諧擾動的振幅值，ω為角頻率.

斷層的摩擦強度τf可表示為：

τf=μσ,

(5)

公式(5)中σ為有效正應(yīng)力；μ為摩擦系數(shù)，可由RSF定律表達(Dieterich，1979a，b；Ruina，1983)：

(6)

公式(6)中，μ為摩擦系數(shù)，V為滑移速率，θ為狀態(tài)變量，V*是參考速率，μ0是當V=V*時的摩擦系數(shù).a是直接影響系數(shù)，決定由速度變化引起的摩擦變化；b是演化影響系數(shù)，決定由狀態(tài)變化引起的摩擦變化；DC是特征滑移距離.在斷層的演化過程中，加載剪切應(yīng)力τd與摩擦剪切強度τf須保持相等.

僅用公式(6)描述斷層摩擦是不完整的，它只說明了摩擦同速度和狀態(tài)的關(guān)系，還必須規(guī)定狀態(tài)變量如何隨時間或滑移速率變化，即所謂的狀態(tài)演化定律.目前，應(yīng)用最廣泛的兩個演化方程為(Dieterich，1979a，b；Ruina，1983)：

(7)

(8)

公式(7)稱為Dieterich演化定律；而公式(8)則稱為Ruina演化定律.它們能夠定量地描述斷層面上的摩擦性質(zhì)，從物理上闡明斷層內(nèi)部摩擦的復(fù)雜性.然而，二者都不能與巖石實驗結(jié)果完全一致(Beeler et al.，1994；Marone et al.，1995；Kato and Tullis，2001；Bayart et al.，2006；Nagata et al.，2012).

近期，基于實驗數(shù)據(jù)，Nagata等(2012)通過將實驗獲取的剪切應(yīng)力和滑動速度代入本構(gòu)關(guān)系的方法，計算了RSF定律中狀態(tài)變量的值，其結(jié)果與上述兩種演化定律的預(yù)測結(jié)果存在系統(tǒng)偏差.為此，他們提出了一個修正的演化定律，該定律考慮了剪切應(yīng)力變化對斷層摩擦的弱化作用，在Dieterich定律的基礎(chǔ)上，通過常數(shù)c引入了應(yīng)力速率項：

(9)

c為決定應(yīng)力速率項弱化作用大小的系數(shù)，為一大于零的常數(shù).將公式(6)和(9)結(jié)合，本文稱之為Nagata定律.顯然，當c=0時，Nagata定律退化為Dieterich定律；而當c趨于+∞時，Nagata定律則等效于Ruina定律(Bhattacharya and Rubin，2014).需要強調(diào)的是，在目前的研究中，Nagata定律中的摩擦參數(shù)a，b和DC的大小依賴于c的取值，其大小與Dieterich定律和Ruina定律的摩擦參數(shù)滿足如下比例關(guān)系：

(1+c)DC|N=DC|DR,

(10a)

(10b)

(a-b)|N=(a-b)|DR,

(10c)

上述三個等式左邊的a，b和DC為Nagata定律的摩擦參數(shù)，右邊的a，b和DC對應(yīng)Dieterich和Ruina定律.

線性穩(wěn)定性分析是1D彈簧-滑塊模型研究中的一個重要方面.結(jié)合Dieterich定律或Ruina定律，對1D彈簧-滑塊系統(tǒng)的穩(wěn)定性分析結(jié)果表明，彈簧的臨界剛度可以表示為(Ruina，1983；Segall，2010)：

(11)

(12)

2 結(jié)果

基于Nagata模型，在1D彈簧-滑塊模型的近似下，結(jié)合公式(1)(5)(6)和(9)，采用四階變步長的Runge-Kutta算法，我們實現(xiàn)了在無擾動、靜態(tài)擾動、周期性擾動三種情況下，對斷層摩擦過程的數(shù)值模擬，由此獲得不同應(yīng)力擾動對斷層演化進程的影響.為了便于對比模擬結(jié)果，本文采用Kato(2001)模型所設(shè)定的參數(shù)值，如表1所示.這組模型參數(shù)和初始條件可以模擬加州圣安德烈亞斯斷層大地震的震前滑動的狀態(tài).

表1 模型參數(shù)的含義和取值Table 1 Parameter definitions and values of model

2.1 靜態(tài)擾動下的斷層演化

假定t=0時刻施加應(yīng)力擾動Δτ0H(t)，擾動前系統(tǒng)處于穩(wěn)定滑移狀態(tài)，擾動瞬間應(yīng)力和速率發(fā)生突變，擾動后系統(tǒng)以恒定的速率演化至新的穩(wěn)態(tài).圖1展示了不同c值時Nagata模型的應(yīng)力Δτ隨滑動位移的演化過程，這里Δτ指的是擾動后某一時刻應(yīng)力相對于最終達到新的穩(wěn)定滑移狀態(tài)時應(yīng)力的水平.

由圖1可以看出，當c較小時，Nagata模型對正、負兩個方向擾動的響應(yīng)是完全不對稱的，對于相同大小的負向擾動，它能在更小的位移上達到穩(wěn)態(tài)水平，這同Dieterich模型給出的結(jié)果一致(Kato and Tullis，2001；Bhattacharya and Rubin，2014)；當c較大時，Nagata模型對正、負兩個方向擾動的響應(yīng)則趨于對稱，如圖1c和圖1d中c取10和100時所給出的圖像，這同Ruina模型給出的結(jié)果接近(Kato and Tullis，2001；Bhattacharya and Rubin，2014).實際上，當c=0時，Nagata模型退化為Dieterich模型；而當c趨于+∞時，Nagata模型等價于Ruina模型；隨著c由0增大到 +∞，Nagata模型對正、負擾動的響應(yīng)逐漸由不對稱向?qū)ΨQ轉(zhuǎn)變.Nagata模型的這一特征，本質(zhì)上源于其演化方程(9)中引入了應(yīng)力速率項，圖1給出的結(jié)果同Bhattacharya和Rubin(2014)的結(jié)果一致，也同Dieterich模型和Ruina模型相似(Kato and Tullis，2001).

圖1 靜態(tài)擾動下，Nagata模型的剪切應(yīng)力隨滑動位移的演化不同子圖對應(yīng)了不同c值的Nagata模型的剪切應(yīng)力演化，(a)(b)(c)和(d)分別對應(yīng)c = 0，2，10和100；不同顏色的曲線對應(yīng)不同方向(正向或負向)和大小的靜態(tài)應(yīng)力擾動.Fig.1 The evolution of shear stress with slip displacement by Nagata model under a static changeEach subplot shows shear stress evolution of Nagata model with a given c value. (a), (b), (c), and (d) correspond to c=0, 2, 10 and 100, respectively; and curves with different colors correspond to static stress perturbations with different sign (positive or negative) and magnitude.

上述結(jié)果也說明我們采用Nagata模型來模擬斷層演化的一個重要原因.Dieterich模型和Ruina模型不能同時解釋斷層摩擦的時間依賴性和速率依賴性：前者可以表達滑移速率為0時的斷層愈合過程，但是對正向擾動和負向擾動的響應(yīng)是不對稱的；后者雖然對正、負向擾動的響應(yīng)是對稱的，可它不允許滑移速率為0，即不能表達斷層愈合過程.然而，研究表明，Nagata模型則綜合了Dieterich模型和Ruina模型的優(yōu)點，在具有正負擾動響應(yīng)對稱性的同時又可以表達斷層的愈合過程.因此，本文選取了Nagata模型來模擬斷層演化過程，重點在于探討應(yīng)力擾動對斷層失穩(wěn)時間、地震周期等方面的影響.

(13)

公式(13)中H=b/DC-(1+c)k/σ，ta為一特征時間尺度，表示斷層成核過程的時間尺度(Dieterich，1994；Beeler and Locker，2003).在靜態(tài)擾動作用下，斷層失穩(wěn)時間ts則更新為(Kame et al.，2013b)：

(14)

從公式(13)和(14)可以看出，在靜態(tài)應(yīng)力擾動下，斷層失穩(wěn)時刻與擾動方向和大小有關(guān).對公式(14)進行近似分析，當負向應(yīng)力擾動Δτ0?-aσ/(1+c)時，可以簡化為：

(15)

從公式(15)可以看出，較大的負向靜態(tài)擾動下，失穩(wěn)時刻ts與擾動大小Δτ0近似成正比.而當正向擾動Δτ0?aσ/(1+c)時，我們可將公式(14)簡化為：

(16)

圖2 靜態(tài)擾動下，Nagata模型的滑移速率和位移隨時間的演化不同子圖對應(yīng)了不同c值的Nagata模型的滑移速率或位移演化，(a)(b)(c)和(d)分別對應(yīng)c = 0，2，10和100；不同顏色的曲線對應(yīng)不同方向(正向或負向)和大小的靜態(tài)應(yīng)力擾動.Fig.2 The evolution of slip rate and displacement with time by Nagata model under a static changeEach subplot shows slip rate or displacement evolution of Nagata model with a given c value. (a), (b), (c), and (d) correspond to c=0, 2, 10 and 100, respectively; and curves with different colors correspond to static stress perturbations with different sign (positive or negative) and magnitude.

從公式(16)可以看出，較大的正向靜態(tài)擾動下，失穩(wěn)時刻ts與擾動大小Δτ0的指數(shù)成正比.實際上，Perfettini等(2003a)也得到了靜態(tài)應(yīng)力擾動下失穩(wěn)時刻的類似結(jié)果.

綜合上述關(guān)于靜態(tài)應(yīng)力擾動的研究，可以得到結(jié)論：在負向擾動下(Δτ0<0)，斷層失穩(wěn)比無擾動情形更推遲發(fā)生，且擾動量越大則推后量越大，失穩(wěn)時刻ts與擾動大小Δτ0成正比；在正向擾動下(Δτ0>0)，斷層失穩(wěn)比無擾動情形更提早發(fā)生，失穩(wěn)時刻ts與擾動大小Δτ0的指數(shù)成正比.總之，負向應(yīng)力擾動在一定程度上能減緩斷層演化進程，而正向應(yīng)力擾動則能促進斷層失穩(wěn).

2.2 周期性擾動下的斷層演化

在本文中，周期性應(yīng)力擾動函數(shù)f(t)=Δτ0sin(ωt)，因此擾動量的變化由幅值Δτ0和角頻率ω所決定，而ω與擾動周期Tr相關(guān)(Tr=2π/ω).考慮到季節(jié)性氣候(例如河流或湖水上漲)造成的應(yīng)力變化周期約為一年，那么不妨取Tr=1 a，則Nagata模型在不同幅值的周期性擾動下，速率和位移隨時間的演化過程如圖3所示.其中，擾動幅值Δτ0分別取0.1 MPa，0.5 MPa，1 MPa，2 MPa.

由圖3可以看出，與無擾動情形(Δτ0=0)相比，周期性應(yīng)力擾動下的斷層失穩(wěn)總是提前發(fā)生，且擾動幅值越大，失穩(wěn)時間提前量越大.造成此現(xiàn)象的原因，可以從累積滑動位移角度得到了解.從圖3可以看出，雖然滑移速率隨時間的增加存在一定波動(a1，b1，c1，d1)，但是累積滑動位移隨時間的增加是單調(diào)遞增的(a2，b2，c2，d2).當應(yīng)力擾動量大于0時，滑移速率相應(yīng)增大，累積滑動位移增加量也較大；而當應(yīng)力擾動量小于0時，滑移速率是下降的，從而使得累積滑動位移增加量上升變緩.事實上，擾動加載的持續(xù)時間為Tc，累積滑動位移量約為Δu≈V0Tc·I0(X)，I0(X)為零階修正的Bessel函數(shù)(I0(X)≥1，X=(1+c)Δτ0/(aσ)).因此，sin(ωt)形式下的周期性擾動促進了斷層滑移，故失穩(wěn)時間總是提前.

除了季節(jié)性氣候變化的因素，自然界還存在更多種頻率的應(yīng)力擾動，例如地震應(yīng)力波和人類活動的干擾(地?zé)衢_采，油頁巖開采等).為此，我們進一步探討了擾動頻率對斷層演化進程的影響.圖4展示了Nagata模型在不同頻率的周期性擾動下，速率和位移隨時間的演化過程.其中，固定擾動幅值Δτ0=1 MPa，調(diào)整擾動周期Tr=0.5 a，1 a，2 a.觀察圖像可以看出，對于相同的擾動幅值，不同擾動頻率下，斷層在相近的時刻發(fā)生失穩(wěn).這一現(xiàn)象說明，周期性擾動下的斷層失穩(wěn)時刻主要受擾動幅值影響，而擾動頻率對失穩(wěn)時間的影響很小.

圖4中的現(xiàn)象和結(jié)論也可以從近似解析解的角度進行分析和解釋.類似于Perfettini等(2003b)和Ader等(2014)得到的結(jié)果，當擾動周期Tr?ta時，斷層的失穩(wěn)時間td可近似為(詳見后文附錄A)：

(17)

當X>0時，I0(X)>1.對比公式(13)和(17)可以看出，無論Δτ0取值如何，總有td

圖4 不同頻率的周期性擾動下，Nagata模型的滑移速率和位移隨時間的演化不同子圖對應(yīng)了不同c值的Nagata模型的滑移速率或位移演化，(a)(b)(c)和(d)分別對應(yīng)c=0，2，10和100；不同顏色的曲線對應(yīng)不同頻率的周期性應(yīng)力擾動.Fig.4 The evolution of slip rate and displacement with time by Nagata model under different periodic disturbance frequencyEach subplot shows slip rate or displacement evolution of Nagata model with a given c value. (a), (b), (c), and (d) correspond to c=0, 2, 10 and 100, respectively; and curves with different colors correspond to periodic stress perturbations with different frequency.

綜合上述關(guān)于周期性應(yīng)力擾動的研究，可以得到結(jié)論：對于簡諧函數(shù)所描述的周期性應(yīng)力擾動，無論擾動量的大小、頻率和持續(xù)時間如何，周期性觸發(fā)下的斷層失穩(wěn)總是提前發(fā)生，且擾動幅值越大，則失穩(wěn)提前量越大.進一步講，如果地震應(yīng)力波都可以由傅里葉級函數(shù)表達，那么地震造成的周期性應(yīng)力擾動總是能導(dǎo)致斷層失穩(wěn)時間的提前，且提前量與該地震的震級有關(guān).

2.3 周期性擾動對地震復(fù)發(fā)周期的影響

前面兩部分研究的都是斷層的單次失穩(wěn)，若加入輻射衰減項，即公式(1)右邊的最后一項，我們可以實現(xiàn)斷層在時間域上演化的多次循環(huán)解，從而了解周期性應(yīng)力擾動對地震復(fù)發(fā)周期的影響.例如，采用應(yīng)力擾動函數(shù)f(t)= Δτ0sin(2πt/Tr)，其中，Δτ0=1 MPa，擾動周期Tr=1 a，可得如圖5、圖6和表2所示的模擬結(jié)果.圖5給出了演化過程中的斷層滑移速率和累積滑動位移，圖6和表2則分別給出了前21次和前11次失穩(wěn)中，每相鄰兩次失穩(wěn)之間的時間間隔，圖6橫坐標和表2第一行的數(shù)字表示失穩(wěn)間隔的序號.從圖5可以看出，當c值較小時(c=0，2)，無應(yīng)力擾動加入(紅色曲線)的斷層演化呈現(xiàn)出確定的周期性，如c=0時的循環(huán)周期約為41.2 a，c=2時的循環(huán)周期為38.5 a(圖6和表2).當c?1時，斷層循環(huán)的重復(fù)發(fā)生時間間隔會出現(xiàn)極小的波動，但仍然保持了固有的周期性，如c=10時的循環(huán)周期約為37.8 a，c=100時的循環(huán)周期約為40.5 a(圖6和表2).圖5中黑色實線則分別給出了Δτ0=1 MPa時Nagata模型多次循環(huán)下的滑移速率和滑動位移的變化.從滑移速率隨時間變化的特征來看，應(yīng)力擾動主要影響斷層震間階段的演化，會出現(xiàn)高頻波動的現(xiàn)象，波動周期與擾動周期Tr一致，均遠小于震間階段的持續(xù)時間.顯然，從圖5、圖6和表2可以看出，斷層失穩(wěn)的時間間隔出現(xiàn)了大于或小于固有周期(無應(yīng)力擾動情形的周期)的現(xiàn)象，斷層演化不再具有明確的周期性.當擾動量Δτ0增加時，斷層失穩(wěn)的時間間隔會變得更不規(guī)則，而這種不規(guī)則性在c=0時較為明顯.

圖5 周期性擾動下，Nagata模型的滑移速率和位移關(guān)于時間的周期解不同子圖對應(yīng)了不同c值的Nagata模型的滑移速率或位移演化，(a)(b)(c)和(d)分別對應(yīng)c=0，2，10和100；黑色曲線為周期性擾動情形，紅色曲線為無擾動情形.Fig.5 Periodic solutions of slip rate and displacement with time by Nagata model under periodic disturbanceEach subplot shows slip rate or displacement evolution of Nagata model with a given c value. (a), (b), (c), and (d) correspond to c=0, 2, 10 and 100, respectively; and the black curve corresponds to the periodic disturbance case, and the red curve corresponds to the undisturbed case.

圖6 周期性動態(tài)擾動下，Nagata模型的失穩(wěn)時間間隔不同子圖對應(yīng)了不同c值的Nagata模型的失穩(wěn)時間間隔，(a)(b)(c)和(d)分別對應(yīng)c=0，2，10和100；不同顏色的曲線對應(yīng)不同幅值的周期性應(yīng)力擾動.Fig.6 Instability time interval of Nagata model under periodic dynamic disturbanceEach subplot shows instability time interval of Nagata model with a given c value. (a), (b), (c), and (d) correspond to c=0, 2, 10 and 100, respectively; and curves with different colors correspond to periodic stress perturbations with different magnitude.

表2 無擾動和周期性擾動情形下，Nagata模型的失穩(wěn)時間間隔對比(單位：a)Table 2 Comparison of instability time intervals by Nagata model in the case of no disturbance and periodic disturbance (unit: a)

3 討論

本文通過數(shù)值方法，實現(xiàn)了靜態(tài)和周期性應(yīng)力擾動下，不同c值條件下Nagata模型對斷層演化過程的模擬.研究表明，正向靜態(tài)擾動和所有周期性擾動都可以使失穩(wěn)時間提前(圖2，圖3)，但是相同Δτ0條件下，兩種擾動造成的時間提前量不同.當Δτ0分別取0.1 MPa，0.5 MPa，1 MPa，2 MPa時，與無擾動(Δτ0=0 MPa)相比，Nagata模型在靜態(tài)擾動和周期性擾動下的失穩(wěn)時刻如表3所示(結(jié)果保留小數(shù)點后兩位).

表3數(shù)據(jù)顯示，對于相同的Δτ0，與周期性擾動相比，靜態(tài)擾動下的斷層失穩(wěn)更早、時間提前量更大.這也揭示了周期性擾動和靜態(tài)擾動的不同之處.公式(14)和(17)可以從模型原理角度解釋對于相同Δτ0，靜態(tài)擾動比周期性擾動下的斷層失穩(wěn)更早發(fā)生的原因.雖然周期性擾動下的時間提前量小于靜態(tài)擾動的結(jié)果，但周期性擾動對失穩(wěn)時刻的影響還是十分明顯，尤其當Δτ0很大時.

表3 靜態(tài)和周期性擾動下，Nagata模型的失穩(wěn)時刻(單位：a)Table 3 Instability time by Nagata model in the case of static and periodic disturbance (unit: a)

另外，從表3還可以發(fā)現(xiàn)，在相同的擾動條件下，c值越大，斷層失穩(wěn)所需時間越短.實際上，當c=0時，Nagata模型的失穩(wěn)時刻最晚，同Dieterich模型一致；隨著c的增大，Nagata模型失穩(wěn)越來越早；當c趨于+∞時，Nagata模型的失穩(wěn)時刻最早，接近于Ruina模型.

采用Nagata模型實現(xiàn)本文給出的數(shù)值模擬結(jié)果，其模型參數(shù)的取值也十分重要.對于目前的研究，Nagata模型的摩擦參數(shù)a，b和DC的取值與c值大小有關(guān)，在給定a，b和DC原始值的基礎(chǔ)上，按照公式(10)給出的比例關(guān)系取值.公式(10)并非隨意得到的，而是通過嚴格的物理與數(shù)學(xué)推導(dǎo)所得到(Bhattacharya and Rubin，2014)，但由于文章的篇幅有限，在此不再展開討論.

即使Nagata模型能兼顧Dieterich模型和Ruina模型的優(yōu)點，同時解決擾動響應(yīng)對稱性的問題和斷層愈合過程的問題，但是常參數(shù)c的取值仍然是一個尚未解決的問題.雖然應(yīng)力速率項的引入使Nagata模型在c較小時具有Dieterich模型的特征，在c較大時具有Ruina模型的特征，但是也很難找到一個合適的c值，將兩種模型的特征同時融合在某一個Nagata模型中.若c值過小，則Nagata模型對正、負向擾動的響應(yīng)不對稱；若c值過大，則按公式(10)給出的摩擦參數(shù)不在合理范圍內(nèi).前人最常用的c值為2，也有人認為c的上限為10(Bhattacharya et al., 2015)，但目前沒人能給出c的最佳值.

4 結(jié)論

本文從Nagata定律出發(fā)，結(jié)合1D彈簧-滑塊模型，通過四階變步長的Runge-Kutta算法，實現(xiàn)了對靜態(tài)和周期性應(yīng)力擾動下斷層演化過程的數(shù)值模擬.研究的重點集中在Nagata模型在不同擾動下的失穩(wěn)時刻和地震周期，探討應(yīng)力擾動對斷層失穩(wěn)機制的影響.主要發(fā)現(xiàn)和結(jié)論如下：

(1)當c=0時，Nagata模型退化為Dieterich模型，對正負靜態(tài)擾動的響應(yīng)是不對稱的；當c趨于+∞時，Nagata模型等價于Ruina模型，對正負靜態(tài)擾動的響應(yīng)是對稱的；隨著c的增大，Nagata模型下的靜態(tài)應(yīng)力擾動響應(yīng)由不對稱逐漸趨于對稱.

(2)在靜態(tài)擾動下，斷層失穩(wěn)時刻與擾動的方向和大小有關(guān)：與無擾動情形相比，對于負向擾動，失穩(wěn)總是推遲發(fā)生，且時間推后量與擾動量大小成正比；對于正向擾動，失穩(wěn)總是提早發(fā)生，且時間提前量與擾動量的指數(shù)成正比.

(3)簡諧應(yīng)力擾動總是使斷層失穩(wěn)提早發(fā)生，且擾動振幅越大，失穩(wěn)提前量越大；但是，擾動頻率或擾動周期對失穩(wěn)時刻并無顯著影響.對于相同的Δτ0，正向靜態(tài)擾動下的斷層失穩(wěn)更早，時間提前量比周期性擾動時的更大.

(4)當遠場加載速率給定且為常量時，若無應(yīng)力擾動，則斷層的失穩(wěn)時間間隔基本相同，斷層演化過程顯示出確定的周期性；但在周期性應(yīng)力的擾動下(簡諧應(yīng)力變化)，斷層失穩(wěn)時間間隔則會出現(xiàn)明顯的非周期性，從而使得斷層在時間域上的演化變得更為復(fù)雜.

附錄A 不同擾動下的斷層失穩(wěn)時刻

在斷層演化過程中，摩擦剪應(yīng)力τf始終等于加載剪應(yīng)力τd，結(jié)合正文原理部分的公式(1)(5)(6)，關(guān)于時間求導(dǎo)，整理可得：

(A1)

其中，H=b/DC-(1+c)k/σ.公式(A1)兩邊同時做不定積分，并用初始條件確定積分常數(shù)，經(jīng)過一系列整理和簡化可得：

(A2)

對于不同的情形，τload的具體形式不同，下面就分別從無擾動、靜態(tài)擾動、周期性擾動三種情形下出發(fā)，分析推導(dǎo)斷層失穩(wěn)時間tf、ts、td.

無擾動情形

對于無擾動情形，外界加載應(yīng)力τload可由正文的公式(2)表達，將其代入公式(A2)，積分后整理可得：

(A3)

tf就是無擾動情形下的斷層失穩(wěn)時間.

靜態(tài)擾動情形

對于靜態(tài)擾動情形，外界加載應(yīng)力τload可由正文的公式(3)表達，將其代入公式(A2)，積分后整理可得：

(A4)

ts就是靜態(tài)擾動情形下的斷層失穩(wěn)時間.

周期性擾動情形

對于周期性擾動情形，外界加載應(yīng)力τload可由正文的公式(4)表達，將其代入公式(A2)：

(A5)

以2π/ω為間隔，將積分總區(qū)間[0,td]劃分為n個子區(qū)間(n=ωtd/2π)，若擾動頻率足夠大，滿足ω?2π/td，則公式(A5)可以簡記為：

(A6)

(A7)

g(Δτ0)表示公式(A5)第二部分周期函數(shù)在一個周期上的積分值：

(A8)

令Y=ωy，利用換元法，再利用正余弦函數(shù)的周期性和對稱性，可將公式(A8)改寫為：

(A9)

(A10)

td就是周期性擾動情形下的斷層失穩(wěn)時間.