馬聰慧, 錢峰, 張海明

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

0 引言

2013年4月20日四川省雅安市蘆山縣發(fā)生了MS7.0地震，這是繼汶川地震以來(lái)，在龍門山斷裂帶上的又一次強(qiáng)震，引起了嚴(yán)重的山體崩塌、滑坡、泥石流、砂土液化等次生災(zāi)害，造成了重大的人員傷亡和財(cái)產(chǎn)損失(Zhang et al., 2013).震后野外應(yīng)急科考和淺層人工地震勘探調(diào)查發(fā)現(xiàn)，蘆山地震極震區(qū)和附近斷裂帶沿線均沒(méi)有出現(xiàn)構(gòu)造成因的地震地表破裂帶，推斷此次地震發(fā)震斷層屬于隱伏斷層(李傳友等, 2013; 徐錫偉等, 2013; 雷生學(xué)等, 2014).不同研究者采用近震、遠(yuǎn)震數(shù)據(jù)反演蘆山地震震源機(jī)制解，結(jié)果表明，主震發(fā)震斷層走向約210°，滑動(dòng)角近90°，可視為純逆沖型地震；由于使用的地震數(shù)據(jù)或?yàn)V波頻段不同，斷層傾角和發(fā)震深度結(jié)果有較大差異，整體而言，斷層傾角約為42°，主震深度約15 km(劉杰等, 2013; 呂堅(jiān)等, 2013；謝祖軍等, 2013;許力生等，2013; 曾祥方等,2013).蘆山地震余震重定位結(jié)果顯示，余震主要分布在10～20 km深度上，淺部余震稀少，在垂直于龍門山斷裂帶的剖面上，余震集中在傾向相反、相互交叉的兩個(gè)條帶上，構(gòu)成“ Y ”字型分布(Fang et al., 2013; Long et al., 2015; 趙榮濤等, 2015; Lu et al., 2017).此外，深地震反射剖面進(jìn)一步顯示了蘆山地震區(qū)復(fù)雜的地殼結(jié)構(gòu)和斷裂構(gòu)造特征，震區(qū)地殼厚度約40 km，上部地殼分層特征明顯，褶皺和斷裂構(gòu)造清晰.由于受到北西-南東向擠壓應(yīng)力的作用，上部地殼形成了一系列傾角上陡下緩的鏟形斷層.結(jié)合余震重定位和震源機(jī)制解等研究結(jié)果，推斷蘆山地震發(fā)震斷層位于新開(kāi)店斷層和大邑?cái)鄬又g，伴有反沖斷層，斷層系統(tǒng)呈“ Y ”字型分布，深部歸并到約16 km左右的滑脫面上，斷層具有一定埋深(王夫運(yùn)等, 2015; 馮楊洋等, 2016).綜合地表地質(zhì)調(diào)查、余震定位、震源機(jī)制解以及深地震反射剖面等結(jié)果，推斷蘆山地震可能是一次發(fā)生在復(fù)雜斷層系統(tǒng)中的盲逆破裂事件(徐錫偉等, 2013; 王夫運(yùn)等, 2015).

這次地震的震源破裂過(guò)程是受到廣泛關(guān)注的問(wèn)題.不同的研究者利用地震波、GPS和InSAR等觀測(cè)資料，采用有限斷層模型對(duì)蘆山地震的破裂過(guò)程和同震滑動(dòng)分布進(jìn)行了反演研究(Hao et al., 2013; Jiang et al., 2013; Liu et al., 2013; 王衛(wèi)民等, 2013; 張勇等, 2013; Zhao et al., 2013; Zhang et al., 2014; 劉琦等, 2016; 鄭緒君等, 2018).結(jié)果顯示，此次地震以單個(gè)破裂事件為主，破裂方向性不明顯；破裂持續(xù)約10 s；同震滑動(dòng)主要集中在震源附近，分布在斷層走向方向約20 km范圍內(nèi)，最大滑動(dòng)量約為1.3 m，斷層淺部沒(méi)有明顯滑動(dòng)量分布，破裂沒(méi)有突破地表；整體破裂方向接近90°，此次地震可視為純逆沖型地震；地震震級(jí)約MW6.6；部分結(jié)果顯示震源東北方向最終滑動(dòng)量較小(Hao et al., 2013; Zhang et al., 2014; 鄭緒君等, 2018).而GPS、InSAR等反演手段可能受發(fā)震斷層模型影響較大，得到的平均最大滑動(dòng)量約為0.8 m(Jiang et al., 2013; 譚凱等, 2015; 劉琦等, 2016).

