汶川地震失穩機制的數值模擬研究

李平恩，廖力，奉建州

1 中國地震局地球物理研究所，北京 100081 2 北京白家疃國家地球科學野外科學觀測研究站，北京 100095

0 引言

強烈地震是一種嚴重的自然災害，給生命財產造成巨大損失.地震學家很早就認識到，大陸構造地震是斷層和圍巖系統在構造作用下的一種力學失穩過程(Nur, 1978; Rice, 1979; Stuart, 1979a)，其中，斷層的軟化特性是造成地震失穩的最重要原因(殷有泉和張宏，1984; Reches and Lockner, 2010).研究斷層地震的孕育發生過程必須考慮包含斷層和圍巖在內的整個巖石力學系統.Stuart(1979a)首先提出用準靜態力學方法研究地震不穩定性，并通過建立一個考慮應變軟化不穩定性的平面理論模型研究了1971年美國加州San Fernando 6.4級逆沖地震(Stuart, 1979b)，計算得到的震前和震后的地表抬升量與觀測資料吻合較好.在研究中，斷層軟化本構規律可采用在峰值應力后的應變軟化或應變率軟化形式，以及應力-滑動或應力-滑動率的摩擦規律來表示(Stuart and Mavko, 1979; Rice, 1983).考慮到Stuart(1979b)提出的應變軟化本構方程與空間坐標有關，不符合本構性質的客觀性原理，在理論上有瑕疵.殷有泉和張宏(1982，1984)建立了嚴謹的應變軟化的斷層面本構關系，用間斷面單元模擬斷層，在此基礎上給出了地震非穩定模型的一般數學表述，并從理論上指出地震失穩發生在峰值應力之后.當斷層本構關系的軟化行為采用理想化的速率無關的滑移弱化或速率狀態相關的摩擦關系進行描述時，地震失穩的位移應力曲線(Rudnicki,1988)形態與殷有泉和張宏(1984)的結果完全一致.這從理論上進一步肯定了斷層本構的軟化特性是導致地震失穩的必要條件，它與本構關系的具體描述形式無關.最近，李平恩和殷有泉(2009，2014)采用穩定性理論進一步討論了直立走滑型和傾斜逆沖型斷層地震的失穩機制，其中斷層軟化本構關系采用理論推導的負指數形式的軟化模式，它假設斷層局部微元強度遵循Weibull概率分布，采用統計力學方法推導得到宏觀本構關系.研究表明，當圍巖切線剛度小于斷層峰后軟化階段切線剛度絕對值的最大值時，整個系統的切線剛度為負，會發生地震失穩，且伴隨應力突跳和圍巖應變能釋放，其余情況僅僅是斷層無震滑動，不會發生地震.

以上理論模型均為單自由度問題，接近實際情況的地震模型都是多自由度的，需要采用數值方法進行模擬研究.王連捷等(2009，2011)將圍巖看作彈性體，斷層看成具有應變軟化特性的彈塑性體，斷層和圍巖組成統一的地質介質系統，采用有限差分FLAC-3D軟件分別模擬了汶川MS8.0地震以及玉樹MS7.1地震的發震過程，計算了應力降、能量釋放量、斷層錯動量、同震位移、震前位移、地震復發周期等重要參數，結果與野外調查資料具有較好的一致性.由于問題的復雜性，研究中所建立的數值模型均是高度簡化和概念性的.

2008年5月12日汶川MS8.0地震發生在青藏高原東緣的龍門山斷裂帶上，其發震動力學機制引起廣泛關注.龍門山斷裂帶在汶川地震前并沒有7.0級以上地震記錄(鄧起東等，2011; Ran et al., 2013).震前GPS觀測和地震地質研究結果都表明龍門山斷裂帶自晚更新世以來具有比較弱的活動性(Gan et al., 2007)，橫跨整個斷裂帶的滑動速率不超過2 mm·a-1，單條斷裂的滑動速率只有不到1 mm·a-1(張培震等，2008).汶川震后，科學家從構造變形和地震地質、數值模擬等不同角度對汶川地震的發震機制進行研究.張培震等(2009，2012)提出了汶川地震的多單元組合模式，認為地震的孕育和發生是川西高原變形單元、龍門山斷裂帶閉鎖單元和四川盆地支撐單元共同作用的結果.張竹琪等(2010)的研究結果顯示龍門山斷裂帶深部緩傾角、淺部高傾角逆沖斷層結構對汶川地震有促進作用.深部結構探測的結果表明汶川地震初始破裂段上地殼呈現高速異常，并且速度結構和物性結構差異是控制地震破裂分布的主要因素(鄧文澤等，2014).朱守彪和張培震(2009)、Zhu和Zhang(2010)、陳棋福等(2015)、馬林飛等(2018)采用數值模擬方法，通過建立橫跨龍門山斷裂帶的二維平面應變模型，研究了汶川地震的動力學機制，討論了摩擦系數和斷層傾角等對地震發生的影響.他們在研究中均采用速率-狀態相依的摩擦本構關系模擬斷層面上的摩擦機制.

