高雅琪， 史保平

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

0 引言

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

由于地震活動(dòng)的復(fù)雜性，當(dāng)?shù)貧?nèi)部出現(xiàn)不同時(shí)空尺度的應(yīng)力變化時(shí)，應(yīng)力躍變和依賴于時(shí)間的應(yīng)力對(duì)發(fā)震斷層的加載，處于不同演化時(shí)期的斷層滑移和滑移速率將會(huì)受到不同程度的擾動(dòng)，進(jìn)而可改變?cè)瓉?lái)的斷層演化進(jìn)程.影響區(qū)域地震活動(dòng)性出現(xiàn)變化的應(yīng)力擾動(dòng)主要可分為兩種：1)靜態(tài)應(yīng)力擾動(dòng)；2)動(dòng)態(tài)應(yīng)力擾動(dòng).針對(duì)靜態(tài)擾動(dòng)模型，前人已建立了比較完善的物理模型并用于對(duì)區(qū)域地震活動(dòng)性變化的闡述(Gross and Kisslinger，1997；Gross and Bürgmann，1998；Toda et al.，1998，2012；Dieterich et al.，2000；Perfettini et al.，2003a；孫楠等，2014).針對(duì)動(dòng)態(tài)擾動(dòng)模型，主要的研究則側(cè)重于周期/非周期應(yīng)力擾動(dòng)對(duì)斷層失穩(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).近期的研究則擴(kuò)展到由地表及近地表應(yīng)力加載變化同區(qū)域地震活動(dòng)相關(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)力變化同地震活動(dòng)的關(guān)聯(lián)性.喜馬拉雅山地區(qū)季風(fēng)與溫度的變化會(huì)所引起當(dāng)?shù)厮牡淖兓?，如湖水增加和減少相當(dāng)于地表應(yīng)力的加載和卸載作用.Ader和Avouac(2013)則進(jìn)一步分析了這樣的準(zhǔn)周期性載荷同區(qū)域地震活動(dòng)的關(guān)聯(lián)性.此外，人類(lèi)活動(dòng)如油頁(yè)巖開(kāi)采過(guò)程中地下注水引發(fā)的地下孔隙壓力的變化，從而形成近地表應(yīng)力加載的準(zhǔn)周期性變化，亦可誘發(fā)地震活動(dòng)(Shapiro，2015).同樣，潮汐、固體潮等周期性載荷對(duì)地震觸發(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)，可以用于解釋巖石塊體接觸面摩擦過(guò)程的復(fù)雜行為.由于在實(shí)際應(yīng)用時(shí)所具有的優(yōu)缺點(diǎn)(見(jiàn)下文原理部分的闡述)，許多學(xué)者都對(duì)這兩個(gè)方程進(jìn)行了必要的修正.但是，一致的狀態(tài)變量演化的數(shù)學(xué)形式仍然是一個(gè)尚未解決的問(wèn)題.實(shí)際上，現(xiàn)有的狀態(tài)演化方程沒(méi)有一個(gè)能與實(shí)驗(yàn)數(shù)據(jù)完全吻合(Beeler et al.，1994；Marone et al.，1995；Bayart et al.，2006；He and Wong，2014).Nagata等(2012)在狀態(tài)演化方程中引入了應(yīng)力速率項(xiàng)，給出了修正后的Dieterich演化定律.近期的研究結(jié)果表明，Nagata定律綜合了Dieterich和Ruina定律的優(yōu)點(diǎn)，使斷層模型在對(duì)正、負(fù)向擾動(dòng)的響應(yīng)具有對(duì)稱(chēng)性的同時(shí)，又能表達(dá)斷層愈合過(guò)程(Bhattacharya and Rubin，2014).Nagata演化定律已應(yīng)用于對(duì)斷層成核過(guò)程、斷層演化及余震發(fā)生率模型等多個(gè)領(lǐng)域的數(shù)值分析研究(Kame et al.，2013a,b，2015； Bhattacharya and Rubin，2014)，Yoshida等(Yoshida，2018；Yoshida et al., 2020)則采用該定律探討了高頻地震波作用下的地震動(dòng)態(tài)觸發(fā)機(jī)制.本文采用1D彈簧-滑塊模型，從Nagata演化定律出發(fā)，進(jìn)一步探討靜態(tài)應(yīng)力擾動(dòng)和周期性應(yīng)力擾動(dòng)，尤其極高頻應(yīng)力擾動(dòng)對(duì)斷層演化過(guò)程的影響.我們將介紹當(dāng)前研究的基本原理，展示靜態(tài)應(yīng)力擾動(dòng)和周期性應(yīng)力擾動(dòng)下的斷層演化過(guò)程的區(qū)別，分析應(yīng)力擾動(dòng)量大小及擾動(dòng)頻率等因素對(duì)斷層失穩(wěn)時(shí)刻、地震重復(fù)發(fā)生周期的影響.