目前，針對(duì)蘆山地震破裂過(guò)程和同震滑動(dòng)分布的研究多通過(guò)波形擬合等反演手段進(jìn)行，沒(méi)有考慮導(dǎo)致破裂過(guò)程出現(xiàn)的原因——應(yīng)力因素.為了從力學(xué)角度解釋此次地震震源破裂過(guò)程的原因，需要進(jìn)行震源動(dòng)力學(xué)研究.本文基于蘆山地震野外地質(zhì)調(diào)查、震源機(jī)制解反演、余震定位、深地震反射剖面等結(jié)果，構(gòu)建蘆山地震發(fā)震斷層模型，在此基礎(chǔ)上，結(jié)合巖石力學(xué)實(shí)驗(yàn)得到的結(jié)果，以震源運(yùn)動(dòng)學(xué)反演結(jié)果作為約束，通過(guò)試錯(cuò)法調(diào)整震源動(dòng)力學(xué)計(jì)算參數(shù)，對(duì)蘆山地震進(jìn)行大量動(dòng)力學(xué)模擬研究，得到蘆山地震破裂傳播的可能情況，進(jìn)而分析討論不同動(dòng)力學(xué)計(jì)算參數(shù)對(duì)蘆山地震破裂過(guò)程和同震滑動(dòng)分布的影響.

1 三維斷層模型的構(gòu)建

目前對(duì)于蘆山地震破裂過(guò)程和同震滑動(dòng)分布的研究，大多將發(fā)震斷層簡(jiǎn)化為具有固定傾角的單一平面(Hao et al., 2013; Liu et al., 2013; 王衛(wèi)民等, 2013; Zhang et al., 2014; 鄭緒君等, 2018).但是蘆山地震余震定位、深地震反射剖面等結(jié)果表明，此次地震實(shí)際的發(fā)震斷層可能是一個(gè)上陡下緩的鏟型斷層，淺部連接近水平的滑脫帶，幾何結(jié)構(gòu)比較復(fù)雜(Fang et al., 2013; Long et al., 2015; 王夫運(yùn)等, 2015; 馮楊洋等, 2016; Lu et al., 2017).由于斷層的幾何結(jié)構(gòu)對(duì)地震破裂傳播有著重要影響(Aochi et al., 2000a, b; Aochi and Fukuyama, 2002)，因此本文結(jié)合蘆山地震余震定位、深地震反射剖面等結(jié)果，構(gòu)建蘆山地震三維鏟型斷層模型，該模型包括NW傾向的主斷層及其連接的深度約16.6 km的基底滑脫帶，如圖1所示.為了更好地刻畫鏟型，本研究采用非結(jié)構(gòu)化的三角形網(wǎng)格離散方案(錢峰等, 2019)，在計(jì)算量允許的情況下對(duì)模型進(jìn)行離散，如圖1所示.斷層模型走向方向長(zhǎng)30 km，傾向方向長(zhǎng)22 km，共1640個(gè)三角形單元，三角形網(wǎng)格的尺度Δs=0.245 km.

2 研究方法和參數(shù)設(shè)置

本研究采用能夠靈活處理非平面斷層破裂問(wèn)題的邊界積分方程法(BIEM)，對(duì)蘆山地震破裂過(guò)程進(jìn)行數(shù)值模擬.由于蘆山地震發(fā)震斷層上斷點(diǎn)埋深大于1 km，自由表面對(duì)破裂過(guò)程的影響可以忽略不計(jì)(Zhang and Chen, 2006)，因此可以采用基于全空間模型的BIEM.這種方法利用基本解建立相應(yīng)的邊界積分方程，將域問(wèn)題轉(zhuǎn)化為邊界問(wèn)題，因此不僅降低了求解維度，而且便于求解具有復(fù)雜幾何形狀的斷層破裂問(wèn)題，比如彎折斷層、鏟型斷層等.由于僅與斷層幾何形狀有關(guān)的積分核計(jì)算與破裂過(guò)程的計(jì)算是分離的，因此非常適合幾何形態(tài)固定，需要通過(guò)試錯(cuò)法探究動(dòng)力學(xué)參數(shù)對(duì)破裂過(guò)程影響的研究.