因此，考慮到斷層的軟化特性對地震失穩的重要性，本文通過建立橫跨龍門山斷裂帶并包含川西高原和四川盆地在內的汶川二維平面地震力學模型，采用有限元方法數值計算得到了包含斷層和圍巖在內的整個巖石力學系統的平衡路徑曲線.從地震不穩定性角度詳細探討了汶川地震的孕育發生過程和失穩機制，分析了斷層傾角以及斷層材料參數對系統穩定性狀態和地震失穩過程的影響.

1 汶川地震不穩定性模型

1.1 模型的建立

汶川地震發生在龍門山斷裂帶中段，震中位于映秀鎮西南，同震破裂沿中央斷裂和前山斷裂以朝北東方向破裂為主(鄧起東等，2011).在地震震源破裂過程中，南段以逆沖為主，向北逐漸轉變為以走滑為主.如圖1所示，取垂直于龍門山斷裂帶并過汶川地震初始破裂點的剖面建立模型.由于剖面所在斷裂的同震破裂主要由垂直于斷層走向的水平擠壓導致，平面應變假設可近似成立(馬林飛等，2018).

圖1 汶川地震構造背景及模型剖面位置

過圖1中剖面建立垂直于龍門山斷裂帶，包含川西高原和四川盆地在內的汶川地震不穩定性二維平面應變地震力學模型，如圖2所示.模型取從地表到Moho面深度的縱剖面，四川盆地與川西高原以龍門山斷裂帶為界.龍門山斷層帶、川西高原和四川盆地共同組成統一的斷層圍巖巖石力學系統.川西高原地形高度取為4 km，四川盆地地形高度取500 m，龍門山斷裂帶高度取為3 km，向東逐漸降低，10 km后降到四川盆地高度.龍門山斷裂帶兩盤在地表位置的川西高原一側為點A，四川盆地一側為點A′.模型右側的四川盆地Moho面深度取為40 km，上地殼厚度取為22.5 km；左側的川西高原Moho面深度為70 km，上地殼和下地殼厚度分別為26 km和48 km.考慮到四川盆地所在的華南地塊相對川西高原運動速度很低，可視為固定不動，川西高原一側受構造加載的推擠作用.因此，在圖2右端面和四川盆地Moho面取法向約束，在川西高原左端面施加遠場位移a，Moho面取自由邊界.這樣的邊界條件模擬了下地殼物質向上運動的趨勢，有利于斷層解鎖.

圖2 汶川地震不穩定性地震力學模型示意圖

考慮到龍門山斷裂帶在空間幾何結構上具有高角度鏟形結構特性，淺地表傾角達60°～70°，甚至更陡(張培震等，2008)，因此在模型中將其設置成一條淺部高傾角、深部低傾角的鏟形斷層.為研究斷層傾角變化對汶川地震失穩的影響，設置了龍門山斷裂帶深淺部傾角不同的3組模型，幾何參數見表1所示.3組有限元模型如圖3所示.

表1 3組模型中龍門山斷裂帶幾何參數

圖3 龍門山斷裂帶深淺部不同傾角的3組有限元模型

1.2 本構關系和材料參數

巖石圈介質在漫長的地質演化歷史中具有流變性質，而在地震發生瞬時，則表現為彈性性質，本文采用Maxwell黏彈性本構模型模擬巖石圈介質在不同時間尺度的力學性質，本構關系為

(1)

