基于斷層傾角分段變化的瑪多地震發(fā)震斷層構(gòu)造應力場演化數(shù)值模擬分析

吳立新， 盧菁琛， 毛文飛*， 胡俊， 周子龍， 李志偉， 齊源， 姚汝冰

1 中南大學地球科學與信息物理學院， 長沙 410083 2 中南大學地災感知認知預知實驗室， 長沙 410083 3 中南大學資源與安全工程學院， 長沙 410083

0 引言

2021年5月22日2時4分，青海省果洛藏族自治州瑪多縣發(fā)生了MS7.4強震，震中位于東昆侖斷裂帶以南約70 km(34.59°N，98.34°E)、鮮水河斷裂的北側(cè)，震源深度17 km(中國地震臺網(wǎng)中心，2021).震后大量研究人員開展了野外地質(zhì)調(diào)查研究(潘家偉等，2021；李智敏等，2021)、地球物理反演(王未來等，2021；尹欣欣等，2021；詹艷等，2021；楊君妍等，2021)、衛(wèi)星遙感分析(祝愛玉等，2021；李志才等，2021；劉雷等, 2021; Qi et al., 2021；Wang et al.,2022).現(xiàn)場地質(zhì)調(diào)查資料顯示，瑪多地震形成了包括張剪裂隙、擠壓鼓包等多種破壞類型的復雜地表破裂帶，為典型左行走滑位錯性質(zhì).該破裂帶主要沿東昆侖斷裂帶南部的江錯斷裂分布，整體呈N105°E走向，全長約151 km，最大水平位錯量量約2.9 m，位于江錯斷裂西段末端(潘家偉等，2021).

瑪多地震發(fā)生在青藏高原中部巴顏喀拉塊體上，該區(qū)域是中國大陸的地震活躍區(qū)，其周緣孕育了包括東昆侖斷裂、鮮水河斷裂、甘孜—玉樹斷裂、馬爾蓋查卡—若拉崗日斷裂、龍門山斷裂和阿爾金斷裂在內(nèi)的多條活躍斷裂(張培震等，2003；鄧起東等，2002，2010).在印度板塊的擠壓作用下，巴顏喀拉塊體呈NE向運動，且其滑動速率自西向東逐漸減小(Kirby et al.，2007；李陳俠等，2011；李煜航等，2015).近二十多年來，在巴顏喀拉地塊周緣已經(jīng)發(fā)生了9次7.0級以上的大地震，如2008年四川汶川MS8.0地震(徐錫偉等，2008b)、2010年青海玉樹MS7.1地震(陳立春等，2010).

瑪多地震的發(fā)震斷層為江錯斷裂，是一條與達日斷裂、瑪多—甘德斷裂和東昆侖斷裂近平行展布的深大斷裂(潘家偉等，2021，王未來等，2021)，總長約700 km，整體走向為NWW，斷層傾向以NE為主，傾角近垂直，為左旋走滑性質(zhì)(王未來等，2021).從活動斷裂劃分上來看，江錯斷裂屬于東昆侖斷裂南側(cè)的次級斷裂(鄧起東等，2002；徐錫偉等，2008a)，瑪多地震震中大體位于地表破裂帶的中段，且在破裂帶兩端均出現(xiàn)分枝狀破裂現(xiàn)象(潘家偉等，2021).

據(jù)已有研究，地震序列間大多具有觸發(fā)和延緩效應(程佳等，2011；楊興悅等，2013；李平恩等，2019；劉雷等，2021)；而本次瑪多地震沿發(fā)震斷裂帶兩端形成了分叉式地表破裂現(xiàn)象，這在青藏高原的走滑地震中首次顯現(xiàn)(潘家偉等，2021)；而且，斷層傾角的空間展布變化相較前序地震更為復雜，呈現(xiàn)出一種凹凸交替的分段變化形態(tài)(徐志國等，2021；王未來等，2021).因此，有必要詳細分析瑪多地震發(fā)震斷層上構(gòu)造應力的演化特征及未來破裂發(fā)震的危險性.