2.1 研究方法

BIEM基于格林函數(shù)求解彈性波動(dòng)方程，建立斷層上的牽引力與滑動(dòng)速率函數(shù)之間的關(guān)系(Fukuyama and Madariaga, 1998; Tada et al., 2000)：

(1)

圖1 蘆山地震斷層模型 x表示斷層傾向方向，y表示斷層走向方向，z表示深度.Fig.1 The fault model of the Lushan earthquake x is the direction of fault inclination, y is the direction of fault strike, and z is the depth.

為了求解自發(fā)破裂問(wèn)題，需要結(jié)合一定的摩擦本構(gòu)關(guān)系.摩擦本構(gòu)關(guān)系決定了破裂行為，控制了破裂過(guò)程.在同震破裂計(jì)算問(wèn)題中，通常采用滑動(dòng)弱化準(zhǔn)則(Ida, 1972)，其表示形式為：

(2)

其中，T(D)表示滑動(dòng)量為D時(shí)的剪切應(yīng)力，Tu=σu-σf表示應(yīng)力降，σu表示破裂強(qiáng)度，σf表示剩余應(yīng)力，Dc表示臨界滑動(dòng)弱化位移，H(Dc-D)表示階躍函數(shù)，當(dāng)Dc>D時(shí)值為1.對(duì)于斷層面上的某一點(diǎn)，在構(gòu)造應(yīng)力場(chǎng)作用下，應(yīng)力逐漸積累，直到剪應(yīng)力值超過(guò)σu發(fā)生滑動(dòng).隨后剪應(yīng)力隨著滑動(dòng)量的增加而線性地減?。划?dāng)D超過(guò)Dc后，剪應(yīng)力保持為恒定的常數(shù)σf.

結(jié)合邊界積分方程和滑動(dòng)弱化準(zhǔn)則，定義初始應(yīng)力分布、臨界滑動(dòng)位移分布、成核區(qū)分布、破裂強(qiáng)度分布等動(dòng)力學(xué)參數(shù)后就可以進(jìn)行破裂過(guò)程模擬.

2.2 震源動(dòng)力學(xué)參數(shù)的設(shè)置

對(duì)動(dòng)力學(xué)破裂過(guò)程計(jì)算的結(jié)果依賴于多個(gè)動(dòng)力學(xué)參數(shù)的共同作用，因此有必要謹(jǐn)慎設(shè)置每一個(gè)計(jì)算參數(shù).

首先是區(qū)域初始應(yīng)力場(chǎng).由于蘆山地震可以視為純逆沖型地震，因此本文不考慮沿?cái)鄬幼呦蚍较虻膽?yīng)力.雖然很少有研究直接給出地下深部應(yīng)力情況，但是由于溫度和圍壓的變化，應(yīng)力在很大程度上與深度有關(guān)(如, Aochi and Madariaga, 2003; Zhang et al., 2019)，因此本文考慮隨深度變化的應(yīng)力設(shè)置.為了使破裂能夠在震源附近成核傳播，需要優(yōu)先考慮震源所在深度的應(yīng)力分布情況. 假設(shè)參考平面為同一深度下的水平面，其受到豎直方向應(yīng)力σV和水平方向應(yīng)力σH、τH，水平方向應(yīng)力垂直于斷層走向方向，沿著x方向?yàn)檎较?根據(jù)子斷層單元傾角的大小，通過(guò)應(yīng)力變換可以計(jì)算出子斷層單元上的初始正應(yīng)力σ與剪應(yīng)力τ.由于蘆山地震動(dòng)力學(xué)模擬受運(yùn)動(dòng)學(xué)反演結(jié)果約束，因此當(dāng)進(jìn)行試錯(cuò)實(shí)驗(yàn)時(shí)得到的最大滑動(dòng)量與運(yùn)動(dòng)學(xué)反演結(jié)果基本符合時(shí)，即可得到合適的區(qū)域初始應(yīng)力設(shè)置.本文選取的應(yīng)力設(shè)置如圖2a所示.