其中，K(t)和G(t)分別為Maxwell模型的體積模量和剪切模量，η為黏滯系數，K和G分別為彈性體積模量和彈性介質模量，可以由線彈性介質的楊氏模量E和泊松比ν換算得到(王敏中等，2011).E和ν可根據研究區深部反演得到的三維P波、S波速度模型和密度模型計算得到(李平恩等，2019).地震波速結構(Li et al., 2013)和巖石圈三維流變結構反演結果(孫玉軍等，2013)顯示，相比青藏高原東部的川西高原地區，四川盆地的地殼介質較為堅硬，彈性模量更大，黏滯系數更高.綜合以上結果并參考其他學者的模型參數取值(朱守彪和張培震，2009；陳棋福等，2015；劉盼等，2017；馬林飛等，2018)，計算中模型材料參數見表2，地殼泊松比統一取為0.25.

表2 模型中川西高原和四川盆地分層材料參數

研究中考慮斷層的軟化特性，并且注意到遠場剛施加位移載荷a時，斷層并沒有錯動，只有當斷層受力達到足夠大數值時才會發生破裂，因此取龍門山斷裂帶具有初始破裂強度，只有當斷層面內剪應力τf達到一定數值時斷層才發生破裂進而開始錯動.斷層破裂面通常是壓剪性破裂，一般采用具有壓力相關性的Coulomb型破裂準則，它的剪切強度是

τf=c-μσn=c+μ|σn|,

(2)

式中σn是垂直斷層面的正壓力，按彈性力學約定，以拉應力為正，因此正壓力本身取負值.μ和c是斷層的材料參數，分別稱為斷層面的內摩擦系數和黏聚力，且有關系φ=arctanμ，φ是內摩擦角.如果斷層面內剪應力用τ表示，那么斷層的破裂準則為

f=τ-τf=τ-c+μσn=0,

(3)

它是一個雙參數(c和μ)的破裂準則.將斷層兩盤切向相對位移稱為斷層錯距u，為研究斷層的穩定性，至少需要假設一個強度參數(c或者μ)是斷層錯距u的減函數，可采用負指數形式的軟化模式(李平恩和殷有泉，2014)

(4)

或

(5)

式中c0和μ0分別是斷層的初始(當u=0時)黏聚力和內摩擦系數，m是強度曲線的形狀參數，有m≥1，u0是平均錯距的一種度量.由于u0的選取也能反映強度曲線的形狀，稱其為曲線的“胖度”(李平恩和殷有泉，2009).為研究的方便，同時取c和μ分別為如式(4)和式(5)所示的軟化模式，因此，描述斷層軟化特性的參數共有4個，分別是c0、μ0、u0和m.以大型通用有限元計算程序ADINA為計算平臺，采用二次開發技術開發了斷層軟化本構模塊.斷層采用間斷面單元進行模擬，除斷層外的其他介質采用Maxwell黏彈性本構模型描述.遠場加載速度取為2 mm·a-1，總加載量a等于加載速度乘以按年計算的加載時間，總加載時間取2000年.

2 地震不穩定性分析

在以上模型參數和加載條件下，給定不同的斷層軟化參數，可以計算得到遠場位移a與四川盆地一側右端面約束反力P的關系，它們在能量上是共軛的.對于包含斷層和圍巖在內的整個巖石力學系統，可以將a看作是輸入的控制變量，P是系統的響應.在以a為橫軸，P為縱軸的坐標平面內，每一個點(a,P)代表一個平衡狀態，整個a-P曲線代表所有的平衡狀態，稱為平衡路徑曲線，可以用來研究整個系統的穩定性(李平恩和殷有泉，2014).下面具體討論斷層傾角和斷層軟化參數對地震不穩定性的影響.

2.1 初始內摩擦系數的影響