目前針對區(qū)域構(gòu)造應力場的演化模擬研究中，陳連旺等(2008)模擬了在川滇地區(qū)主要活躍斷裂上發(fā)生的強震對潛在強震庫侖破裂應力的加卸載作用，發(fā)現(xiàn)活躍斷裂同震時引起的破裂庫侖應力變化大多呈增加趨勢.楊興悅等(2013)采用單元降剛法模擬了自1900年以來巴顏喀拉塊體及周緣地區(qū)發(fā)生的MS7.0以上的地震，計算了該地區(qū)先前地震對后續(xù)地震的觸發(fā)效應，結(jié)果表明：歷史地震的誘發(fā)效應并不是完全連續(xù)的，而是具有一定的時間滯后性.吳萍萍等(2014)針對1893年以來強活躍性的左行走滑斷裂——鮮水河斷裂帶發(fā)生的MS6.7以上地震，計算了各次地震對于后續(xù)地震的同震庫侖破裂應力的變化，以及庫侖應力的變化對后續(xù)地震的觸發(fā)作用.李平恩等(2019)以GNSS水平運動速率觀測值和最大主壓應力方向測量值為約束條件，獲得了構(gòu)造背景的初始應力場；進而模擬自1900年以來巴顏喀拉塊體周緣MS7.0級以上強震序列，研究了強震與庫侖破裂應力變化的關系；結(jié)果表明：巴顏喀拉活躍塊體上的強震對后續(xù)地震的發(fā)生既有觸發(fā)效應也有延緩作用.而針對2021MS7.4瑪多地震，祝愛玉等(2021)基于InSAR地表形變觀測結(jié)果開展了區(qū)域應力場數(shù)值模擬計算，分析認為發(fā)震斷層西南端和東北端應力場變化有利于未來余震的發(fā)生，且瑪多地震的發(fā)生對部分臨近斷裂帶上的應力分布具有顯著影響，增大了其未來發(fā)生滑動的可能性.馮淦等(2021)基于地震破裂模型和均勻彈性半空間模型，計算和分析了瑪多地震對其周圍地區(qū)應力場和位移場的影響，發(fā)現(xiàn)瑪多地震的發(fā)生，增大了發(fā)震斷層西段末端以及東段末端西側(cè)區(qū)域的庫侖破裂應力.

鑒于模型復雜度會在很大程度上影響數(shù)值模擬計算效率以及結(jié)果的收斂性，且從尺度效應上來說適當簡化模型并不會對模擬結(jié)果產(chǎn)生過大的影響，過去在研究塊體邊界或某一斷裂帶上發(fā)生的地震時，大多采用簡化模型——將斷層處理為單一傾角斷層或直立斷層.但對于2021年瑪多地震而言，震后余震叢發(fā)，斷層空間的分段與傾角變化十分復雜，若依舊采用簡化模型進行模擬，則難以從應力演化的角度來解釋其余震位置、地表破裂的分段性特征(潘家偉等，2021；王未來等，2021).因此，本文基于瑪多地震發(fā)震斷層面結(jié)構(gòu)產(chǎn)狀的分段特征(王未來等，2021)，精細化重構(gòu)了發(fā)震斷層面模型，以此對瑪多地震發(fā)震斷層面上構(gòu)造應力場的時空演化過程進行了數(shù)值模擬分析.

1 三維走滑斷層有限元模型

1.1 模型構(gòu)建

為模擬瑪多MS7.4左行走滑地震的應力演變過程，根據(jù)已有的地質(zhì)資料以及震后中國地震局臺網(wǎng)提供的信息，構(gòu)建了33°51′N—35°30′N，96°41′E—100°33′E區(qū)域范圍內(nèi)的有限元模型.模型覆蓋范圍內(nèi)的構(gòu)造斷裂包括：達日斷裂、江錯斷裂、甘德南緣斷裂、西藏大溝—昌麻河斷裂，如圖1所示.

震后基于余震序列的精確定位研究發(fā)現(xiàn)(王未來等，2021；徐志國等，2021)，江錯斷層(發(fā)震斷層)并不是整體上單調(diào)南傾或北傾，其走向和傾角均呈現(xiàn)分段變化的特征，故該斷層在空間上呈現(xiàn)出凹凸體形態(tài).因此，在構(gòu)建江錯斷層模型時，根據(jù)地震序列重定位所獲得的各分段走向及傾角結(jié)果(王未來等，2021)，對江錯斷層的走向和傾角進行了細化處理.如圖2a所示，江錯斷層的走向可分為13段，由西向東依次為AB段、BC段、…、MN段，各分段走向介于280°～298°之間，而斷層整體平均走向為285°(王未來等，2021).相應地，發(fā)震斷層各分段具有不同的傾向和傾角(圖2b)，斷層中間段區(qū)域(F→K)整體為北傾，傾向為11°～27°，傾角為76.6°～86.4°，東、西兩段(A→F和K→N)則為南傾，傾向為190°～208°，傾角為76.6°～89.3°；為保證不同分段間的連續(xù)性，在構(gòu)建幾何模型時相鄰分段連接處的斷層傾角進行過渡式緩變處理.考慮到本文主要是針對發(fā)震斷層(江錯斷層)的模擬分析，圖2c所示區(qū)域內(nèi)除發(fā)震斷層外，其余各斷層均簡化為直立斷層.基于上述考慮，構(gòu)建如圖2c所示有限元網(wǎng)格模型.模型中各斷層幾何參數(shù)如表1所示.模型中X軸為東西向，Y軸為南北向，Z軸為深度變化方向；模型WE長233 km，NS長100 km，深度20 km；經(jīng)自適應網(wǎng)格劃分后的節(jié)點數(shù)為135886個，單元數(shù)為37191個.采用有限元計算軟件ANSYS對數(shù)值模型進行計算.