斷層的破裂強(qiáng)度一般是正應(yīng)力與靜摩擦系數(shù)的乘積.當(dāng)正應(yīng)力已知，設(shè)置靜摩擦系數(shù)即可計(jì)算得到破裂強(qiáng)度，本文中靜摩擦系數(shù)μs取0.6，動(dòng)摩擦系數(shù)μd取0.1(如, Aochi and Madariaga, 2003; Zhang et al., 2019).當(dāng)剪應(yīng)力超過(guò)破裂強(qiáng)度時(shí)，斷層開(kāi)始發(fā)生滑動(dòng).

對(duì)于臨界滑動(dòng)弱化位移Dc，大量的實(shí)驗(yàn)室實(shí)驗(yàn)以及野外觀測(cè)資料分析表明，在不同地震中Dc取值自由度較大.通常在脆性巖層中Dc為一常數(shù)，在地殼淺部和韌性巖層中Dc較大(Ohnaka，1992；張麗芬和姚運(yùn)生, 2013).在實(shí)際的地震動(dòng)力學(xué)模擬中，還應(yīng)根據(jù)具體情況不斷調(diào)整Dc取值，直至破裂與運(yùn)動(dòng)學(xué)反演結(jié)果大致符合.本文假設(shè)Dc隨深度變化(如, Aochi and Fukuyama, 2002; Zhang et al., 2019)，滿足公式(3) :

(3)

其中h1為4 km，h2為16 km，d0為0.3 m，見(jiàn)圖2b.

圖2 動(dòng)力學(xué)計(jì)算參數(shù)隨深度的變化 (a) 區(qū)域初始應(yīng)力:σV是豎直方向應(yīng)力，σH、τH為水平方向應(yīng)力; (b)臨界滑動(dòng)弱化位移Dc.Fig.2 Variation of dynamic parameters with depth (a) Regional initial stress:σV is vertical stress, σH and τH are horizontal stress; (b) Critical slip weakening displacement Dc.

h1、h2的取值是參考余震重定位時(shí)采用的速度模型設(shè)置的(Fang et al., 2013; 趙榮濤等, 2015).

為了使破裂在開(kāi)始時(shí)刻被觸發(fā)，在當(dāng)前網(wǎng)格條件下，成核區(qū)半徑小于3 km無(wú)法孕育破裂，因此設(shè)置一個(gè)半徑為R=3 km的成核區(qū)，對(duì)成核區(qū)內(nèi)部初始剪應(yīng)力Ti重新賦值，使其略大于破裂強(qiáng)度Tu，也就是假定成核區(qū)已發(fā)生破裂Ti=1.01Tu.

在動(dòng)力學(xué)破裂模擬中，若沒(méi)有其他障礙體的情況下，破裂永不停止.但是根據(jù)動(dòng)力學(xué)反演結(jié)果來(lái)看(如, Hao et al., 2013; Liu et al., 2013; Zhang et al., 2014; 鄭緒君等, 2018)，滑動(dòng)量主要分布在震源附近約20 km范圍內(nèi)，即破裂被限制在一定區(qū)域內(nèi).因此本研究假設(shè)震源20 km范圍以外存在一近圓形的障礙體，使得斷層的破裂強(qiáng)度是原來(lái)的2倍，當(dāng)破裂接觸到障礙體后，會(huì)自動(dòng)減慢甚至停止，以達(dá)到限制破裂區(qū)域的目的.

由于本研究主要分析初始應(yīng)力分布、臨界滑動(dòng)位移分布、成核區(qū)分布、破裂強(qiáng)度分布對(duì)破裂過(guò)程的影響，其余參數(shù)均采用固定值，為降低復(fù)雜度，本文采用均勻介質(zhì)模型，如表1所示，其中VP,VS,ρ分別為介質(zhì)的P波、S波速度以及密度, CFL為穩(wěn)定性系數(shù)(Fukuyama and Madariaga,1998)，Δs為計(jì)算步長(zhǎng)，Δt為時(shí)間步長(zhǎng)，Step為計(jì)算步數(shù)，H為震源深度.

表1 介質(zhì)參數(shù)與斷層模型參數(shù)設(shè)置Table 1 Parameters of medium and fault model

3 不同動(dòng)力學(xué)參數(shù)的破裂傳播過(guò)程模擬

為了分析初始應(yīng)力分布、臨界滑動(dòng)位移分布、成核區(qū)分布、破裂強(qiáng)度分布對(duì)破裂過(guò)程和最終滑動(dòng)量分布影響，本文采用控制變量法(張麗芬等，2016)進(jìn)行大量數(shù)值模擬實(shí)驗(yàn)，通過(guò)對(duì)比模擬的結(jié)果探究不同震源動(dòng)力學(xué)參數(shù)對(duì)蘆山地震破裂過(guò)程和最終滑動(dòng)量分布的影響.

