彎折斷層系統上超剪切破裂傳播的產生條件

鄭玲瓏, 錢峰, 張海明

北京大學地球與空間科學學院地球物理系, 北京 100871

0 引言

地震發生過程中斷層的破裂速度是控制地震波形態的重要因素，對于近場的強地面運動有顯著的影響.隨著觀測記錄的增多，在越來越多的地震中發現了超剪切傳播的證據.已被發現出現超剪切破裂的地震包括1979年Imperial Valley地震(Archuleta, 1984)，1999年Duzce地震(Bouchon et al., 2001)，1999年Kocaeli地震(Bouchon et al., 2001)，2001年昆侖山口西地震(Bouchon and Vallée, 2003)，2002年Denali地震(Dunham and Archuleta, 2004)，2010年玉樹地震(Wang and Mori, 2012)，以及2013年Craig地震(Yue et al., 2013)等.在超剪切破裂過程中，斷層不同時刻輻射的剪切波互相干涉、疊加，形成馬赫波，使得地震波振幅更大，輻射范圍更廣，往往造成更嚴重的地震災害(Bernard and Baumont, 2005; Dunham and Archuleta, 2005; Dunham and Bhat, 2008).多數破壞性的地震往往都發生在幾何形態復雜的斷層系統上，由于斷層空間延展尺度較大，斷層的幾何形狀將會對破裂過程產生重要的影響(馬瑾等，1996).研究復雜斷層系統上的超剪切破裂轉化的條件，對于震源過程的科學認識和防震減災工作都具有重要的意義.

B-A機制是針對二維模型提出的，認為無量綱比值S臨界值為1.77，只有滿足一定的背景應力和摩擦條件時，才能夠產生超剪切破裂(Andrews, 1985).針對三維均勻介質模型，Madariaga和Olsen(2000)提出了另一種無量綱化參數κ作為確定破裂速度的指標，定義為破裂應力釋放能與斷裂能之比，包含了更多破裂信息，適用于高應力場斷層破裂情況；而S則需要小于更小的臨界值1.19，才能觸發超剪切破裂(Dunham, 2007).考慮包含地表的均勻半空間介質，S值大于1.19也產生超剪切破裂(Bizzarri, 2010; Bizzarri and Das, 2012; Zhang and Chen, 2006).Olsen等(1997)、Kaneko和Lapusta(2010)認為自由界面SV-P衍射波的相變轉換所產生的附加應力場是超剪切破裂形成的主要原因，而走滑斷層自由表面所誘發的超剪切破裂也可能會返回至次瑞利速度(Hu et al., 2019).即使不考慮自由界面的效應，三維全空間斷層介質的非均勻性使得在不滿足B-A機制所預測的背景應力水平下仍然有可能出現超剪切破裂(Passelègue et al., 2013)，但也可能使得當背景預應力超過所預測的臨界水平時仍然以亞瑞利波速穩定傳播，無法出現超剪切破裂(Fukuyama and Madariaga, 2000).

影響超剪切破裂出現的復雜因素還包括，非成核區介質摩擦性質和背景預應力分布的不均勻性，以及斷層幾何形狀的復雜性等.一般以包含障礙體(barrier)和凹凸體(asperity)的斷層面的粗糙度和準靜態裂縫表示非成核區介質摩擦性質和背景預應力的非均勻性.Fukuyama和Olesn(2002)研究了不同凹凸體寬度和不同背景預應力對超剪切破裂的影響，發現在一定條件下，含有凹凸體的非均勻介質相對于均勻介質，更容易激發超剪切破裂；對于寬度較小的凹凸體需要更高的背景預應力來觸發超剪切破裂，當凹凸體寬度大于某一臨界值后，所需要的背景預應力不再降低，保持穩定后仍然可以產生超剪切破裂.當非成核區存在準靜態二次擴張裂縫時，根據B-A機制，二次裂縫將會受到主破裂前S波和P波應力場的影響開始破裂，當主破裂前鋒與二次裂縫前端相遇時速度會突然超過剪切波速，轉換成超剪切破裂(Liu and Lapusta, 2008).針對高剪切強度的障礙體，Dunham等(2003)提出在障礙體與成核區之間具有一定距離時，破裂前鋒繞過障礙體，其產生的凹形聚焦效應能夠使破裂傳播速度迅速增加到超剪切狀態.之后的研究發現障礙體誘導形成超剪切破裂需要滿足一定的條件(Duan, 2012; Yang et al., 2012, 2013; Weng et al., 2015)，障礙體埋深越淺、尺寸越大和強度越高, 越容易產生超剪切破裂，當障礙體位于自由界面，與成核區在自由界面上的距離大于一定長度時，能夠擴大超剪切破裂的傳播距離(Xu et al., 2016).當斷層面本身大面積凹凸不平，粗糙度越大越容易產生局部超剪切破裂，但持續的超剪切往往出現在粗糙度相對中等的位置，過大的粗糙度會損耗能量，使超剪切破裂無法持續(Bruhat et al., 2016).