圖2 三維有限元模型及斷裂分布特征(a) 發(fā)震斷層各分段平均走向，正北走向為0°； (b) 發(fā)震斷層各分段平均傾向及傾角，如NE12°∠87.2°表示該段斷層傾向為北東12°，傾角為87.2°； (c) 有限元模型及斷裂分布(紅色為發(fā)震斷裂，黑色為非發(fā)震斷裂，虛線框范圍為本次精細化傾角區(qū)域). Fig.2 Three dimentional finite element model and faults distribution(a) Average strike of each fault segmentation, where north is 0°; (b) The average dip and dip angle of each fault segmentation,where NE12°∠87.2° illustrates that the fault segmentation dips to northeast (12°) with a dip angle 87.2°; (c) Finite element model and faults distribution (red represents seismogenic fault, black stands for non-seismogenic fault, and the dashed frame is sophisticated dip angle zone).

表1 模型幾何及介質(zhì)參數(shù)Table 1 Geometric and physical properties of model

1.2 本構(gòu)關系及材料參數(shù)

從構(gòu)造應力長期演變來看，地殼介質(zhì)應當是黏性的，體現(xiàn)出在構(gòu)造運動長期流變環(huán)境下的應力積累特性；從短期來看地殼介質(zhì)又是彈性的.為了反映這種長期為黏性、短期為彈性的變形特征，本研究采用了黏彈性本構(gòu)關系.根據(jù)研究區(qū)域的巖性分析結(jié)果(劉大明等，2018；白國典等，2018)，選取了模型主體介質(zhì)參數(shù)的密度為2700 kg·m-3，彈性模量為1800 MPa，泊松比為0.28.同時參考已有圍繞巴顏喀拉塊體地震模擬中所采用的物性參數(shù)(李平恩等，2019)，選取介質(zhì)的黏滯系數(shù)為1023Pa·s.考慮到斷層介質(zhì)的性質(zhì)要弱于周圍活動地塊介質(zhì)，故參照前人研究方法將斷層考慮為具有一定厚度的軟弱結(jié)構(gòu)帶，即主體介質(zhì)的楊氏模量參數(shù)乘以弱化系數(shù)，以體現(xiàn)斷層與周圍塊體間的差異，在本研究中所采用的弱化系數(shù)為1/3(楊興悅等，2013；李平恩等，2019).因江錯斷裂南側(cè)的斷層活動程度遠小于江錯斷層，故本研究將南側(cè)甘德南緣斷裂和達日斷裂的強度進行適度提升、減弱其活躍性，以便進行相似性模擬.各斷層的介質(zhì)參數(shù)具體見表1.

1.3 邊界條件及模型合理性驗證

本文采用GNSS年均速度場(Zheng et al., 2017)作為有限元模型的邊界約束，對GNSS觀測值進行插值從而獲得邊界上的速度場.由于模型的實際地理范圍為33°51′N—35°30′N，96°41′E—100°33′E，區(qū)內(nèi)的GNSS觀測站點數(shù)量有限(僅有9個)，若直接通過這些站點獲得的運動速度值進行內(nèi)插，會影響插值結(jié)果的可靠性，故將測站的范圍適當放大至32°N—37°N，96°E—101°E(均在巴顏喀拉塊體內(nèi)部).通過對空間范圍放大后的區(qū)域內(nèi)測站數(shù)據(jù)(共35個測站)進行內(nèi)插，獲得原區(qū)域邊界上的年均運動速度值，并與加載時間(10萬年)相乘，從而獲得模型邊界上總的位移加載量.由于在深度方向上GNSS觀測值的差異暫無法獲知，因此本研究設定：在模型的不同深度下運動速度保持一致.此外，考慮到研究區(qū)域東邊界受龍門山斷裂帶(近南北向)影響，其運動速率遠小于其他三個邊界，故將模型東邊界視為固定邊界.同時模型底部水平向自由，垂向固定，上表面為自由邊界.

為了檢驗模型參數(shù)的合理性，在上述模型幾何參數(shù)、物性參數(shù)及邊界約束條件下進行加載模擬，對比分析模擬得到的速度場與GNSS觀測獲得的速度場.結(jié)果如圖3所示，靠近模型四周邊界上的模擬值和觀測值在大小和方向上偏離程度均較?。欢Ｐ椭胁繀^(qū)域，二者在大小上雖然在局部區(qū)域存在一定差異(受斷層阻擋導致其運動速率減小)，但從整體的運動方向上來看，模擬結(jié)果很好地體現(xiàn)了受印度板塊推擠，巴顏喀拉塊體內(nèi)部區(qū)域呈自西向東運動的特征.

圖3 GNSS觀測值與數(shù)值模擬值的對比Fig.3 Comprison between GNSS observations and numerical simulations

2 模擬結(jié)果

2.1 庫侖破裂應力變化速率

根據(jù)庫侖破裂準則，當巖石趨近破裂時庫侖破裂應力σCFS為

σCFS=τ+μ(σn-p)，

(1)

式中μ為摩擦系數(shù)，τ為剪切應力(沿斷層面方向)，σn為正應力(壓應力為正)，p為孔隙壓力.地震模擬中通常采用庫侖破裂應力的變化來表征應力的積累速率(Harris，1998)，當μ為常數(shù)時，可得庫侖破裂應力的變化ΔσCFS為