3.1 初始應(yīng)力的影響

為了研究斷層在不同初始應(yīng)力狀態(tài)下的破裂情況，在不改變其他動(dòng)力學(xué)參數(shù)的情況下，選取三組不同的初始應(yīng)力設(shè)置，以震源深度15 km為例，圖3顯示了A、B和C三種應(yīng)力狀態(tài)在莫爾圓中的位置，橫坐標(biāo)代表斷層面上正應(yīng)力，縱坐標(biāo)代表斷層面上剪應(yīng)力，τp=μsσ代表破裂線.其中，A狀態(tài)代表圖2a區(qū)域應(yīng)力設(shè)置中震源位置的受力情況.A狀態(tài)下震源處所受正應(yīng)力為7.53 MPa，剪應(yīng)力為3.2 MPa，B狀態(tài)下震源處所受正應(yīng)力為7.71 MPa，剪應(yīng)力為3.59 MPa，C狀態(tài)下震源處所受正應(yīng)力為7.13 MPa，剪應(yīng)力為2.89 MPa.從圖3可以看出，C點(diǎn)距離破裂線最遠(yuǎn)，B點(diǎn)距離破裂線最近.

圖4展示了不同初始應(yīng)力狀態(tài)下的滑動(dòng)速率快照(a—c)和最終滑動(dòng)量分布圖(d—f).可以看到，在A、B應(yīng)力狀態(tài)下，破裂在成核區(qū)孕育，隨著剪應(yīng)力的逐漸積累，破裂沿?cái)鄬幼呦蚍较虺蕶E圓形向四周傳播，在接觸到障礙體后，破裂逐漸停止；而在C應(yīng)力狀態(tài)下，剪應(yīng)力的積累沒(méi)能克服斷層的破裂強(qiáng)度，破裂在此應(yīng)力狀態(tài)下無(wú)法穩(wěn)定地向外傳播.根據(jù)滑動(dòng)弱化摩擦準(zhǔn)則，當(dāng)剪應(yīng)力超過(guò)斷層的破裂強(qiáng)度時(shí)，斷層兩側(cè)才開(kāi)始相對(duì)移動(dòng)，由此可知，初始應(yīng)力是決定斷層是否發(fā)生錯(cuò)動(dòng)的關(guān)鍵因素之一.

相比較而言，在A應(yīng)力狀態(tài)下破裂傳播較慢，沿走向方向約8 s到達(dá)障礙體產(chǎn)生停止震相，最大滑動(dòng)速率約0.36 m·s-1，沿?cái)鄬觾A向方向破裂未能傳播到斷層頂部，最終滑動(dòng)量分布范圍比較小，最大滑動(dòng)量約為1.1 m；而在B應(yīng)力狀態(tài)下，破裂快速向外傳播，沿走向方向約6 s到達(dá)障礙體產(chǎn)生停止震相，最大滑動(dòng)速率約0.68 m·s-1，沿?cái)鄬觾A向方向破裂傳播接近斷層頂部，最終滑動(dòng)量分布范圍較大，最大滑動(dòng)量約為1.37 m.這意味著，初始應(yīng)力設(shè)置越接近破裂線，破裂傳播越快，最大滑動(dòng)速率越大；同時(shí)，破裂沿著傾向方向的“爬坡能力”更強(qiáng)，最終滑動(dòng)量分布范圍越大，最大滑動(dòng)量數(shù)值越大.

圖3 初始應(yīng)力分布 圓表示應(yīng)力莫爾圓，τp=μsσ和τr=μdσ分別表示破裂強(qiáng)度線和殘余應(yīng)力線，A、B和C表示三種不同應(yīng)力狀態(tài)，σ、τ分別表示正應(yīng)力 和剪應(yīng)力，μs、μd分別表示靜摩擦系數(shù)和動(dòng)摩擦系數(shù).Fig.3 Initial stress distribution Circle is the Mohr circle of stress. τp=μsσ and τr=μdσ are the breaking strength line and the residual stress line, respectively. A, B and C represent three different stress states. σ and τ are the normal and shear stresses, μs and μd are the static and kinematic friction coefficients, respectively.