斷層的幾何形狀是影響超剪切破裂產生的另一個重要因素.目前涉及到斷層幾何形狀的研究主要集中于幾種端元的復雜斷層模型：階躍斷層、彎折斷層和分叉斷層.關于階躍斷層上超剪切破裂的影響因素研究較為完善，當滿足一定范圍內的斷層間隔和足夠長度的平行重疊距離后，階躍位置的S值會發生變化降低到臨界值之下，積蓄能量在次級斷層上產生超剪切破裂(Kase and Kuge, 1998; Ryan and Oglesby, 2014; Hu et al., 2016; 袁杰和朱守彪，2017).因為彎折和分叉斷層系統的復雜性，在數值模擬中相對較難實現，故而對其產生超剪切破裂的研究更加困難.Aochi等(2000a,b)基于邊界積分方程法(BIEM)對簡單彎折和分叉斷層進行了數值模擬，發現當光滑彎曲斷層面的最終角度為π/18時，破裂能夠以超剪切速度傳播，初始剪切應力和摩擦參數之間的關系對于簡單分叉斷層系統動態破裂傳播的斷層選擇是最重要的，法向應力的變化通常較小.不同的彎折角度對破裂的傳播和停止也有著重要影響(張麗芬等，2016).Kame等(2003)通過邊界積分方法研究了二維分叉斷層破裂速度影響條件，發現只有當分叉斷層的S值從3突降到臨界值以下的0.8時，才能在斷層上產生超剪切破裂.因此對于彎折和分叉斷層來說，B-A機制所提出的影響因素仍占據主導地位，針對不同形狀的彎曲斷層可能有著不同的臨界狀態.對于分叉斷層系統，在一定條件下分叉間的相互作用也有可能使其中一支分叉斷層產生超剪切破裂(Bhat et al., 2004).

目前有關復雜斷層系統上的超剪切破裂的研究較少.本文將基于三維空間中非結構化的BIEM，針對不同彎折角度的彎折斷層進行大量高效的數值模擬，在參數空間中通過破裂前鋒傳播圖和破裂速度相圖研究彎折斷層的破裂速度變化，分析彎折斷層的破裂行為和超剪切破裂的產生條件.

1 自發破裂模擬方法及破裂相圖算法

1.1 邊界積分方程方法(BIEM)

本研究采用的自發破裂數值模擬方法為BIEM，它建立了斷層應力與滑動速率的關系(Das and Aki, 1977; Aochi et al., 2000a,b)：

(1)

(2)

對于斷層自發破裂過程而言，必須結合控制破裂行為的摩擦準則.常用的摩擦準則有滑動弱化摩擦準則(Ida, 1972)和速率-狀態摩擦準則(Dieterich, 1979).對于同震破裂過程的模擬，通常采用滑動弱化摩擦準則：

(3)

其中，τ是當前時刻單元的剪切應力值，τu和τf分別是破裂強度和殘余應力，Dc為臨界滑動弱化位移，D為當前時刻的單元累計滑動量.當斷層單元上的剪切應力累積到了破裂強度時，單元開始破裂，其剪切應力隨著滑動量線性減小，直至達到穩定滑動后保持不變的摩擦力.將BIEM和滑動弱化準則聯立，給定初始應力后就可計算得到破裂過程.

1.2 破裂速度自動提取

完整破裂相圖的參數空間包含幾千個破裂結果，需要逐個判斷每個算例的破裂速度.為了提高效率，本文采用了適用于全空間均勻走滑斷層模型的破裂速度自動提取方法，能夠在破裂傳播計算過程中得到算例的破裂速度，步驟如下：