為分析斷層初始內摩擦系數對地震失穩的影響，對于模型1，取c0=5000 Pa，u0=0.7 m，m=2時，分別計算了初始內摩擦角φ0=arctanμ0取不同值時的平衡路徑曲線，如圖4所示.為分析地震從孕育到發生過程中斷層兩盤的垂直位移變化特征，同時給出了龍門山斷裂帶兩盤地表處(川西高原一側為點A，四川盆地一側為點A′，如圖2所示)垂直位移隨時間的變化曲線，如圖5所示.由圖4可知，當初始內摩擦角為23°時，從平衡路徑曲線來看，隨遠場位移a的增加，約束反力P也緩慢變化，沒有發生突跳，表明斷層處于緩慢無震滑動狀態，沒有地震失穩發生.平衡路徑曲線有4個關鍵點，即C、E、F和G點，它們將平衡路徑曲線分成不同的段，各段對應的斷層運動狀態、應力和應變能變化特征不同.結合圖5可知，在點C之前，斷層上下盤垂向位移增速緩慢且增速相同，表明斷層處于閉鎖狀態，系統處于整體加載狀態，約束反力P隨遠場位移a不斷增大.隨遠場位移a繼續增加，在點C開始，斷層上盤垂直位移增速明顯加快，表明斷層克服閉鎖狀態開始啟動并緩慢無震滑動.在平衡路徑曲線應力峰值點E之前的CE段，約束反力P繼續增加，但增速明顯降低，表明在斷層緩慢無震滑動過程中儲存的應變能開始不斷釋放，使得應力增速降低，但釋放的應變能速度小于遠場位移加載增加的應變能速度，因此系統總應變能和應力是持續增加的.峰值應力點E也是平衡路徑曲線的轉向點，在轉向點E之后的EFG段，約束反力P開始降低，表明其間應變能釋放速度大于遠場位移加載增加應變能的速度，因此系統總應變能和應力持續減小.在E點之前，斷層上下盤垂直位移方向相同，均是垂直向上.過E點后，斷層下盤垂直位移開始變為垂直向下運動，斷層上下兩盤垂直運動方向相反.點F將EFG段進一步分為EF段和FG段兩部分.在EF段，應力和斷層兩盤垂直位移變化速率相對較緩，表明在該階段應變能釋放速度相對較慢.而在過點F之后的FG段，應力減小速度和斷層上盤垂直位移增加速度先顯著加快后逐漸變慢，表明應變能釋放速度先明顯加快后逐漸變慢，斷層無震滑動速度也先明顯加快后逐漸變慢.最后，在點G，約束反力P開始由減小狀態重新轉變為增加狀態，斷層上盤垂直位移增速降低，斷層下盤開始逐漸重新轉變為垂直向上運動.過點G后，系統應力和總應變能又開始增加，斷層上盤或下盤的垂直位移增速逐漸趨于一個穩定值，斷層進入勻速緩慢無震滑動狀態.

由圖4和圖5a可知，隨著初始內摩擦角從23°增大到30°，初始內摩擦系數增大，斷層面摩擦系數越大，斷層逐漸由緩慢無震滑動狀態轉變為閉鎖狀態.在此情況下，隨著遠場位移的不斷增加，斷層面內的剪應力越來越大，最終超過克服斷層閉鎖狀態所需的臨界剪應力，導致斷層突然錯動，伴隨應力突跳，系統失穩從而發生地震.當內摩擦角為30°時，從系統平衡路徑曲線來看，D是應力峰值點，即力的轉向點.在該點之前，系統的狀態是穩定的.D也是斷層啟動點，超過此點，斷層錯動加速，最終在A點失穩，同時地震發生.在A點失穩時，應力由A點突跳到B點，伴隨應力降、斷層錯距和應變能的急劇釋放.在失穩前的DA段，斷層錯動加速，應力開始降低，應變能開始釋放，是失穩的前兆.在失穩臨界點A之前，屬于地震孕育階段，其間圍巖應變能大量儲集，變形是漸變的.在臨界點A發生失穩時，系統的應力、斷層位錯和能量均發生突變.這反映了地震從孕育到發生是一個從漸變到突變的物理過程.地震失穩后過B點，應力和應變能又開始增加，斷層逐漸進入震后無震勻速緩慢滑動狀態.隨著初始內摩擦角從30°增加到35°，初始內摩擦系數進一步增大，地震失穩時對應的遠場位移更大，應力降和斷層錯距更大，體現出克服斷層閉鎖狀態需要的臨界剪應力更大，對應失穩時釋放的應變能更多，地震的震級更高.

圖4 初始內摩擦角取不同值時的平衡路徑曲線

圖5 初始內摩擦角取不同值時龍門山斷裂帶兩盤地表處垂直位移隨時間的變化