ΔσCFS=Δτ+μ(Δσn-Δp)，

(2)

對于各向同性均勻介質(zhì)，孔隙壓力對摩擦系數(shù)的影響可用視摩擦系數(shù)μ′=μ(1-B)來表示，其中B為Skempton系數(shù)(取決于巖石膨脹系數(shù)和流體所占比例)，取值范圍為0～1(Rice，1992)，則

ΔσCFS=Δτ+μ′Δσn.

(3)

視摩擦系數(shù)μ′取值通常為0～0.8，0～0.4為大變形，0.4～0.8為小變形，平均值約為0.4(Freed，2005).視摩擦系數(shù)取值大小對于庫侖破裂應力的極性分布影響較小，僅對應力大小有一定的影響(萬永革等，2000)，故在本研究采用平均視摩擦系數(shù)0.4進行計算.考慮到研究區(qū)構(gòu)造運動的時間相似性(李平恩等，2019)，將獲取的三維模型邊界處運動速度場乘以10萬年的時間尺度，從而獲得模擬計算所需總位移加載量；為了便于分析應力演化過程，整個計算過程分為40個時間步進行位移加載.

由圖4可知，江錯斷層面上庫侖破裂應力年變化速率較大的位置有3處，分別為斷層西段末端的中上部區(qū)域、震源附近和斷層東段的末端區(qū)域(圖4).A、B、C三處的年庫侖破裂應力變化速率均為正值，即庫侖破裂應力呈增加趨勢，且各處最大庫侖破裂應力積累速率分別為717.5 Pa·a-1、388.6 Pa·a-1、2691.2 Pa·a-1，分別是江錯斷層面上整體平均庫侖應力積累速率(19.6 Pa·a-1)的36.6倍、19.8倍、137.3倍.

圖4 江錯斷層面上庫侖破裂應力的年變化速率Fig.4 Coulomb stress changing velocity per year on the Jiangcuo fault

從區(qū)域的差異性來看，震源以東斷層面的平均庫侖破裂應力積累速率為39.3 Pa·a-1，顯著大于震源以西斷層面的平均應力積累速率14.6 Pa·a-1.西段A區(qū)靠近鄂陵湖，深度范圍為0～10 km；其中包含有庫侖破裂應力變化速率為負值的子區(qū)域(圖4中藍色區(qū)域)，變化速率最大為-269.3 Pa·a-1，表明該負值子區(qū)域在構(gòu)造運動作用下應力逐漸釋放.中段震源附近(B區(qū)域)的庫侖破裂應力積累速率呈現(xiàn)以震源為中心向四周輻射狀減小的趨勢，其中震源位置的庫侖破裂應力積累速率為64.8 Pa·a-1.斷層東段末端C區(qū)是整個江錯斷層上應力積累速率最大的區(qū)域，應力積累集中分布在0～20 km深度范圍內(nèi)，同時在8 km、12 km和20 km深度處存在3個積累速率大于1000 Pa·a-1的子區(qū)域；在4 km、10 km和7 km深度處，則有3個庫侖破裂應力變化速率為負值的環(huán)狀區(qū)域，其最值為-244.7 Pa·a-1.本研究在震源附近的模擬計算結(jié)果與李平恩等(2019)計算所得的庫侖破裂應力積累速率(約為50 Pa·a-1)基本一致.

前人的地震模擬研究表明，庫侖破裂應力積累速率快與斷裂帶歷史強震復發(fā)周期短具有一致性，較低水平的應力積累速率對應強震發(fā)生的應力積累周期更長(劉雷等，2021).本次瑪多地震所在江錯斷層面上A、B、C處的庫侖破裂應力積累速率處于高水平狀態(tài)，即這些區(qū)域更易形成應力集中區(qū)，意味著這些區(qū)域發(fā)生地震(或強震)的可能性要遠高于其他區(qū)域，這或與瑪多主震、余震的空間分布特征具有一定的相關性.

2.2 江錯斷層面上的應力演化

(4)

式中σx、σy、σz、τxy、τyz、τzx為應力張量分量.等效應力增大會加快介質(zhì)失效而促進孕震過程，等效應力減小會阻礙介質(zhì)屈服而延緩孕震過程.模擬得到的江錯斷層面上的應力演化結(jié)果如圖5所示.

由圖5可知，江錯斷層面上的應力演化過程為：首先在江錯斷層面的東段末端4～20 km深度范圍內(nèi)形成最值為4 MPa的應力集中區(qū)域(圖5a)；隨后，在江錯斷層面西段末端(鄂陵湖南側(cè)區(qū)域)0～8 km深度范圍出現(xiàn)了約為2 MPa的傾斜帶狀應力集中區(qū)(圖5b)；之后，在震源附近出現(xiàn)了應力集中現(xiàn)象(圖5c，呈向東、西兩側(cè)擴展的趨勢)；此后，這3處的應力進一步積累、擴展和強化.當數(shù)值模擬的加載時間結(jié)束時，江錯斷層東段末端、江錯斷層西段末端和瑪多地震主震震源處的最大等效應力分別為296 MPa、109 MPa和40 MPa.