圖1 破裂速度拾取示意圖 (a) 平面斷層模型，黃色高亮單元為破裂前鋒中心點傳播路徑的網格單元, 綠色單元為成核區單元. (b)平面斷層滑動速率時空分布，圖中VP和VS兩條白線分別代表P波和S波的傳播速度，紅圈為破裂前鋒隨時間變化的x方向位置，綠 線為紅圈的擬合線條，其斜率近似反映破裂的傳播速度.Fig.1 Illustration of the rupture speed picking (a) The planar fault model. The highlighted yellow elements are the grids on the propagation path of the midpoint of the rupture front, the green element is a grid inside the nucleation zone. (b) The spatiotemporal evolution of the slip velocity of the planar fault model. The white lines denote the reference P and S wave velocity respectively. The red circles are the locations of rupture fronts along the fault varying with time, whereas the slope of the fitted connecting green line approximates the rupture propagation speed.

(1)圖1a中顯示為綠色單元的成核區(初始破裂區域)近似為圓心坐標為(22.5,7.5)，半徑為3 km的圓，根據成核區的位置，標記破裂前鋒的中心點在斷層上的傳播路徑中黃線所示，對路徑上單元按照x方向坐標從小到大排列；

(2)在算例計算過程的每個時間步循環內，對步驟1中已排序的傳播路徑單元的滑動速率值進行自動篩選，找到成核區沿傳播方向上滑動速率降為0的最小x坐標單元，即破裂前鋒對應單元，并記錄時間及該單元x方向坐標，用圖1b中的紅圈標記；

(3)每隔一定時間間隔對步驟2中記錄單元的坐標和時間進行x-t線性擬合，擬合結果可示意為圖1b中紅圈擬合的綠線，計算擬合線的斜率得到破裂傳播速度，當破裂的時間循環計算結束時就可自動得到一系列實時變化的破裂速度.

1.3 平面斷層的相圖對比

為了檢驗破裂速度自動拾取方案的正確性，我們使用與Xu等(2015)近似相同的參數(表1)和相同的無量綱化公式(4)計算得到平面斷層破裂相圖：

(4)

圖2顯示了本文與前人的平面斷層破裂相圖，對比發現平面斷層破裂相圖整體趨勢相同，驗證了破裂過程自動拾取的正確性.

圖2 平面斷層破裂相圖 左(本文)，右(Xu et al., 2015).0～0.2區間內過小，導致計算結果不可靠，故整體在破裂相圖中以灰色表示， Supershear區域表示算例為超剪切破裂， SubRayleigh區域表示算例破裂速度為亞瑞利波速，Self-arresting區域表示算例破裂未持續傳播自發停止.Fig.2 Phase diagram of rupture on a planar fault (Left) Result in this study, and (right) that in Xu et al. (2015). Since the small value in the range 0～0.2 would produce unreliable results, this part of parameter space is shaded by gray. Region in which super-shear ruptures occur is marked in “Supershear”. Rupture speed below Rayleigh wave speed, “SubRayleigh”, and rupture ceases spontaneously, “Self-arresting”.

表1 平面斷層模型參數Table 1 Parameters of the planar fault models

2 彎折斷層破裂模擬結果與討論

利用前面所描述的方法，本文分別對彎折角度為15°、25° 和40° 的彎折斷層模型進行破裂相圖計算.彎折斷層模型如圖3所示，斷層模型從x方向45 km處開始彎折，整個斷層模型由平面部分15 km×45 km的矩形平面和彎折平面部分為15 km×15 km的正方形平面以彎折角α拼接而成，成核區位于矩形平面的中央.

圖3 彎折角α=15°的彎折斷層模型 黃色三角形單元為圖6所取彎折面上靠近彎折交界線單元.Fig.3 The fault model with a bending angle α 15° The highlighted yellow triangle element is the grid on the bending plane chosen in Fig.6 that is close to the bending boundary.

(5)

其中：τa=σ33(μs-μd)，σ33為彎折斷層平面部分初始正應力.

表2 彎折斷層模型參數Table 2 Parameters of the fault models with a bending angle