由圖5可知，對于初始內摩擦角為23°時的無震滑動狀態，上盤和下盤的垂直位移變化特征有所不同.具體是，在點C之前，斷層處于閉鎖狀態，上下盤垂直位移方向一致，均是垂直向上增加，且速率相等.在斷層克服閉鎖狀態開始無震滑動的CE段，上盤垂直位移增速相比閉鎖狀態時是增加的，而下盤雖然垂直位移增速相比閉鎖狀態時減小，但總的垂直位移仍然是增加的，表明在此階段，上下盤之間的黏結作用導致斷層兩盤垂直位移方向仍然保持一致，整體均處于垂直向上運動的狀態.過轉向點E之后，隨著應變能的加速釋放，上下盤之間的黏結作用降低，上盤和下盤的垂直運動方向開始反向，表現為上盤垂直向上運動，下盤垂直向下運動.過點G后，應力和應變能又開始增加，斷層上盤垂直位移增速逐漸降低并趨于穩定，斷層下盤由垂直向下運動逐漸轉變為垂直向上緩慢勻速運動，上下盤的垂直運動方向均為垂直向上，趨于一致.最后斷層保持緩慢勻速滑動.對于初始內摩擦角為30°時的地震失穩狀態，斷層兩盤在峰值應力點D之前一直處于閉鎖狀態，兩盤的垂直位移增速相同，方向垂直向上.過D點后，斷層克服閉鎖狀態開始啟動，并且加速錯動，斷層上下兩盤的垂直運動方向相反，上盤垂直向上、下盤垂直向下運動.最后，斷層兩盤垂直位移在失穩點A發生突跳，上盤川西高原一側垂直向上，下盤四川盆地一側垂直向下，體現出逆沖地震的特征，上盤垂直位移變化量明顯高于下盤.地震失穩后，斷層上盤很快恢復到垂直向上勻速運動狀態，斷層下盤也由垂直向下運動逐漸轉變到垂直向上勻速運動狀態，斷層進入震后緩慢勻速無震滑動狀態.隨著初始內摩擦角進一步增加，地震失穩突跳的時間延后，失穩時斷層錯距增大，失穩前的斷層加速錯動階段變短.由圖4和圖5進一步可知，初始內摩擦系數變化時，系統無論是斷層緩慢無震滑動狀態還是地震失穩狀態，系統重新進入應力和應變能增加狀態后，約束反力P隨遠場位移a的增速是一致的，并且斷層上盤或下盤的垂直位移增速也是一致的.

根據前述分析，汶川地震從孕育到發生的過程與圖4和圖5中對應的不穩定的地震失穩情況是相同的，考慮到其震級達到MS8.0，龍門山斷裂帶初始內摩擦系數可能較大.

2.2 斷裂帶傾角的影響

為分析龍門山斷裂帶深淺部傾角對地震失穩的影響，計算參數取c0=5000 Pa，u0=0.7 m，m=2，φ0=15°時，給出模型1至模型3的平衡路徑曲線，如圖6所示.由圖可知，當斷層深淺部總體平均傾角由50°增加到70°時，系統逐漸由穩定的無震滑動狀態轉變為地震失穩狀態，并且地震失穩時對應的遠場位移比無震滑動時斷層啟動點對應的遠場位移更大.表明在相同的孕震條件下，斷層傾角越大，越具備發生地震的條件，但需要更長的孕震時間；傾角越小，圍巖儲存的應變能通過斷層無震滑動的方式釋放掉，不容易發生地震.對于斷層無震滑動或地震失穩情況，系統重新進入應力和應變能增加狀態后，約束反力P隨遠場位移a的增速是不同的，與斷層傾角有關，總體而言，傾角越大，約束反力P隨遠場位移a的增速越大.

圖6 龍門山斷裂帶深淺部傾角取不同值時的平衡路徑曲線

當發生地震失穩時，為研究斷裂帶傾角的影響，取計算參數c0=10000 Pa，u0=0.7 m，m=2，φ0=45°時，給出模型1至模型3的平衡路徑曲線，如圖7所示.由圖可知，當斷層深淺部總體平均傾角由50°增加到70°時，失穩時對應的遠場位移增大，表明克服斷層閉鎖狀態需要的臨界剪應力更大，孕震時間更長；但突跳時的應力降先增加再減小.圖8給出了斷層兩盤地表處垂直位移隨時間的變化.由圖可知，斷層傾角越大，地震失穩發生的時間越晚.隨斷層傾角增加，對于斷層上盤，地震失穩時垂直位移變化量越大；對于斷層下盤，地震失穩時垂直位移變化量先減小再增大.斷層傾角越大，地震失穩時斷層錯距越大.地震失穩斷層進入緩慢勻速無震滑動狀態后，應力增速以及斷層上盤或下盤的垂直位移增速與斷層傾角有關，總體而言，斷層傾角越大，應力增速以及斷層兩盤垂直位移增速也稍大.

圖7 當地震失穩時龍門山斷裂帶深淺部傾角取不同值時的平衡路徑曲線