圖5 江錯斷層面上的應力演化結(jié)果(a) 加載時間為2500年； (b) 加載時間為7700年； (c) 加載時間為1.3萬年； (d) 加載時間為2.1萬年； (e) 加載時間為5.4萬年； (f) 對應的加載時間為10萬年(圓圈代表2021年5月22日至2021年11月1日內(nèi)在研究區(qū)發(fā)生的震級≥3.0余震序列，其中圓圈顏色代表震源深度，圓圈大小對應不同震級，五角星為主震)； (g) 主震、余震震源位置與江錯斷層面的空間關系(水平投影).Fig.5 Tectonic stress evolution on the Jiangcuo fault(a) 2,500 years after loading; (b) 7700 years after loading; (c) 13000 years after loading; (d) 21,000 years after loading; (e) 54000 years after loading; (f) 100000 years after loading (circles are earthquake series which greater than or equal to MS3.0, from May 22, 2021 to November 1, 2021, in which colors stand for the depth of earthquake and sizes represent earthquake magnitude, and star means main earthquake); (g) Horizontal projection of epicenters of main shock, aftershocks and Jiangcuo fault.

瑪多地震余震重定位的研究結(jié)果顯示，瑪多地震的地震序列呈現(xiàn)出明顯的空間分段特征，在發(fā)震斷層東段上存在著余震活動較少的稀疏帶(王未來等，2021；尹欣欣等，2021)，這與圖5g中江錯斷層面上東段的低應力區(qū)(K→L)較為相符.由此推測認為：斷層面上局部區(qū)域應力積累的空間差異導致了余震的空間分段性.結(jié)合震級大小和地震數(shù)量來看，發(fā)震斷層面上的地震序列基本可分為4段，其中有3處位于應力集中區(qū)，分別對應發(fā)震斷層面的東、西兩段末端及震源區(qū)域，有一處位于應力水平較低處(圖5f).整體來看，地震序列分布與發(fā)震斷層面上的應力積累情況具有較好的空間對應關系，積累應力較大位置對應的地震數(shù)量叢集、強度大.此外，發(fā)震斷層面西段與震中之間的地震序列密集段對應的是應力水平較低處；考慮到模型的構(gòu)建情況，東、西兩段以及震源處的應力集中區(qū)均位于斷層傾角變化的過渡處，在構(gòu)造運動自西向東推擠作用下，易于在傾角變化處積累應力，從而誘發(fā)地震；而低應力處位于兩分段的中間位置，并未出現(xiàn)明顯的應力積累，這也說明，本研究的模型構(gòu)建在分段劃分上還可進一步細化.

2.3 區(qū)域構(gòu)造應力演化

由圖6可知，地表構(gòu)造應力場演化過程為：整體構(gòu)造應力場的背景值為0.1～1 MPa；首先在江錯斷裂西段末端地表破裂分叉處出現(xiàn)了較小的應力集中區(qū)(圖6a)，大小約為1～2 MPa; 隨后斷層西段末端應力集中區(qū)向南北兩側(cè)擴展，且同時在江錯斷裂東段末端分叉處和西藏大溝—昌麻河斷裂也出現(xiàn)了約為1～2 MPa的應力集中區(qū)(圖6b).隨著構(gòu)造運動持續(xù)加載，東、西段末端區(qū)域的構(gòu)造應力進一步積累并向外擴展，其中在西段分叉處構(gòu)造應力增大至10～100 MPa(圖6c)，同時，在鄂陵湖區(qū)西南側(cè)也出現(xiàn)了局部的應力集中現(xiàn)象；此外，東、西段末端的應力集中現(xiàn)象，也開始沿著地表斷層線向震中方向擴展.當10萬年的加載過程結(jié)束時，斷層東、西段末端的應力集中區(qū)都進一步擴展并相互連通(圖6d)，整個斷層上最大應力集中區(qū)為西端分叉處; 震中東、西兩側(cè)斷層線上的構(gòu)造應力約1～2 MPa，但震中位置卻無顯著的構(gòu)造應力集中現(xiàn)象.

圖6 地表構(gòu)造應力場演化結(jié)果(a) 加載時間為3.1萬年； (b) 加載時間為5.1萬年； (c) 加載時間為8.2萬年； (d) 加載時間為10萬年.Fig.6 Evolution of tectonic stress on ground surface(a) 31000 years after loading; (b) 51000 years after loading; (c) 82000 years after loading; (d) 100000 years after loading.