1 原理

本文采用1D彈簧-滑塊模型來(lái)類(lèi)比真實(shí)斷層，該模型由彈性的彈簧和剛性的滑塊組成.根據(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í)間的函數(shù)；δ和V分別是滑塊的滑動(dòng)位移和速率(V=dδ/dt)；k=ηG/l為剛度系數(shù)，其中G為斷層面的剪切模量，l為斷層的特征尺度，η為描述斷層或裂紋幾何狀態(tài)的參數(shù)，與真實(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表示剪切波速度.對(duì)于地殼內(nèi)部發(fā)震構(gòu)造而言，最簡(jiǎn)單的設(shè)定是板塊驅(qū)動(dòng)所形成的恒定剪切應(yīng)力速率加載，即外部剪切應(yīng)力加載僅由遠(yuǎn)場(chǎng)板塊運(yùn)動(dòng)引起：

(2)

(3)

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

(4)

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

斷層的摩擦強(qiáng)度τf可表示為：

τf=μσ,

(5)

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

(6)

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

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

(7)

(8)

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

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

(9)

c為決定應(yīng)力速率項(xiàng)弱化作用大小的系數(shù)，為一大于零的常數(shù).將公式(6)和(9)結(jié)合，本文稱(chēng)之為Nagata定律.顯然，當(dāng)c=0時(shí)，Nagata定律退化為Dieterich定律；而當(dāng)c趨于+∞時(shí)，Nagata定律則等效于Ruina定律(Bhattacharya and Rubin，2014).需要強(qiáng)調(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)

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

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

(11)

(12)

2 結(jié)果

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

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

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

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

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

圖1 靜態(tài)擾動(dòng)下，Nagata模型的剪切應(yīng)力隨滑動(dòng)位移的演化不同子圖對(duì)應(yīng)了不同c值的Nagata模型的剪切應(yīng)力演化，(a)(b)(c)和(d)分別對(duì)應(yīng)c = 0，2，10和100；不同顏色的曲線對(duì)應(yīng)不同方向(正向或負(fù)向)和大小的靜態(tài)應(yīng)力擾動(dò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é)果也說(shuō)明我們采用Nagata模型來(lái)模擬斷層演化的一個(gè)重要原因.Dieterich模型和Ruina模型不能同時(shí)解釋斷層摩擦的時(shí)間依賴性和速率依賴性：前者可以表達(dá)滑移速率為0時(shí)的斷層愈合過(guò)程，但是對(duì)正向擾動(dòng)和負(fù)向擾動(dòng)的響應(yīng)是不對(duì)稱(chēng)的；后者雖然對(duì)正、負(fù)向擾動(dòng)的響應(yīng)是對(duì)稱(chēng)的，可它不允許滑移速率為0，即不能表達(dá)斷層愈合過(guò)程.然而，研究表明，Nagata模型則綜合了Dieterich模型和Ruina模型的優(yōu)點(diǎn)，在具有正負(fù)擾動(dòng)響應(yīng)對(duì)稱(chēng)性的同時(shí)又可以表達(dá)斷層的愈合過(guò)程.因此，本文選取了Nagata模型來(lái)模擬斷層演化過(guò)程，重點(diǎn)在于探討應(yīng)力擾動(dòng)對(duì)斷層失穩(wěn)時(shí)間、地震周期等方面的影響.

(13)

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

(14)

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

(15)

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

(16)

圖2 靜態(tài)擾動(dòng)下，Nagata模型的滑移速率和位移隨時(shí)間的演化不同子圖對(duì)應(yīng)了不同c值的Nagata模型的滑移速率或位移演化，(a)(b)(c)和(d)分別對(duì)應(yīng)c = 0，2，10和100；不同顏色的曲線對(duì)應(yīng)不同方向(正向或負(fù)向)和大小的靜態(tài)應(yīng)力擾動(dò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)擾動(dòng)下，失穩(wěn)時(shí)刻ts與擾動(dòng)大小Δτ0的指數(shù)成正比.實(shí)際上，Perfettini等(2003a)也得到了靜態(tài)應(yīng)力擾動(dòng)下失穩(wěn)時(shí)刻的類(lèi)似結(jié)果.

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

2.2 周期性擾動(dòng)下的斷層演化

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

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

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

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

(17)

當(dāng)X>0時(shí)，I0(X)>1.對(duì)比公式(13)和(17)可以看出，無(wú)論Δτ0取值如何，總有td

圖4 不同頻率的周期性擾動(dòng)下，Nagata模型的滑移速率和位移隨時(shí)間的演化不同子圖對(duì)應(yīng)了不同c值的Nagata模型的滑移速率或位移演化，(a)(b)(c)和(d)分別對(duì)應(yīng)c=0，2，10和100；不同顏色的曲線對(duì)應(yīng)不同頻率的周期性應(yīng)力擾動(dò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)力擾動(dòng)的研究，可以得到結(jié)論：對(duì)于簡(jiǎn)諧函數(shù)所描述的周期性應(yīng)力擾動(dòng)，無(wú)論擾動(dòng)量的大小、頻率和持續(xù)時(shí)間如何，周期性觸發(fā)下的斷層失穩(wěn)總是提前發(fā)生，且擾動(dòng)幅值越大，則失穩(wěn)提前量越大.進(jìn)一步講，如果地震應(yīng)力波都可以由傅里葉級(jí)函數(shù)表達(dá)，那么地震造成的周期性應(yīng)力擾動(dòng)總是能導(dǎo)致斷層失穩(wěn)時(shí)間的提前，且提前量與該地震的震級(jí)有關(guān).

2.3 周期性擾動(dòng)對(duì)地震復(fù)發(fā)周期的影響

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

圖5 周期性擾動(dòng)下，Nagata模型的滑移速率和位移關(guān)于時(shí)間的周期解不同子圖對(duì)應(yīng)了不同c值的Nagata模型的滑移速率或位移演化，(a)(b)(c)和(d)分別對(duì)應(yīng)c=0，2，10和100；黑色曲線為周期性擾動(dòng)情形，紅色曲線為無(wú)擾動(dòng)情形.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 周期性動(dòng)態(tài)擾動(dòng)下，Nagata模型的失穩(wěn)時(shí)間間隔不同子圖對(duì)應(yīng)了不同c值的Nagata模型的失穩(wěn)時(shí)間間隔，(a)(b)(c)和(d)分別對(duì)應(yīng)c=0，2，10和100；不同顏色的曲線對(duì)應(yīng)不同幅值的周期性應(yīng)力擾動(dò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 無(wú)擾動(dòng)和周期性擾動(dòng)情形下，Nagata模型的失穩(wěn)時(shí)間間隔對(duì)比(單位：a)Table 2 Comparison of instability time intervals by Nagata model in the case of no disturbance and periodic disturbance (unit: a)

3 討論

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

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

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

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

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

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

4 結(jié)論

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

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

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

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

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

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

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

(A1)

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

(A2)

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

無(wú)擾動(dòng)情形

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

(A3)

tf就是無(wú)擾動(dòng)情形下的斷層失穩(wěn)時(shí)間.

靜態(tài)擾動(dòng)情形

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

(A4)

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

周期性擾動(dòng)情形

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

(A5)

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

(A6)

(A7)

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

(A8)

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

(A9)

(A10)

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