圖8 當地震失穩時龍門山斷裂帶深淺部傾角取不同值時斷層兩盤地表處垂直位移隨時間的變化

2.3 初始黏聚力的影響

為分析斷層初始黏聚力對地震失穩的影響，對于模型1，取計算參數為u0=0.7 m，m=2，φ0=25°時，分別計算了初始黏聚力c0取不同值時的平衡路徑曲線，如圖9所示.由圖可知，當初始黏聚力增加時，斷層逐漸從穩定的緩慢無震滑動狀態轉變到失穩狀態.并且初始黏聚力越大，地震失穩突跳時對應的遠場位移越大，應力降也越大.這表明斷層初始黏聚力越大，斷層更容易處于閉鎖狀態，閉鎖時間越長，克服閉鎖狀態需要的剪應力也越大.在斷層閉鎖期間，隨遠場位移的增加，斷層附近圍巖持續變形，將有利于積累更多的應變能.當斷層面內剪應力超過克服閉鎖狀態需要的臨界剪應力時，圍巖儲存的大量應變能容易在極短的時間內得到釋放，伴隨應力降和斷層突然錯動，導致地震發生.初始黏聚力越小，斷層克服閉鎖狀態需要的剪應力越小，斷層容易滑動，這導致斷層附近圍巖不容易發生持續變形，進而積累更多的應變能.圍巖存儲的應變能更容易在斷層緩慢無震滑動中釋放掉，不會發生地震.對于汶川MS8.0特大地震而言，斷層的初始黏聚力可能非常大，這樣導致地震的孕育時間非常長，能積累極大的應變能.

圖9 初始黏聚力取不同值時的平衡路徑曲線

2.4 強度曲線胖度參數的影響

為分析強度曲線胖度參數對地震失穩的影響，對于模型1，取計算參數為c0=5000 Pa，m=2，φ0=25°時，分別計算了胖度參數u0取不同值時的平衡路徑曲線，如圖10所示.圖11給出胖度參數取不同值時的斷層兩盤地表處垂直位移隨時間的變化.由圖10可知，胖度參數增加時，斷層從不穩定的失穩狀態逐漸轉變到穩定的無震滑動狀態，并且斷層無震滑動的啟動點對應的遠場位移增大.表明胖度參數增大會導致斷層趨于無震滑動.由圖10和圖11可知，盡管胖度參數越小，越容易發生地震失穩，但斷層閉鎖狀態的時間也越短，失穩點越提前，地震發生時的應力降和斷層錯距越小，釋放的應變能也越小，相應的孕震時間也越短.考慮到汶川地震的震級達到MS8.0，不同研究結果均顯示其復發周期長達千年量級(張培震等，2008；朱守彪和張培震，2009)，因此龍門山斷裂帶的胖度參數可能較大.

圖10 強度曲線胖度參數取不同值時的平衡路徑曲線

圖11 強度曲線胖度參數取不同值時龍門山斷裂帶兩盤地表處垂直位移隨時間的變化

2.5 強度曲線形狀參數的影響

為分析強度曲線形狀參數對地震失穩的影響，對于模型1，取計算參數為c0=5000 Pa，u0=0.7 m，φ0=25°時，分別計算了形狀參數m取不同值時的平衡路徑曲線，如圖12所示.由圖可知，隨著形狀參數從1增加到7，斷層從穩定的緩慢無震滑動狀態逐漸轉變到不穩定的失穩狀態.當發生地震失穩時，隨形狀參數增大，失穩點對應的時間和地震發生時的應力降也有所增加，但增加量逐漸減小.形狀參數m越大，材料脆性越大，當m趨于無窮大時，對應理想脆性材料(李平恩和殷有泉，2009).因此，當m增大到一定值時，繼續增加m值，地震失穩時對應的時間和應力降不再發生顯著變化，如圖12所示.

圖12 強度曲線形狀參數取不同值時的平衡路徑曲線

3 討論

考慮斷層的軟化特性時，斷層傾角以及斷層材料參數對地震穩定性均有影響.當這些參數變化時，整個巖石力學系統會表現為穩定和不穩定2種不同狀態.