作者在瑪多地震后的當年6月份開展了為期一周的震區(qū)野外科考，沿發(fā)震斷裂及環(huán)鄂陵湖追蹤調(diào)查了地表破裂與災情(圖7)；結(jié)合潘家偉等(2021)報告的地表破裂情況分析，發(fā)現(xiàn)本研究揭示的地表構(gòu)造應力集中區(qū)的分布情況與地表破裂帶位置基本一致.從整體的分布形態(tài)上來看，地表破裂帶在東、西兩端均形成了大規(guī)模分叉的“帚狀”結(jié)構(gòu)，但震中位置地表破裂發(fā)育并不突出；地表破裂由震中向東、西兩段逐漸發(fā)育并與東、西兩段末端連續(xù)；與之對應，圖5所示地表應力集中區(qū)分布情況也是在東西兩段的末端出現(xiàn)分叉，在震中并無明顯應力積累，且由震中位置向東、西兩段逐漸積累.此外，從局部上來看，西段的地表破裂北支向北西消失于鄂陵湖，南支向西延伸至鄂陵湖南側(cè)后逐漸消失，對應斷裂帶西段末端積累的應力，也是呈南、北兩支消失于鄂陵湖；而斷裂帶東段的地表破裂北支破裂程度明顯大于南支，這與東段末端的應力集中區(qū)分布差異也是一致的.

圖7 瑪多地震地表破裂現(xiàn)場調(diào)查結(jié)果(a—f分別對應圖6d中的1—6處)(a) 沙土液化； (b) 間斷式拉擠破裂； (c) 野馬灘大橋垮塌； (d) 地表拉張開裂； (e) 黃河鄉(xiāng)大橋橋墩破損； (f) 無人機探視.Fig.7 Surface ruptures investigation of Madoi earthquake (a—f correspond to 1—6 in Fig.6d respectively)(a) Liquefaction of sand layer; (b) Extrusion rupture on the ground surface; (c) Bridge on the Wild-horse wetland; (d) Tessive fissure on the ground surface; (e) Damaged bridge at the Yellow River town; (f) Unmanned aerial vehicle (UAV) visitation.

考慮到江錯斷層面的深度及其傾角分段性特征，計算震源深度處(17 km)水平面上的構(gòu)造應力場演化，結(jié)果如圖8所示.震源深度水平面內(nèi)的構(gòu)造應力集中現(xiàn)象主要出現(xiàn)在江錯斷層東段末端以及震源處.從應力演化過程來說，首先在東段末端與分支斷層的交匯處出現(xiàn)了約1～2 MPa的應力集中現(xiàn)象(圖8a).隨后，東段末端的應力集中區(qū)域呈南北向擴展態(tài)勢，其中心區(qū)域應力增大至10 MPa左右；同時，震源區(qū)域應力集中開始顯現(xiàn)(圖8b).當東段末端應力集中擴展至西藏大溝—昌麻河斷裂和甘德南緣斷裂時，震源區(qū)域應力集中沿斷層線向東、西兩側(cè)擴展(圖8c).當10萬年加載過程結(jié)束時，東段末端的應力集中區(qū)域已顯著影響到北側(cè)的西藏大溝—昌麻河斷裂和南側(cè)的甘德南緣斷裂，區(qū)域中心處最大應力增至100～300 MPa(圖8d); 同時，整個發(fā)震斷層都處于相對的高應力狀態(tài)，震源位置應力也增至10～100 MPa，可能是瑪多地震發(fā)震的根本原因.此外，從瑪多震源深度平面內(nèi)江錯斷層的走向分布來看，震源恰好位于兩分段的傾角變化處，因此在構(gòu)造運動的不斷加載下，易于在該位置形成應力閉鎖區(qū)；同時，整個巴顏喀拉塊體在周緣塊體的推擠下自西向東運動，這就導致震源處的應力沿斷層面向東、向西擴展，由此形成了一個條帶狀的應力集中區(qū).

圖8 震源所在深度水平面內(nèi)的構(gòu)造應力場演化結(jié)果(虛線為地表斷層線)(a) 加載時間為5000年； (b) 加載時間為1.6萬年； (c) 加載時間為3.6萬年； (d) 加載時間為10萬年.Fig.8 Evolution of tectonic stress in the horizontal plane at hypocenter depth (dashed line represents fault line on ground surface)(a) 5000 years after loading; (b) 16000 years after loading; (c) 36000 years after loading; (d) 100000 years after loading.

對比地表和震源所在平面的構(gòu)造應力時空演化過程(圖6和圖8)可知：首先，瑪多地震主震震源所在水平面內(nèi)應力集中現(xiàn)象主要體現(xiàn)在東段末端區(qū)域以及整個江錯斷層上，而地表的應力集中主要體現(xiàn)在東、西段末端區(qū)域(震中位置未見應力集中現(xiàn)象)，說明斷層走向及傾角的變化對整個斷層面上構(gòu)造應力的積累與集中發(fā)育過程以及三維空間展布具有顯著影響.

其次，在時間上，江錯斷層在瑪多地震震中區(qū)域及斷層東段末端區(qū)域應力演化要比同區(qū)域震源深度處的應力演化滯后，即當加載時間為3.6萬年時，在震源所在水平面上，震源位置斷層面和斷層東段均出現(xiàn)了明顯的應力集中現(xiàn)象(等效應力>6 MPa, 圖8c)，而當加載時間為5.1萬年時，地表才在江錯斷層東、西段的末端出現(xiàn)了一些應力集中現(xiàn)象(等效應力<6 MPa, 圖6b).