圖4 不同初始應(yīng)力狀態(tài)下的滑動(dòng)速率快照?qǐng)D(a—c)和最終滑動(dòng)量分布圖(d—f)Fig.4 Snapshots of slip rate (a—c) and final slip distribution (d—f) under different initial stress states

根據(jù)Hao等(2013)對(duì)蘆山地震運(yùn)動(dòng)學(xué)反演結(jié)果，破裂持續(xù)時(shí)間約10 s左右，最大滑動(dòng)量約1.2 m，淺于8 km的斷層面上滑動(dòng)量分布接近于0，破裂并未表現(xiàn)明顯的“爬坡能力”，因此在當(dāng)前參數(shù)設(shè)置下，A狀態(tài)下的應(yīng)力取值較為合適.

3.2 Dc的影響

臨界滑動(dòng)弱化位移Dc是摩擦準(zhǔn)則中的重要參數(shù).由于本研究采用分段函數(shù)表示Dc隨深度變化分布，如圖2b，因此設(shè)置不同的d0即可代表Dc對(duì)破裂情況的影響，本文選取d0=0.3 m、d0=0.2 m和d0=0.4 m三組情況值進(jìn)行破裂模擬.

圖5顯示了d0=0.3 m和d0=0.2 m的破裂滑動(dòng)速率快照?qǐng)D(a,b)和最終滑動(dòng)量分布圖(c,d).結(jié)果顯示，當(dāng)d0=0.3 m時(shí)，破裂傳播較慢，最大滑動(dòng)速率約0.36 m·s-1，破裂前鋒比較寬，破裂傳播穩(wěn)定，沿?cái)鄬觾A向方向破裂未能傳播到斷層頂部，最終滑動(dòng)量分布范圍比較小，最大滑動(dòng)量約為1.1 m；當(dāng)d0=0.2 m時(shí)，破裂傳播較快，最大滑動(dòng)速率約1.22 m·s-1，破裂前鋒非常窄，沿?cái)鄬觾A向方向破裂傳播接近斷層頂部，最終滑動(dòng)量分布范圍較大，最大滑動(dòng)量約為1.45 m.然而，在當(dāng)前動(dòng)力學(xué)參數(shù)設(shè)置下，當(dāng)d0=0.4 m時(shí)，破裂在成核區(qū)孕育失敗，無(wú)法向外傳播.這意味著，Dc對(duì)破裂滑動(dòng)速率有著很大的影響.Dc越小，破裂越容易發(fā)生，破裂傳播越快，最大滑動(dòng)速率越大；同時(shí)，破裂沿著傾向方向的“爬坡能力”更強(qiáng)，最終滑動(dòng)量分布范圍越大，最大滑動(dòng)量數(shù)值越大；但是相對(duì)的破裂前鋒比較窄，在數(shù)值模擬中容易積累誤差.

當(dāng)給定斷層面應(yīng)力情況后，在運(yùn)動(dòng)學(xué)反演最大滑動(dòng)量結(jié)果的約束下，Dc就只能在較窄的范圍內(nèi)取值.根據(jù)鄭緒君等(2018)反演得到的蘆山地震滑動(dòng)速率最大為0.48 m·s-1可知，蘆山地震破裂滑動(dòng)速率并未超過(guò)0.5 m·s-1，因此在當(dāng)前參數(shù)設(shè)置下，取d0=0.3 m較為合適.

3.3 成核區(qū)半徑的影響

由于破裂主要發(fā)生在震源附近20 km范圍內(nèi)，因此在設(shè)置成核區(qū)的時(shí)候不宜太大.為了研究成核區(qū)的大小對(duì)破裂情況的影響，本文選取成核區(qū)半徑R=3 km、R=3.5 km和R=2.5 km三組情況值進(jìn)行破裂模擬.