觀察三個破裂相圖的自發停止破裂、亞瑞利破裂和超剪切破裂之間的兩條分界線I和II可以看出，自發停止破裂區域范圍在三個不同角度的彎折斷層相圖中完全一致，不受彎折角度改變影響；而三個相圖中分界線II的不同，則顯示當破裂能夠穩定傳播時，彎折面的角度對破裂速度有明顯影響.隨著彎折角度增加，超剪切區域范圍減小，亞瑞利破裂范圍增大，15° 與25° 彎折斷層破裂相圖的差異較小，彎折角度為40° 時亞瑞利波破裂范圍明顯增加，且超剪切破裂邊界線變化趨勢與平面破裂相圖較為相似.

圖4 彎折斷層破裂相圖對比圖 (a)、(b)、(c)分別為彎折角α=15°，25° 和40° 的結果.綠色虛線I為自發停止破裂和亞瑞利破裂之間的分界線，藍色虛線II為亞瑞利破裂和超剪切破裂之間的分界線，Sub-Rayleigh 1代表算例破裂速度為亞瑞利波速且破裂未傳播至彎折面，Sub-Rayleigh 2代表破裂速度 接近但不超過剪切波速的算例，圖中A1—I6所圈參數的算例依次對應圖6相同編號的子圖.Fig.4 Comparison of phase diagrams of rupture on faults with different bending angles (a), (b) and (c) are results with the bending angle of 15°, 25° and 40°, respectively. The dashed green line I separates the Self-arresting and SubRayleigh. The dashed blue line II separates the SubRayleigh and Supershear. Examples marked in “SubRayleigh 1” represent ruptures whose speed are below the Rayleigh wave speed. Examples marked in “SubRayleigh 2” represent ruptures whose speed are approaching, but not beyond the shear wave velocity. The examples circled and annotated with A1—I6 correspond to each subplot in order in Fig.6.

圖5a、b、c分別為平面初始正應力(σ33,σ11)設置為(-10, -26) MPa，(-10, -15) MPa, (-10, -10) MPa的不同彎折面應力莫爾圓示意圖，每張子圖展示了給定同一平面初始應力時，不同彎折角的斷層面初始應力和S.可以看出初始狀態下對于15°、25°和40°彎折斷層面，α越小的斷層面越靠近破裂線τp，參數S越小，意味著破裂越容易傳播，也就越容易產生超剪切破裂.圖5a為圖4相圖初始應力設置，與圖4的彎折角越小超剪切破裂區域越大規律一致.三幅子圖中圖5a即本文所采用應力設置，初始狀態下各個彎折面對應S值最小，作為對照組的圖5b、c的應力設置使得平面和彎折斷層在莫爾圓上位置整體左移，逐漸偏離破裂線.

圖5 斷層面不同初始正應力莫爾圓示意圖 (a)、(b)、(c)分別為平面初始正應力(σ33,σ11)設置為(-10, -26) MPa，(-10, -15) MPa, (-10, -10) MPa的不同彎折面應力莫爾圓示意圖，三個彎折模型除彎折角α外的其他參數相同，所取為0.85，S=(τp-τo)/(τo-τr)，τp為破裂強度，τo為初始剪應力，τr為剩余剪應力，σn為正應力.圖中彎折角α=0°的灰色正方形所標記的橫縱坐標展示平面斷層的初始正應力和剪應力， 其他形狀標記對應了 不同彎折角的彎折面初始正應力和剪應力.Fig.5 Diagram of Mohr′s circle under different initial normal stresses on the fault model The planar initial normal stress (σ33,σ11) for the Mohr′s circle diagram (a), (b) and (c) are set to be (-10, -26) MPa, (-10, -15) MPa and (-10, -10) MPa, respectively. All parameters but the bending angle α are the same for three models. is set to be 0.85 and the ratio S=(τp-τo)/(τo-τr), where τp is the rupture strength, τo is the initial shear stress, τr is the residual shear stress, and σn is the normal stress. The case α=0° is denoted by the gray square, whose horizontal and vertical coordinate represent the initial normal stress and shear stress on the planar part, respectively. Other symbols correspond to the initial normal stress and shear stress on the bending part， with different bending angles.