2.4 討論

前人針對巴顏喀拉塊體及周緣地區(qū)的地震模擬中，主要構(gòu)造斷層或發(fā)震斷層的產(chǎn)狀大多數(shù)呈單一傾角或直立傾角(楊興悅等，2013；李平恩等，2019；劉雷等，2021).瑪多地震發(fā)生在巴顏喀拉塊體內(nèi)部的歷史地震空白區(qū)，余震叢發(fā)且形成了分叉式地表破裂帶，采用單一傾角進行模擬難以解釋其孕震過程的特征性變化.為此，對圖2b所示模型進行簡化，進而對比分析斷層幾何模型的復雜度對計算結(jié)果的影響.

圖9a所示模型不考慮斷層走向及傾角的分段性特征，計算結(jié)果僅體現(xiàn)了構(gòu)造運動自西向東不斷推擠所導致的應力積累現(xiàn)象，斷層面上積累的最大等效應力為60 Pa.圖9b所示模型僅考慮斷層走向分段而無傾角變化(均為直立)，計算結(jié)果雖然體現(xiàn)了一定的分段性特征，但與圖5f所示計算結(jié)果相比，在深度方向上應力分布無顯著差異，斷層面上積累的最大等效應力僅為200 Pa.

圖9 簡化模型斷層面上的應力積累(a) 發(fā)震斷層直立且無走向分段； (b) 發(fā)震斷層直立且有與圖2a一致的走向分段.Fig.9 Accumulative stress on the fault of simplified model(a) The vertical seismogenic fault without strike sectioning; (b) The vertical seismogenic fault with strike sectioning as in Fig.2a.

顯然，在以上兩種模擬方式下，發(fā)震位置應出現(xiàn)在高等效應力的斷層東段側(cè)末端處，且不會在斷層兩段末端出現(xiàn)分叉式地表破裂帶.因此，精細的斷層模型方能更好地體現(xiàn)江錯斷裂東、西兩段的應力分布差異以及瑪多地震震源處的應力集中現(xiàn)象.綜合對比圖9與圖5f的計算結(jié)果，從斷層模型的應力積累大小及分布來看，對于瑪多地震這種震后依舊持續(xù)活躍且具有特征性地表破壞的情況，采用精細化的斷層模型進行有限元模擬分析是更為合適和科學合理的.

基于江錯斷層面上的庫侖破裂應力變化(圖4)及等效應力演化(圖5)、地表及瑪多地震主震震源深度水平面內(nèi)構(gòu)造應力場演化(圖6和圖8)、地表破裂發(fā)育(圖7)、主震及余震序列空間分布(圖6)及潘家偉等(2021)所做的地表破裂調(diào)查結(jié)果，綜合分析認為：1)受江錯斷裂走向及傾角的分段變化影響，瑪多地震震源區(qū)及斷層東、西段末端為構(gòu)造應力變化活躍區(qū)或高應力分布區(qū)，與瑪多地震主震及余震序列的空間分布相對應；2)地表破裂帶位置與地表構(gòu)造應力集中區(qū)具有空間相關性，體現(xiàn)了瑪多MS7.4主震會對地表高應力集中區(qū)域產(chǎn)生擾動，這也是導致地表分段式破裂的重要原因；3)顧及構(gòu)造應力、地表破裂、余震活動的空間分布差異，可進一步將沿江錯斷裂的瑪多地震孕震—發(fā)震影響劃分為5段，即由西向東分別為A-E段、E-G段、G-K段、K-L段、L-N段(分段標識號見圖2)：A-E段構(gòu)造應力主要在10 km以淺局部積累和變化，地表可見形變及破裂(潘家偉等，2021)，且地表高應力區(qū)向斷層南北兩側(cè)的擴展態(tài)勢明顯，余震在斷層兩側(cè)均有發(fā)育，尤其向南側(cè)延伸較遠；E-G段構(gòu)造應力在全深度上近似均勻分布和同步變化，地表形變與破裂發(fā)育顯著，是余震叢發(fā)分布區(qū)；G-K段構(gòu)造應力主要在深部(10～20 km)發(fā)育，應力閉鎖區(qū)靠近主震震源位置，而5 km以淺及地表應力積累很弱，地表形變及破裂發(fā)育也不明顯；K-L段與A-E段類似，應力積累也發(fā)生在10 km以淺，地表形變破裂顯著(潘家偉等，2021)，但因其構(gòu)造應力值小于A-E段，所以余震發(fā)育較少；L-N段應力積累作用明顯、在深度上近似均勻分布，地表形變破裂及余震活動均顯著發(fā)育，且在震源所在水平面上應力集中顯著并向斷層面南北兩側(cè)大范圍擴展.