對于穩定狀態情況，如圖4至圖5，斷層初始處于閉鎖狀態，隨著遠場位移的增大，斷層會克服閉鎖狀態開始緩慢無震滑動，其間圍巖中儲存的應變能會不斷釋放，但釋放的應變能速度小于遠場位移加載增加的應變能速度，因此雖然應力增速降低，但應力是增加的，系統總應變能也是增加的.當達到應力峰值點之后，斷層無震滑動速度逐漸加快，應變能釋放速度也逐漸增大，且應變能釋放速度超過遠場加載增加的應變能速度，系統總應變能不斷減少，應力也不斷降低，表明存儲的應變能通過斷層無震滑動的方式快速釋放，這個能量大部分轉化為斷層摩擦滑動產生的熱.隨著應變能的不斷釋放，應力減小速度和斷層滑動速度不斷降低.最后，在遠場位移作用下，應力和總應變能又開始勻速增加，斷層進入緩慢勻速無震滑動狀態.

對于不穩定狀態情況，如圖4至圖5，初始斷層一直處于閉鎖狀態，當達到應力峰值點之后，斷層開始克服閉鎖狀態啟動，并且錯動加速，應力開始降低，應變能開始釋放.應力峰值點是應力增加和減小的分界點，過該點后，隨著斷層錯動加速，最終在失穩點發生地震，儲存在圍巖中的應變能在極短的時間內急劇釋放，同時發生應力突跳，伴隨應力降和斷層突然錯動產生斷層錯距.在失穩點之前屬于地震孕育階段，應力、變形和應變能的變化是漸變的.在失穩點發生地震失穩時，應力、變形和應變能的變化是突變的，是地震的發生階段.系統穩定性分析充分反映了地震從孕育到發生是一種從漸變到突變的物理過程.地震發生后，應力和總應變能又開始增加，斷層進入震后緩慢無震勻速滑動狀態.研究得到一個地震周期內整個巖石力學系統的應變能從積累到突然釋放的全過程，與馬林飛等(2018)采用速率-狀態相依斷層摩擦本構的模擬結果吻合.

對于地震失穩情況的平衡路徑曲線，在應力峰值點和失穩點之間是斷層錯動加速階段，是失穩的前兆，其間應力開始連續降低，應變能開始連續釋放.斷層傾角和斷層材料參數取不同值時的計算結果顯示，在大部分情況下，該階段存在的時間是短暫的，在有些情況下應力峰值點和失穩點非常接近，這表明斷層剛克服閉鎖狀態進入加速錯動階段后就很快發生地震失穩，導致失穩前兆的存在時間異常短暫.

當斷層材料參數，即初始內摩擦系數、初始黏聚力、強度曲線胖度和形狀參數變化時，可能影響整個巖石力學系統的穩定性狀態，如圖4至圖5、圖9至圖12所示，但無論是穩定的斷層無震滑動狀態還是不穩定的地震失穩狀態，當應力和總應變能重新進入增加狀態后，約束反力的增速是相同的，斷層上盤或下盤垂直位移的增速也是相同的.而當斷層傾角不同時，如圖6至圖8所示，當應力和總應變能重新進入增加狀態后，約束反力的增速以及斷層兩盤垂直位移的增速是不同的.這表明應力和應變能的增速由遠場加載速度、巖石力學系統的結構和圍巖材料屬性決定，與斷層軟化特性參數無關.研究結果顯示，斷層淺部為高傾角70°時，深部為緩傾角35°時，巖石力學系統也會發生失穩，這從數值模擬角度支持了宮猛等(2020)得到的汶川地震破裂斷層面傾角為35°的結果.

文中沒有討論重力的影響，這是因為本文重點關注的是斷層軟化特性對地震失穩的影響.定性來講，考慮重力時會引起斷層面上的正應力σn和剪應力τ發生變化，進而通過式(3)影響斷層的破裂狀態.但總體而言，重力并沒有直接改變斷層的軟化本構特性，因為斷層黏聚力式(4)和內摩擦系數式(5)并沒有改變，因此重力對地震穩定性狀態的影響是通過改變斷層面上的應力間接實現的，并非導致地震失穩的本質原因，具體影響將在今后的工作中進一步討論.

川西高原和四川盆地的材料參數的不確定性對地震穩定性狀態的影響與重力的情況類似，都是通過改變斷層面上的應力間接影響地震失穩，并不會改變整個巖石力學系統穩定和不穩定的2種狀態，對于不穩定狀態，僅僅會影響地震發生失穩的時間和地震失穩時的斷層錯距等.姚琪等(2012)分析了斷層兩盤巖性差異對汶川地震的影響.實際計算中，我們也對多組川西高原和四川盆地的材料參數模型進行了分析，進一步驗證了以上結論.