圖6顯示了R=3 km和R=3.5 km的破裂滑動(dòng)速率快照?qǐng)D(a,b)和最終滑動(dòng)量分布圖(c,d).結(jié)果顯示，當(dāng)R=3 km時(shí)，破裂傳播較慢，在6 s時(shí)還未到達(dá)障礙體，最大滑動(dòng)速率約0.36 m·s-1，最大滑動(dòng)量約為1.1 m；當(dāng)R=3.5 km時(shí)，破裂傳播較快，最大滑動(dòng)速率約0.5 m·s-1，最大滑動(dòng)量約為1.45 m；兩者的最終滑動(dòng)量分布范圍相差不大.然而，在當(dāng)前動(dòng)力學(xué)參數(shù)設(shè)置和網(wǎng)格條件下，R=2.5 km時(shí)破裂在成核區(qū)孕育失敗，無(wú)法向外傳播.這意味著，在固定的網(wǎng)格劃分情況下，成核區(qū)需要足夠大，蘊(yùn)含足夠多的能量，才能夠使破裂成功向外傳播；同時(shí)，成核區(qū)半徑大小更多影響的是破裂成核的快慢.成核區(qū)半徑越大，應(yīng)力積累的越快，裂紋向外傳播的越早，成核區(qū)最終滑動(dòng)量越大.

圖5 不同d0的滑動(dòng)速率快照?qǐng)D(a,b)和最終滑動(dòng)量分布圖(c,d)Fig.5 Snapshots of slip rate (a,b) and final slip distribution (c,d) under different d0

圖6 不同成核區(qū)半徑的滑動(dòng)速率快照?qǐng)D(a,b)和最終滑動(dòng)量分布圖(c,d)Fig.6 Snapshots of slip rate (a,b) and final slip distribution (c,d) under different radius of nucleation zone

根據(jù)Zhao等(2013)對(duì)蘆山地震運(yùn)動(dòng)學(xué)反演結(jié)果，蘆山地震在前2 s并未見(jiàn)明顯破裂，直到4 s左右破裂才明顯向外傳播，因此在當(dāng)前參數(shù)設(shè)置下，取R=3 km較為合適.

3.4 成核區(qū)初始剪應(yīng)力的影響

為了使破裂能被觸發(fā)，在設(shè)置成核區(qū)初始剪應(yīng)力Ti時(shí)，通常使其略高于斷層的破裂強(qiáng)度Tu.為了研究成核區(qū)的初始剪應(yīng)力大小對(duì)破裂情況的影響，本文選取Ti=1.01Tu、Ti=1.001Tu兩組情況值進(jìn)行破裂模擬.

圖7顯示了不同成核區(qū)初始剪應(yīng)力的破裂滑動(dòng)速率快照?qǐng)D(a,b)和最終滑動(dòng)量分布圖(c,d)，結(jié)果顯示，當(dāng)Ti=1.01Tu時(shí)，破裂在3.96 s時(shí)已經(jīng)明顯地向外傳播，最大滑動(dòng)速率約為0.36 m·s-1,最大滑動(dòng)量為1.11 m；當(dāng)Ti=1.001Tu時(shí)，3.96 s時(shí)破裂仍舊在成核區(qū)積累能量，在5.96 s時(shí)明顯地向外傳播，最大滑動(dòng)速率約0.34 m·s-1，滑動(dòng)量分布范圍相對(duì)較小，最大滑動(dòng)量為0.96 m.這意味著，設(shè)置成核區(qū)初始剪應(yīng)力略大于破裂強(qiáng)度，初步觸發(fā)破裂，破裂起始后需要一定時(shí)間積聚能量才會(huì)顯著地向外傳播.當(dāng)Ti越大，成核區(qū)初始蘊(yùn)含的能量越大，需要積聚能量向外傳播的時(shí)間越短，破裂越容易孕育并向外傳播，最終滑動(dòng)量分布范圍就越大，最大滑動(dòng)量值越大；Ti越小，并不會(huì)特別影響破裂的滑動(dòng)速率.

根據(jù)Hao等(2013)、Zhang等(2014)對(duì)蘆山地震運(yùn)動(dòng)學(xué)反演結(jié)果，破裂持續(xù)時(shí)間約10 s左右，因此在當(dāng)前參數(shù)設(shè)置下，Ti=1.01Tu的取值較為合適.

3.5 局部不均勻破裂強(qiáng)度的影響

當(dāng)斷層局部地區(qū)存在破裂強(qiáng)度不均勻現(xiàn)象，可能導(dǎo)致破裂傳播情況發(fā)生變化.本文假設(shè)震源東北方向有一破裂強(qiáng)度為20 MPa的局部障礙體進(jìn)行破裂模擬.