圖6 15°、25°和40° 彎折斷層分別取0.35, 0.70和0.93時的滑動速率時空分布組圖 每個子圖代表一種參數組合下結果，子圖左下角所標序號對應著圖4相應相圖中所標的序號，紅圈代表破裂前鋒，為了方便判斷破裂 前鋒的速度，標出了P波和S波速度.破裂算例對應參數以標注在白色括號內，白色虛線為彎折交界線x=45 km的標記.Fig.6 The selected spatiotemporal evolution of slip velocity in 15°-, 25°- and 40°- bending fault when is 0.35, 0.70 and 0.93, respectively Each subplot represents the result under specific combination of parameters which is marked in white brackets as The sequential number at the lower left corner of each subplot is the same as that in Fig.4. Red circles denote the rupture front. P and S wave velocities are illustrated for reference. The white dashed line marks the location of bending boundary (x=45 km).

圖7 彎折斷層相同初始正應力的滑動速率時空分布對比圖 (a)、(b)、(c)分別為平面初始正應力(σ33,σ11)設置為(-10, -15) MPa的彎折角α=15°，25° 和40° 的彎折斷層滑動速率時空分布圖，子圖所取為(0.70, 0.85)，圖中VP和VS兩條白線分別代表P波和S波的傳播速度，紅圈示意為破裂前鋒隨時間的變化 趨勢，綠線為紅圈的擬合線條，其斜率近似反映破裂的傳播速度.Fig.7 Comparison of the spatiotemporal evolution of slip velocity under the same initial normal stresses on the bending fault model The planar initial normal stress (σ33,σ11) for the spatiotemporal evolution of slip velocity on the model with a bending angle of 15°, 25° and 40° are set to be (-10,-15) MPa, respectively. The parameters for the subfigures are (0.70,0.85), respectively. The white lines denote the reference P and S wave velocity, respectively. The red circles are the locations of rupture fronts along the fault varying with time, whereas the slope of the fitted connecting green line approximates the rupture propagation speed.

圖8 彎折面上交界線附近一個單元的不同彎折角度斷層模型應力時間演化對比圖 (a)、(b)、(c)分別為平面初始正應力(σ33,σ11)設置為(-10,-26) MPa，(-10,-15) MPa, (-10,-10) MPa的不同彎折角度斷層模型應力時間演化對比圖.每個子圖中三個彎折模型算例除彎折角度外的其他參數相同，為(0.70, 0.85)，Te為 剪切應力，Tu為破裂強度，N為正應力.Fig.8 Comparison of the temporal variation of stress at one element grid close to the bending line on the fault model with different bending angles (a), (b) and (c) correspond to the comparison of the temporal variation of stress on the fault model with different bending angles when the planar initial normal stress (σ33,σ11) are set to be (-10,-26) MPa, (-10,-15) MPa and (-10,-10) MPa, respectively. All parameters but bending angles are the same for three subfigures. is (0.70,0.85), where Te,Tu and N are the shear stress, the fracture strength and the normal stress, respectively.

3 結論

對不同彎折角度的相圖和選定參數組合下的滑動速率時空圖分析，可以發現對于本文所設應力條件下彎折面處于壓縮區(α>0)的彎折斷層模型，有以下結論：

(1)彎折角度對于無法持續傳播的自發停止破裂沒有影響，相圖中自發停止破裂區域所占參數空間完全相同；

(3)體現在相圖上，15° 和25° 的超剪切區域明顯比40° 的大，即破裂能夠越過彎折交界線在彎折面上繼續傳播時，更容易形成超剪切破裂； 15° 的超剪切區域又略大于25°，即彎折角越小，越容易形成超剪切破裂；

(4)對初始剪應力較小的情況，超剪切破裂更容易在彎折面上形成，而當初始剪應力較大時，超剪切破裂在平面上就可形成；

根據滑動速率時空分布圖，可以直觀地分析破裂傳播速度和滑動速率的變化，當破裂傳播至彎折面上時，需要克服彎折交界線上的正應力，從而破裂前鋒在交界線處出現錯位延遲效應，同時在彎折面上的滑動速率相對于平面斷層有了明顯的增大.彎折角度越大，需要克服的時間越長，延遲效應越明顯.

本文考慮的全空間彎折斷層模型，包括平面部分和相對彎折部分，破裂成核區設在平面斷層，彎折斷層部分位于平面斷層的壓縮區，這只是比平面斷層稍復雜一點的斷層系統，本文進行了該類斷層系統的破裂行為研究的初步嘗試.對于其他類型的復雜斷層，比如分叉斷層，甚至更為復雜的斷層系統，其破裂行為的控制因素將更為復雜，這將在進一步的研究中予以考慮.