汶川地震從孕育到發生的過程和失穩機制與圖4和圖5中的地震失穩情況是相同的，僅僅是平衡路徑曲線的具體形態可能存在差異.基于前面的分析，總體而言，龍門山斷裂帶的高傾角鏟形結構有助于孕育高達MS8.0的逆沖型地震.同時，龍門山斷裂帶的初始內摩擦系數μ0、初始黏聚力c0和強度曲線形狀參數m可能都較大.盡管強度曲線胖度參數u0越大時更容易導致斷層無震滑動，考慮到u0增加會導致斷層閉鎖時間增加，有利于積累更多的應變能，因此推測龍門山斷裂帶的u0值可能也較大.龍門山斷裂帶4個軟化材料參數的綜合效果是有利于整個巖石力學系統增加斷層閉鎖時間(地震孕育時間)，積累更多的應變能，以致于最終發生汶川MS8.0地震.

構成斷層的介質均屬于巖石類介質.大量三軸巖石力學實驗結果顯示，巖石在加載受力破壞過程中，當巖石的應力達到峰值屈服強度之后，隨著變形繼續增加，巖石強度會發生劣化或降低，稱為巖石的軟化特性.考慮到巖石介質的軟化特性，馬瑾和郭彥雙(2014)開展了平直斷層失穩的實驗室模擬研究.實驗加載時在X方向壓力保持不變，Y方向按位移控制方式加載，且位移速率保持不變，最后得到差應力-時間的全過程曲線，完整體現了地震從孕育到突然失穩的全過程.本文平衡路徑曲線的橫坐標是遠場位移加載總量，由于是勻速加載，因此橫坐標與時間線性相關.縱坐標是四川盆地一側右端面約束反力P，由于模型的上表面和底面均自由，因此約束反力P可視為垂向和水平向的差應力.因此把平衡路徑曲線轉換為垂向和水平向的差應力-時間曲線時，形態保持不變.本文數值模擬給出的地震失穩狀態下的平衡路徑曲線從形態上與馬瑾和郭彥雙(2014)的差應力-時間曲線完全一致，數值模擬結果與物理實驗結果是吻合的.

最近，有學者采用基于區域多個震源機制解反演定量應力張量的方法對汶川地震前孕震區的應力狀態進行了研究，推斷孕震區在2007年中期達到應力峰值，并且獲得了顯著證據證明在峰值前應力隨時間平穩變化并以積累為主，而在峰值后則轉變為明顯的應力釋放階段(王凱英等，2018).這表明汶川地震前，區域應力存在由積累為主轉為釋放為主的可觀測階段.該研究結果與地震失穩情況下在應力峰值點前后的應力變化特征是吻合的.本文為汶川地震前的應力變化研究結果從機理上給出了一種解釋.

4 結論

本文考慮斷層的軟化特性，通過建立垂直于龍門山斷裂帶的汶川地震平面地震力學模型，從系統穩定性角度研究了汶川地震的孕育發生過程和失穩機制，探討了斷層傾角和斷層材料參數變化對系統穩定狀態和地震失穩過程的影響.研究表明，當斷層傾角或斷層材料參數變化時，系統可能表現為穩定的斷層緩慢無震滑動或不穩定的地震失穩2種狀態.初始內摩擦系數、斷層傾角、初始黏聚力和強度曲線形狀參數增大、或者強度曲線胖度參數減小均會導致整個巖石力學系統更趨向于不穩定的地震失穩狀態.當發生地震失穩時，初始內摩擦系數越大、初始黏聚力越大，失穩時的應力降越大，相應地震震級更高.在相同的孕震條件下，斷層傾角增加有助于地震的發生；傾角越小，儲存的應變能可能通過斷層無震滑動的方式釋放掉，不容易發生地震.地震失穩前在應力峰值點和失穩點之間是斷層加速錯動的階段，在該階段應力開始減小，應變能開始釋放，是失穩的前兆.無論是穩定的斷層無震滑動狀態還是不穩定的地震失穩狀態，當應力和總應變能重新進入增加狀態后，應力和斷層兩盤位移變化速率由遠場加載速度、巖石力學系統的結構和圍巖材料屬性決定，與斷層的軟化特性參數無關.

致謝感謝三位匿名審稿專家為提高本文質量給出的有益意見和建議.