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

1 中南大學(xué)地球科學(xué)與信息物理學(xué)院， 長沙 410083 2 中南大學(xué)地災(zāi)感知認(rèn)知預(yù)知實(shí)驗(yàn)室， 長沙 410083 3 中南大學(xué)資源與安全工程學(xué)院， 長沙 410083

0 引言

2021年5月22日2時(shí)4分，青海省果洛藏族自治州瑪多縣發(fā)生了MS7.4強(qiáng)震，震中位于東昆侖斷裂帶以南約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)查資料顯示，瑪多地震形成了包括張剪裂隙、擠壓鼓包等多種破壞類型的復(fù)雜地表破裂帶，為典型左行走滑位錯(cuò)性質(zhì).該破裂帶主要沿東昆侖斷裂帶南部的江錯(cuò)斷裂分布，整體呈N105°E走向，全長約151 km，最大水平位錯(cuò)量量約2.9 m，位于江錯(cuò)斷裂西段末端(潘家偉等，2021).

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

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

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

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

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

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

1.1 模型構(gòu)建

為模擬瑪多MS7.4左行走滑地震的應(yīng)力演變過程，根據(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)造斷裂包括：達(dá)日斷裂、江錯(cuò)斷裂、甘德南緣斷裂、西藏大溝—昌麻河斷裂，如圖1所示.

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

圖2 三維有限元模型及斷裂分布特征(a) 發(fā)震斷層各分段平均走向，正北走向?yàn)?°； (b) 發(fā)震斷層各分段平均傾向及傾角，如NE12°∠87.2°表示該段斷層傾向?yàn)楸睎|12°，傾角為87.2°； (c) 有限元模型及斷裂分布(紅色為發(fā)震斷裂，黑色為非發(fā)震斷裂，虛線框范圍為本次精細(xì)化傾角區(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)關(guān)系及材料參數(shù)

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

1.3 邊界條件及模型合理性驗(yàn)證

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

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

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

2 模擬結(jié)果

2.1 庫侖破裂應(yīng)力變化速率

根據(jù)庫侖破裂準(zhǔn)則，當(dāng)巖石趨近破裂時(shí)庫侖破裂應(yīng)力σCFS為

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

(1)

式中μ為摩擦系數(shù)，τ為剪切應(yīng)力(沿?cái)鄬用娣较?，σn為正應(yīng)力(壓應(yīng)力為正)，p為孔隙壓力.地震模擬中通常采用庫侖破裂應(yīng)力的變化來表征應(yīng)力的積累速率(Harris，1998)，當(dāng)μ為常數(shù)時(shí)，可得庫侖破裂應(yīng)力的變化Δσ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ù)取值大小對于庫侖破裂應(yīng)力的極性分布影響較小，僅對應(yīng)力大小有一定的影響(萬永革等，2000)，故在本研究采用平均視摩擦系數(shù)0.4進(jìn)行計(jì)算.考慮到研究區(qū)構(gòu)造運(yùn)動(dòng)的時(shí)間相似性(李平恩等，2019)，將獲取的三維模型邊界處運(yùn)動(dòng)速度場乘以10萬年的時(shí)間尺度，從而獲得模擬計(jì)算所需總位移加載量；為了便于分析應(yīng)力演化過程，整個(gè)計(jì)算過程分為40個(gè)時(shí)間步進(jìn)行位移加載.

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

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

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

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

2.2 江錯(cuò)斷層面上的應(yīng)力演化

(4)

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

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

圖5 江錯(cuò)斷層面上的應(yīng)力演化結(jié)果(a) 加載時(shí)間為2500年； (b) 加載時(shí)間為7700年； (c) 加載時(shí)間為1.3萬年； (d) 加載時(shí)間為2.1萬年； (e) 加載時(shí)間為5.4萬年； (f) 對應(yīng)的加載時(shí)間為10萬年(圓圈代表2021年5月22日至2021年11月1日內(nèi)在研究區(qū)發(fā)生的震級≥3.0余震序列，其中圓圈顏色代表震源深度，圓圈大小對應(yīng)不同震級，五角星為主震)； (g) 主震、余震震源位置與江錯(cuò)斷層面的空間關(guān)系(水平投影).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ā)震斷層?xùn)|段上存在著余震活動(dòng)較少的稀疏帶(王未來等，2021；尹欣欣等，2021)，這與圖5g中江錯(cuò)斷層面上東段的低應(yīng)力區(qū)(K→L)較為相符.由此推測認(rèn)為：斷層面上局部區(qū)域應(yīng)力積累的空間差異導(dǎo)致了余震的空間分段性.結(jié)合震級大小和地震數(shù)量來看，發(fā)震斷層面上的地震序列基本可分為4段，其中有3處位于應(yīng)力集中區(qū)，分別對應(yīng)發(fā)震斷層面的東、西兩段末端及震源區(qū)域，有一處位于應(yīng)力水平較低處(圖5f).整體來看，地震序列分布與發(fā)震斷層面上的應(yīng)力積累情況具有較好的空間對應(yīng)關(guān)系，積累應(yīng)力較大位置對應(yīng)的地震數(shù)量叢集、強(qiáng)度大.此外，發(fā)震斷層面西段與震中之間的地震序列密集段對應(yīng)的是應(yīng)力水平較低處；考慮到模型的構(gòu)建情況，東、西兩段以及震源處的應(yīng)力集中區(qū)均位于斷層傾角變化的過渡處，在構(gòu)造運(yùn)動(dòng)自西向東推擠作用下，易于在傾角變化處積累應(yīng)力，從而誘發(fā)地震；而低應(yīng)力處位于兩分段的中間位置，并未出現(xiàn)明顯的應(yīng)力積累，這也說明，本研究的模型構(gòu)建在分段劃分上還可進(jìn)一步細(xì)化.

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

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

圖6 地表構(gòu)造應(yīng)力場演化結(jié)果(a) 加載時(shí)間為3.1萬年； (b) 加載時(shí)間為5.1萬年； (c) 加載時(shí)間為8.2萬年； (d) 加載時(shí)間為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.

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

圖7 瑪多地震地表破裂現(xiàn)場調(diào)查結(jié)果(a—f分別對應(yīng)圖6d中的1—6處)(a) 沙土液化； (b) 間斷式拉擠破裂； (c) 野馬灘大橋垮塌； (d) 地表拉張開裂； (e) 黃河鄉(xiāng)大橋橋墩破損； (f) 無人機(jī)探視.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.

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

圖8 震源所在深度水平面內(nèi)的構(gòu)造應(yīng)力場演化結(jié)果(虛線為地表斷層線)(a) 加載時(shí)間為5000年； (b) 加載時(shí)間為1.6萬年； (c) 加載時(shí)間為3.6萬年； (d) 加載時(shí)間為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)造應(yīng)力時(shí)空演化過程(圖6和圖8)可知：首先，瑪多地震主震震源所在水平面內(nèi)應(yīng)力集中現(xiàn)象主要體現(xiàn)在東段末端區(qū)域以及整個(gè)江錯(cuò)斷層上，而地表的應(yīng)力集中主要體現(xiàn)在東、西段末端區(qū)域(震中位置未見應(yīng)力集中現(xiàn)象)，說明斷層走向及傾角的變化對整個(gè)斷層面上構(gòu)造應(yīng)力的積累與集中發(fā)育過程以及三維空間展布具有顯著影響.

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

2.4 討論

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

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

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

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

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

3 結(jié)論

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

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

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

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

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