圖8顯示了加入局部不均勻破裂強(qiáng)度的破裂滑動(dòng)速率快照?qǐng)D(a)和最終滑動(dòng)量分布圖(b).結(jié)果表明，在未傳播到局部障礙體的時(shí)候，破裂傳播一切正常；在5 s時(shí)破裂接觸到局部障礙體后，由于其破裂強(qiáng)度過(guò)大，破裂無(wú)法突破，從而產(chǎn)生停止震相，最終破裂沒(méi)能繼續(xù)向震源東北方向傳播；最終滑動(dòng)量分布范圍不再呈橢圓形分布，震源東北方向滑動(dòng)量很小.

圖7 不同成核區(qū)初始剪應(yīng)力的滑動(dòng)速率快照?qǐng)D(a,b)和最終滑動(dòng)量分布圖(c,d)Fig.7 Snapshots of slip rate (a,b) and final slip distribution (c,d) under different initial shear stress in nucleation zone

圖8 加入局部不均勻破裂強(qiáng)度的滑動(dòng)速率快照?qǐng)D(a)和最終滑動(dòng)量分布圖(b)Fig.8 Snapshot of slip rate (a) and final slip distribution (b) with introducing local inhomogeneous fracture strength

Hao等(2013)、Zhang等(2014)和鄭緒君等(2018)對(duì)蘆山地震運(yùn)動(dòng)學(xué)反演結(jié)果表明，震源東北方向滑動(dòng)量很小，也就是說(shuō)破裂在這個(gè)區(qū)域基本不發(fā)生傳播，斷層最終滑動(dòng)量分布呈不規(guī)則形狀，這與圖8結(jié)果一致，因此可能是斷層面上存在局部障礙體導(dǎo)致的.當(dāng)破裂傳播遇到無(wú)法突破的障礙體，就會(huì)導(dǎo)致破裂行為發(fā)生改變，最終影響斷層面上滑動(dòng)量分布情況.

4 結(jié)論

本文基于蘆山地震野外地質(zhì)調(diào)查、余震定位、深地震反射剖面等信息，構(gòu)建蘆山地震鏟型斷層模型，在震源運(yùn)動(dòng)學(xué)反演結(jié)果的約束下，利用BIEM計(jì)算蘆山地震動(dòng)力學(xué)破裂過(guò)程，討論不同動(dòng)力學(xué)計(jì)算參數(shù)對(duì)破裂過(guò)程和同震滑動(dòng)分布的影響.結(jié)果顯示，初始應(yīng)力是決定斷層是否發(fā)生錯(cuò)動(dòng)的關(guān)鍵，初始應(yīng)力設(shè)置越接近破裂線，破裂越容易孕育傳播，滑動(dòng)速率越大，滑動(dòng)量越大；Dc對(duì)破裂滑動(dòng)速率有著很大的影響，Dc越小，破裂越容易發(fā)生，滑動(dòng)速率越大，但是破裂前鋒比較窄，在數(shù)值模擬中容易積累誤差；成核區(qū)半徑和成核區(qū)初始應(yīng)力主要影響破裂成核的快慢，成核區(qū)半徑越大，應(yīng)力積累的越快，裂紋向外傳播的越早，成核區(qū)最終滑動(dòng)量越大；成核區(qū)初始應(yīng)力越大，成核區(qū)初始蘊(yùn)含的能量越大，需要積聚能量向外傳播的時(shí)間越短，破裂越容易孕育并向外傳播，最終滑動(dòng)量分布范圍就越大，最大滑動(dòng)量值越大；局部不均勻破裂強(qiáng)度主要影響破裂行為和斷層最終滑動(dòng)量分布，當(dāng)破裂傳播遇到局部障礙體，會(huì)導(dǎo)致破裂行為發(fā)生改變，導(dǎo)致斷層最終滑動(dòng)量分布呈不規(guī)則形狀.

本文的動(dòng)力學(xué)模擬結(jié)果可以很好地再現(xiàn)蘆山地震破裂特征，從震源動(dòng)力學(xué)角度解釋為何斷層面上滑動(dòng)量分布不規(guī)則，對(duì)于了解導(dǎo)致這次地震的力學(xué)原因有一定的參考價(jià)值.但是，由于實(shí)際地震的震源過(guò)程是多個(gè)因素共同作用的結(jié)果，在沒(méi)有其余參數(shù)的確切取值的情況下，很難給出某一個(gè)物理參數(shù)的確切值.這還有賴于未來(lái)的綜合多種手段的研究來(lái)進(jìn)一步約束.