需要注意的是：在震源深度水平面內(nèi)(圖8)，震源周圍的應力主要沿發(fā)震斷層面積累和發(fā)展，并未向斷層面南北兩側(cè)擴展，表明瑪多地震主震及余震的發(fā)生已對震源區(qū)斷層面上積累的應力有了較充分的釋放(He et al., 2021)，體現(xiàn)了震源走滑發(fā)震的應力演化特征.而江錯斷裂東段積累的應力向南、北兩側(cè)分別擴展至甘德南緣斷裂和西藏大溝—昌麻河斷裂，盡管這兩處斷裂吸收了部分構(gòu)造應力，但依舊越過這兩個斷裂繼續(xù)朝南、北向擴展(圖8d)，且其范圍遠遠大于余震序列的空間分布范圍(圖5g)，表明在發(fā)震斷層東段末端發(fā)生的余震并不能將江錯斷裂東段積累的應力全部釋放.受限于本次模擬設計模型的時空范圍，并未進一步獲得模型東南側(cè)的應力擴展情況.根據(jù)目前的應力擴展趨勢以及巴顏喀拉塊體內(nèi)斷裂的走向，我們推測：構(gòu)造應力具有向江錯斷裂東段的南部地區(qū)繼續(xù)積累的可能，巴顏喀拉塊體東南區(qū)域可能是未來發(fā)震的高風險區(qū).此外，考慮到地表構(gòu)造應力在江錯斷層西段末端的顯著積累與擴展(圖6)，結(jié)合該段斷層面上的應力積累變化(圖4)，推測未來在江錯斷層西段附近區(qū)域可能會發(fā)生淺源(深度<10 km)小型余震.本研究雖然未考慮應力演化過程中瑪多地震主震及余震對大區(qū)域構(gòu)造應力的擾動，但就未來潛在發(fā)震區(qū)的空間分析和推測而言，與祝愛玉等(2021)、馮淦等(2021)、He等(2021)分析結(jié)果相符.

3 結(jié)論

2021年瑪多MS7.4地震發(fā)生于巴顏喀拉塊體內(nèi)部的歷史地震空白區(qū), 震后余震叢發(fā)；地震所形成的地表破裂帶在震中區(qū)并不顯著，而在發(fā)震斷層東、西段兩端分叉式發(fā)育.本文基于地震序列精定位研究結(jié)果建立了斷層面的三維精細模型(精細化發(fā)震斷層的走向與傾角分段性變化)，對發(fā)震斷層面及巴顏喀拉塊體局部區(qū)域的構(gòu)造應力場演化過程進行了數(shù)值模擬計算，揭示了模擬區(qū)斷層及周邊應力積累的時空過程，表明震源應力時空變化特征與走滑地震發(fā)震機制吻合、地表等效應力空間分布特征與地表破裂發(fā)育位置相符.本文研究揭示的現(xiàn)象及分析結(jié)論包括：

(1)在板塊構(gòu)造應力作用下，江錯斷裂面東、西兩段末端區(qū)域和瑪多地震主震震源附近區(qū)域的庫侖破裂應力變化明顯，且等效應力積累作用顯著，與主震及余震序列具有較好的空間對應關系，表明江錯斷層面上的高應力積累是誘發(fā)瑪多走滑地震主震及余震的主要原因.

(2)從地表等效應力演化的時空過程來看，地表沿江錯斷層應力積累的空間分布與地表破裂的發(fā)育位置基本一致，均為東、西兩端較為顯著且在末端呈南北向分叉發(fā)展的態(tài)勢，而震中附近不甚發(fā)育，表明瑪多地震發(fā)震過程對地表高應力積累區(qū)的擾動是導致地表形變破裂顯著的主要原因.

(3)在震源深度處的水平面上，震源周圍的應力主要沿發(fā)震斷層面積累和發(fā)展，并未向斷層面南北兩側(cè)擴展，體現(xiàn)了震源走滑發(fā)震的應力演化特征；江錯斷層東段的構(gòu)造應力積累與擴展雖然受到了甘德南緣斷裂和西藏大溝—昌麻河斷裂的阻擋和部分吸收，但構(gòu)造應力積累區(qū)的空間范圍遠大于瑪多地震余震分布范圍，表明本次瑪多地震并沒有將斷層東段積累的構(gòu)造應力全部釋放.結(jié)合斷層東段應力向南、北兩側(cè)差異擴展的趨勢，推測巴顏喀拉塊體東南地區(qū)可能是未來發(fā)震高風險區(qū).

致謝震區(qū)野外科考工作中，得到青海大學史培軍校長、劉峰貴教授的指導和幫助；關于余震發(fā)育及地表破裂的分段性，與國家自然災害防治研究院的徐錫偉院長、申旭輝總工，中國地震局地質(zhì)研究所單新建所長及其團隊成員等進行了有益討論；關于巴顏喀拉塊體未來地震風險趨勢，與史培軍校長、徐錫偉院長進行了有益交流，在此一并表示感謝！同時感謝兩位匿名審稿專家提出的有益意見和建議.