龍門(mén)山斷裂帶中段及鄰區(qū)構(gòu)造變形特征及應(yīng)力場(chǎng)演化的數(shù)值模擬研究*

付 平，張新輝，尹健民，韓曉玉，周春華

(長(zhǎng)江科學(xué)院 水利部巖土力學(xué)與工程重點(diǎn)實(shí)驗(yàn)室，湖北 武漢 430010)

0 引言

巖石圈上部應(yīng)力場(chǎng)是上覆巖層重力、大尺度力源(板塊邊界力、地球動(dòng)力過(guò)程產(chǎn)生的力)和小尺度力源(巖體介質(zhì)各向異性、侵入、開(kāi)挖卸荷等產(chǎn)生的非均勻力)共同作用的結(jié)果(Zobacketal，1989)。影響地應(yīng)力場(chǎng)的因素眾多，有板塊相互作用、地幔熱對(duì)流、地球內(nèi)應(yīng)力、地心引力、地球旋轉(zhuǎn)、巖漿浸入和地殼非均勻擴(kuò)容等，另外，溫度和水壓梯度、地表剝蝕及其他物理化學(xué)變化也可引起相應(yīng)的應(yīng)力場(chǎng)變化，但其中只有重力和構(gòu)造應(yīng)力能達(dá)到較高量級(jí)(劉允芳等，2014；姚瑞等，2017)，即自重和地質(zhì)構(gòu)造作用是影響地應(yīng)力場(chǎng)的主要因素。

龍門(mén)山是青藏高原巴顏喀拉塊體和中國(guó)東部華南地塊體的邊界構(gòu)造帶，地質(zhì)構(gòu)造十分復(fù)雜。目前，對(duì)于龍門(mén)山地區(qū)的變形場(chǎng)及應(yīng)力場(chǎng)已有不少研究成果：如白玉柱和周慶(2013)通過(guò)龍門(mén)山斷裂映秀—青川段的有限元模型探討了在垂直于斷層走向的應(yīng)力作用下逆沖構(gòu)造運(yùn)動(dòng)造成的地表變形和斷層附近的應(yīng)力場(chǎng)變化；柳暢等(2012)和萬(wàn)永魁等(2017)通過(guò)構(gòu)建巖石圈的三維粘彈性有限元模型研究了龍門(mén)山斷裂帶及其周緣區(qū)域的形變運(yùn)動(dòng)和應(yīng)力累積狀態(tài)；陳棋福等(2015)和尹力等(2018)分別基于粘彈性接觸和二維粘彈塑性計(jì)算模擬了強(qiáng)震復(fù)發(fā)過(guò)程中龍門(mén)山斷裂帶周?chē)貐^(qū)地殼地表的變形分布與演化。上述有限元模型大多建立在斷裂面是平面或?yàn)閱蝹€(gè)斷裂面的假定情況下，主要研究地震復(fù)發(fā)周期內(nèi)的時(shí)間演化歷程，對(duì)龍門(mén)山地區(qū)疊瓦狀多斷裂面、內(nèi)鏟型逆沖運(yùn)動(dòng)形成的變形及應(yīng)力演化的研究較少，而這些是歷史持續(xù)構(gòu)造運(yùn)動(dòng)的直接體現(xiàn)，值得深入研究。本文綜合考慮龍門(mén)山斷裂帶中段及鄰區(qū)的地形地貌、地殼厚度差異、深部巖體物理力學(xué)性質(zhì)以及板塊運(yùn)動(dòng)特征，基于最新地質(zhì)構(gòu)造剖面，建立龍門(mén)山斷裂帶中段的二維摩擦接觸有限元模型，模擬其在構(gòu)造擠壓作用下的變形以及應(yīng)力累積過(guò)程，探討板塊運(yùn)動(dòng)與構(gòu)造應(yīng)力之間的關(guān)系，研究斷裂運(yùn)動(dòng)方式和地應(yīng)力場(chǎng)演化規(guī)律。

1 研究區(qū)概況

龍門(mén)山斷裂帶南起瀘定、天全，向東北經(jīng)灌縣、茂汶、北川、廣元后進(jìn)入陜西勉縣一帶，總體呈NE—SW向展布，傾向NW，長(zhǎng)約500 km，寬30～50 km，主要為北西側(cè)的汶川—茂縣斷裂和青川斷裂、中央的映秀—北川斷裂、南東側(cè)的灌縣—安縣斷裂4條近平行斷裂及其控制的構(gòu)造巖片組成的前展式逆沖推覆構(gòu)造帶(岳沖等，2017；李海兵等，2018)，其中映秀—北川斷裂也稱(chēng)主中央斷裂，如圖1所示。該區(qū)域處于中國(guó)西部地質(zhì)、地貌和氣候的陡變帶，地形起伏劇烈，高差達(dá)(3 500±500)m，莫霍面急劇變化，地殼厚度由西北方青藏高原東緣的近70 km減至東南方四川盆地的近40 km(王椿鏞等，2003；Huangetal，2002)。相對(duì)穩(wěn)定的華南地塊阻擋了巴顏喀拉塊體SE向的逃逸，使得在巴顏喀拉塊體東邊界產(chǎn)生NW—SE向或WNW—ESE向的水平擠壓，造成應(yīng)力在龍門(mén)山推覆構(gòu)造帶上高度累積，不僅導(dǎo)致青藏高原的持續(xù)隆升，還使其周緣地區(qū)產(chǎn)生強(qiáng)烈的地殼變形和頻繁的地震活動(dòng)(朱介壽等，2017)，僅2008年以來(lái)，在龍門(mén)山斷裂帶上就發(fā)生了2次7級(jí)以上強(qiáng)震：2008年5月12日汶川MS8.0地震和2013年4月20日蘆山MS7.0地震。

圖1 龍門(mén)山地區(qū)構(gòu)造背景示意圖Fig.1 Simplified tectonic setting in the Longmenshan areas

2 龍門(mén)山斷裂帶中段計(jì)算模型

2.1 有限元模型

巴顏喀拉塊體沿先存主干斷裂向SE方向逃逸，在龍門(mén)山地區(qū)受到四川盆地阻擋，其應(yīng)力在斷裂帶上不斷地累積和調(diào)整，因此，認(rèn)為在這一地球動(dòng)力學(xué)過(guò)程中平面應(yīng)變假定近似成立。汶川地震震源機(jī)制解(鄭勇等，2009；豐成君等，2018)揭示龍門(mén)山斷裂帶除北段表現(xiàn)出明顯的走滑性質(zhì)外，其余部位均表現(xiàn)為純逆沖性質(zhì)，動(dòng)力源為NWW—SEE向的水平擠壓，所以本文將龍門(mén)山中段構(gòu)造環(huán)境簡(jiǎn)化成青藏高原向SE向推擠、受四川盆地阻擋的二維平面應(yīng)變模型。選取如圖1所示的垂直于龍門(mén)山斷裂帶走向的剖面AA′為研究對(duì)象，該剖面跨越巴顏喀拉塊體、龍門(mén)山斷裂帶和四川盆地3個(gè)主要地質(zhì)構(gòu)造單元，剖面沿程地形地貌變化見(jiàn)圖2。圖2顯示地形起伏集中出現(xiàn)在龍日壩與龍門(mén)山斷裂帶之間；映秀—北川斷裂為龍門(mén)山地形陡變帶的邊界；青藏高原東緣的龍日壩斷裂帶具有右旋走滑兼逆沖性質(zhì)，在一定程度上吸收了巴顏喀拉塊體向四川盆地的擠壓作用，龍泉山背斜是川西前陸盆地沉積邊緣，被認(rèn)為是龍門(mén)山造山帶對(duì)前陸的最大影響范圍。

為避免龍日壩斷裂帶活動(dòng)的影響并體現(xiàn)四川盆地的阻擋作用，筆者綜合上述區(qū)域地質(zhì)背景及地球物理研究給出的最新深部構(gòu)造特征(張新彥等，2017)，構(gòu)建了如圖3所示的數(shù)值計(jì)算模型，模型長(zhǎng)度設(shè)定為160 km、深度約40 km。龍門(mén)山斷裂帶出露地表處位于模型中部，汶川—茂縣斷裂(F1)、映秀—北川斷裂(F2)和灌縣—安縣斷裂(F3)以疊瓦式分布，斷裂傾角由淺至深逐漸放緩，F(xiàn)3和F2在深度約16 km處合并。F2向深部延伸至約23 km處與F1合并，采用四節(jié)點(diǎn)平面應(yīng)變單元進(jìn)行有限元網(wǎng)格劃分，模型共有16 213個(gè)節(jié)點(diǎn)、15 709個(gè)單元，在斷裂部位進(jìn)行網(wǎng)格加密。

巖石介質(zhì)的本構(gòu)關(guān)系采用摩爾-庫(kù)侖屈服準(zhǔn)則描述，物理力學(xué)參數(shù)包括彈性模量、內(nèi)聚力、內(nèi)摩擦角、剪脹角和泊松比。龍門(mén)山地區(qū)主要發(fā)育寶興雜巖和彭灌雜巖，該巖體抵抗破裂的強(qiáng)度大，能夠積累較高的彈性應(yīng)變能密度。計(jì)算模型中的彈性模量和密度參考花崗閃長(zhǎng)巖的巖石力學(xué)特性試驗(yàn)參數(shù)(姚琪等，2012；李玉江等，2013)，并根據(jù)各單元質(zhì)心位置坐標(biāo)按照均勻變化的梯度形式進(jìn)行賦值，以體現(xiàn)研究區(qū)地殼速度結(jié)構(gòu)從淺部到深部表現(xiàn)出的“低速—高速”層狀分布特點(diǎn)(柳暢等，2014)，其它參數(shù)為定值。將活動(dòng)斷裂上下兩盤(pán)的相對(duì)運(yùn)動(dòng)看作是摩擦接觸運(yùn)動(dòng)，在斷裂面兩側(cè)設(shè)置接觸對(duì)，進(jìn)行斷層的滑移模擬，接觸對(duì)主要參數(shù)包括摩擦系數(shù)、滑移系數(shù)、法向接觸剛度等。本文模型參數(shù)的設(shè)置參考了前人研究成果(朱廣安等，2016；王紅才等，2012；龍小剛，朱守彪，2015)，巖體和斷層具體參數(shù)見(jiàn)表1。

圖2 龍門(mén)山斷裂帶中段剖面AA′地形變化圖Fig.2 Topography of the profile AA′ across the middle segment of Longmenshan fault zone

圖3 龍門(mén)山斷裂帶中段及鄰區(qū)有限元模型示意圖Fig.3 Schematic of finite element model of the middle segment of Longmenshan fault zone and its adjacent areas

表1 模型巖體和斷層參數(shù)設(shè)置Tab.1 Parameters setting of rocks and fault in the lithosphere

2.2 動(dòng)力學(xué)邊界條件

將龍門(mén)山斷裂帶所處的構(gòu)造環(huán)境簡(jiǎn)化成位移約束，作為驅(qū)動(dòng)計(jì)算模型的動(dòng)力學(xué)邊界條件。考慮到垂直于斷裂帶在SE方向上表現(xiàn)為擠壓連續(xù)變形特征，選取包含剖面AA′在內(nèi)的、長(zhǎng)約600 km、寬約200 km的長(zhǎng)方形區(qū)域，如圖4a所示，將該區(qū)域內(nèi)的GPS觀測(cè)速度結(jié)果(Zhengetal，2017)轉(zhuǎn)換為以華南地塊為基準(zhǔn)，并根據(jù)龍門(mén)山斷裂帶的走向分解為垂直于和平行于斷裂的分量。由圖4b可見(jiàn)，在龍日壩和龍門(mén)山斷裂帶附近存在2處較為明顯的速度降，這是由于龍日壩和龍門(mén)山斷裂先后吸收了巴顏喀拉塊體的部分推擠作用，向SE跨越龍日壩斷裂帶后的構(gòu)造推擠速度降為約2.1 mm/a。已有研究結(jié)果(Densmoreetal，2007；Zhangetal，2004)亦表明，在華南參考框架下橫跨龍門(mén)山斷裂帶中段NW—SE方向上的水平縮短速率約1.0～4.0 mm/a。因此，計(jì)算模型的邊界條件設(shè)置如下：模型右側(cè)水平方向位移固定、垂直方向位移自由，底部水平方向位移自由、垂直方向位移固定，頂部完全自由?？紤]到中下地殼與地表運(yùn)動(dòng)速度的差異(胡幸平等，2011)，并根據(jù)Li等(2011)利用汶川地震前數(shù)字地震臺(tái)網(wǎng)記錄的波形資料估算得出龍門(mén)山斷裂帶周?chē)畈炕瑒?dòng)速率約為GPS和地質(zhì)等淺表觀測(cè)滑動(dòng)速率值的2倍這一結(jié)論，這里假定構(gòu)造推擠速度沿深度線(xiàn)性變化，如圖3所示，在模型左側(cè)施加由地表2.1 mm/a增加到底部4.2 mm/a的水平向相對(duì)速度荷載。

圖4 研究區(qū)相對(duì)于華南地塊的GPS速度場(chǎng)(a)與垂直于龍門(mén)山斷裂帶走向的速度分量(b)Fig.4 GPS velocities relative to south China block in the study area(a)and the components perpendicular to the strike of Longmenshan fault zone(b)

3 龍門(mén)山斷裂帶中段及鄰區(qū)變形特征

合理地給定初始狀態(tài)是地球動(dòng)力學(xué)數(shù)值模擬的關(guān)鍵前提條件(祝愛(ài)玉等，2016)，本文模型的初始應(yīng)力場(chǎng)為靜巖應(yīng)力場(chǎng)，模型中任意一點(diǎn)各方向的正應(yīng)力等于其上覆巖體的重力，且無(wú)剪應(yīng)力，具體如下：施加重力后，再經(jīng)自重應(yīng)力場(chǎng)平衡算法消除掉因自重產(chǎn)生的節(jié)點(diǎn)位移，在這個(gè)應(yīng)力場(chǎng)近似為自重應(yīng)力、位移場(chǎng)近似為零的初始狀態(tài)下對(duì)模型進(jìn)行速度加載。從應(yīng)力計(jì)算結(jié)果看，歷經(jīng)約67萬(wàn)年的持續(xù)運(yùn)動(dòng)加載，應(yīng)力變化不顯著，基本趨于穩(wěn)定，模型進(jìn)入到一個(gè)相對(duì)穩(wěn)定的構(gòu)造狀態(tài)。

速度加載至67萬(wàn)年，垂直位移分布如圖5所示，汶川—茂縣斷裂、映秀—北川斷裂、灌縣—安縣斷裂出露地表處上下兩盤(pán)垂直方向上的相對(duì)位移演化曲線(xiàn)如圖6所示，F(xiàn)1，F(xiàn)2，F(xiàn)3的垂直方向相對(duì)滑移量分別為72.5，271.2，186.4 m，由此得出其平均逆沖速率分別為0.11，0.40，0.28 mm/a。GPS觀測(cè)和地質(zhì)年代學(xué)研究(馬保起等，2005；李勇等，2006；張培震等，2008)顯示龍門(mén)山斷裂帶的逆沖速率如下：汶川—茂縣斷裂為0.3～0.8 mm/a、映秀—北川斷裂為0.4～1.2 mm/a、灌縣—安縣斷裂為0.3 mm/a。該結(jié)果表明，數(shù)值模擬得到的運(yùn)動(dòng)速率與實(shí)際觀測(cè)結(jié)果基本一致。

圖7顯示這3條活動(dòng)斷裂的上盤(pán)抬升速率均大于下盤(pán)，汶川—茂縣斷裂左側(cè)的巴顏喀拉塊體和龍門(mén)山斷裂帶區(qū)域呈水平縮短和不斷抬升趨勢(shì)，灌縣—安縣斷裂右側(cè)的四川盆地在垂直向則保持相對(duì)穩(wěn)定。這3條斷裂的逆沖速率關(guān)系為vF2＞vF3＞vF1，汶川震后調(diào)查表明(何宏林等，2008；張培震，2008)，沿汶川—茂縣沒(méi)有發(fā)生破裂，而在映秀—北川—平通—南壩一帶產(chǎn)生了嚴(yán)重的地表破裂和錯(cuò)動(dòng)，沿灌縣—安縣也形成了較長(zhǎng)的地表破裂，可以看到，本文模擬結(jié)果支持了上述地震地質(zhì)研究結(jié)論。

已有研究認(rèn)為龍門(mén)山地區(qū)自西向東形成系列的逆沖斷層和推覆構(gòu)造，逆沖斷層將青藏高原東部地殼向東逃逸轉(zhuǎn)換為隆升(Xuetal，2009)，中央的映秀—北川斷裂的強(qiáng)逆沖性使此斷裂長(zhǎng)期發(fā)震從而產(chǎn)生較大的垂直位移累積量，盡管與同震累積垂直位移量和地層剝蝕量相比，持續(xù)逆沖作用所產(chǎn)生的垂直位移量較小，但仍然表明映秀—北川斷裂對(duì)龍門(mén)山地區(qū)的隆升起主要貢獻(xiàn)。根據(jù)Slemmons和Depolo(1986)提出的活動(dòng)斷裂分類(lèi)方案，龍門(mén)山斷裂屬于中低活動(dòng)速率斷裂帶，而該區(qū)域地震活動(dòng)頻繁且強(qiáng)度大，這也從側(cè)面說(shuō)明：雖然滑動(dòng)速率是衡量斷裂活動(dòng)性的一個(gè)重要參數(shù)，但僅依靠地表調(diào)查和震前GPS觀測(cè)得出的活動(dòng)速率評(píng)估活動(dòng)斷裂區(qū)的地震危險(xiǎn)程度具有一定的局限性，還需結(jié)合斷裂深部動(dòng)力學(xué)過(guò)程和區(qū)域地應(yīng)力狀態(tài)開(kāi)展系統(tǒng)分析研究。

圖5 龍門(mén)山斷裂帶中段及鄰區(qū)垂直位移云圖

Fig.5 Vertical displacement distribution of the middle segment of Longmenshan fault zone and its adjacent areas

圖6 龍門(mén)山斷裂帶中段逆沖位移演化曲線(xiàn)Fig.6 Thrust displacement curve of the middle segment of the Longmenshan fault zone

圖7 地表水平向和垂直向變形速率Fig.7 Horizontal and vertical deformation rates of the surface

4 構(gòu)造運(yùn)動(dòng)過(guò)程中應(yīng)力場(chǎng)演化規(guī)律

研究人員在龍門(mén)山斷裂帶區(qū)域開(kāi)展了包括應(yīng)力解除、水壓致裂等方法在內(nèi)的大量原地應(yīng)力測(cè)試工作(安其美等，2004；陳群策等，2012；秦向輝等，2013)。受汶川地震同震應(yīng)力效應(yīng)影響，汶川地震后龍門(mén)山斷裂帶實(shí)測(cè)地應(yīng)力沿?cái)嗔炎呦蚓哂忻黠@分段特征，本文選取斷裂帶中段剖面AA′附近的地應(yīng)力實(shí)測(cè)點(diǎn)(測(cè)試鉆孔位置見(jiàn)圖1)，整理其中巖心相對(duì)完整的深部測(cè)試段獲得的主應(yīng)力數(shù)據(jù)列于表2。為探討區(qū)域構(gòu)造運(yùn)動(dòng)對(duì)地應(yīng)力場(chǎng)的影響，在有限元模型F1，F(xiàn)2，F(xiàn)3這3條斷裂上下兩盤(pán)的不同深度位置上設(shè)置若干組主應(yīng)力監(jiān)測(cè)點(diǎn)，提取各監(jiān)測(cè)點(diǎn)主應(yīng)力量值列于表3，其中，中間主應(yīng)力為面外主應(yīng)力。為檢驗(yàn)?zāi)M結(jié)果的可靠性，將監(jiān)測(cè)點(diǎn)的特征值主應(yīng)力比λ=σ1/σ3與實(shí)測(cè)點(diǎn)的特征值kHv=σH/σv進(jìn)行對(duì)比，λ和kHv均表征最大水平主應(yīng)力作用強(qiáng)度。由表2，3可以看到，模型監(jiān)測(cè)點(diǎn)λ的范圍為1.88～3.17，且基本上呈現(xiàn)深度越深比值越小的特點(diǎn)；實(shí)測(cè)點(diǎn)kHv的范圍為1.39～3.03，同樣呈現(xiàn)上述規(guī)律，但數(shù)值離散性較強(qiáng)，原因在于淺部應(yīng)力易受巖體結(jié)構(gòu)的影響，總體上看，模擬計(jì)算得出的特征值與實(shí)測(cè)結(jié)果大致相符。

表2 龍門(mén)山斷裂帶中段及鄰區(qū)地應(yīng)力實(shí)測(cè)結(jié)果Tab.2 Results of in-situ stress measurement in the middle segment of Longmenshan fault zone and its adjacent areas

表3 模型各監(jiān)測(cè)點(diǎn)主應(yīng)力量值Tab.3 Principal stresses of monitoring points in the simulation model

利用2008年汶川MS8.0主震及余震序列龍門(mén)山中段地區(qū)的震源機(jī)制解(Global CMT，2013；崔效鋒等，2011)，對(duì)該地區(qū)的構(gòu)造應(yīng)力場(chǎng)進(jìn)行反演，得到了主應(yīng)力方向和應(yīng)力形因子R的最優(yōu)解(圖8)。圖8顯示，最大主應(yīng)力軸方位角100.7°，傾角2.6°；中間主應(yīng)力軸方位角9.8°，傾角20.8°；最小主應(yīng)力軸方位角196.5°，傾角70.6°，R為0.72，應(yīng)力結(jié)構(gòu)為逆沖型，顯示該地區(qū)主要受近NWW向的水平擠壓應(yīng)力場(chǎng)控制。本文所構(gòu)建二維平面模型得到的最大主應(yīng)力方向?yàn)槠拭娴娜∠騈WW向，這里統(tǒng)計(jì)了3條斷裂監(jiān)測(cè)點(diǎn)處最大主應(yīng)力傾角的平均值，F(xiàn)1上下盤(pán)分別為3.2°和7.3°，F(xiàn)2上下盤(pán)分別為11.1°和8.8°，F(xiàn)3上下盤(pán)分別為12.3°和2.7°，均為小角度近水平，可以看到數(shù)值模擬得到的主應(yīng)力方向和傾角與震源機(jī)制解反演結(jié)果較接近。應(yīng)力形因子R=(σ1-σ2)/(σ1-σ3)反映了3個(gè)主應(yīng)力的相對(duì)大小，表3中各監(jiān)測(cè)點(diǎn)的R值大于反演得到的最佳值，原因在于本文對(duì)構(gòu)造環(huán)境的平面應(yīng)變簡(jiǎn)化使得沿?cái)嗔褞ё呦虻淖冃问艿綇?qiáng)約束，未能考慮NNE向中間主應(yīng)力產(chǎn)生的壓張效應(yīng)。

圖8 汶川地震序列龍門(mén)山中段地區(qū)震源機(jī)制解反演結(jié)果

Fig.8 Inversion results of focal mechanism solutions of Wenchuan earthquake sequence in the middle segment of Longmenshan fault zone

由圖9所示的斷裂帶附近監(jiān)測(cè)點(diǎn)主應(yīng)力隨時(shí)間演化曲線(xiàn)可以看出，隨著持續(xù)的區(qū)域構(gòu)造運(yùn)動(dòng)，主應(yīng)力基本呈不斷增大趨勢(shì)，符合擠壓逆沖型構(gòu)造環(huán)境，主應(yīng)力量值的變化過(guò)程大致可分為3個(gè)階段：較低水平緩慢增加的低態(tài)穩(wěn)定階段、持續(xù)增加的應(yīng)力累積階段和達(dá)到臨界點(diǎn)后小幅度變化的臨界穩(wěn)定階段。歷經(jīng)10～15萬(wàn)年的低態(tài)穩(wěn)定階段構(gòu)造運(yùn)動(dòng)后，各監(jiān)測(cè)點(diǎn)主應(yīng)力量值開(kāi)始顯著增加，但進(jìn)一步達(dá)到臨界狀態(tài)所需時(shí)間不同，深部要晚于淺部進(jìn)入這個(gè)臨界相對(duì)穩(wěn)定階段。F3斷層6.5 km、F1斷層6.7 km、F2斷層8.1 km處監(jiān)測(cè)點(diǎn)達(dá)到臨界穩(wěn)定狀態(tài)的時(shí)間約為45～50萬(wàn)年，F(xiàn)1斷層15.2 km處監(jiān)測(cè)點(diǎn)達(dá)到臨界穩(wěn)定狀態(tài)的時(shí)間約為60萬(wàn)年，F(xiàn)1斷層20.4 km處監(jiān)測(cè)點(diǎn)則看不到明顯的穩(wěn)定階段，60萬(wàn)年后主應(yīng)力量值仍表現(xiàn)出降低—升高—降低的波動(dòng)性，說(shuō)明當(dāng)活動(dòng)斷裂帶系統(tǒng)的構(gòu)造應(yīng)力達(dá)到一個(gè)臨界值后，應(yīng)力將處于不斷調(diào)整與重新分配的過(guò)程，應(yīng)力的累積和釋放會(huì)形成一個(gè)相對(duì)穩(wěn)定的自平衡狀態(tài)，且深度越深，這種自發(fā)式的應(yīng)力調(diào)整活動(dòng)越頻繁、越劇烈。臨界穩(wěn)定階段出現(xiàn)的應(yīng)力波動(dòng)被認(rèn)為是震前的自發(fā)式調(diào)整，是發(fā)震的必要前提條件，若該階段巖體在無(wú)法承受更高水平的應(yīng)力累積時(shí)快速破裂，應(yīng)力迅速降低，這種長(zhǎng)期累積的應(yīng)變能突然釋放即預(yù)示著地震的發(fā)生。已有研究(聞學(xué)澤，2018；孫云強(qiáng)，羅綱，2018)認(rèn)為斷裂構(gòu)造帶在經(jīng)過(guò)長(zhǎng)時(shí)間的應(yīng)力累積期后，將進(jìn)入能量釋放的地震叢集期，而在發(fā)震后應(yīng)力又會(huì)以一定速度繼續(xù)累積，2008年汶川、2013年蘆山以及2017年九寨溝地震事件也初步證明龍門(mén)山地區(qū)處于應(yīng)力累積—地震發(fā)生—應(yīng)力再累積的地震周期。由此看出，數(shù)值模擬得出的低態(tài)穩(wěn)定—應(yīng)力累積—臨界穩(wěn)定的階段式構(gòu)造應(yīng)力演化過(guò)程與該地區(qū)現(xiàn)實(shí)中的地震周期性活動(dòng)特征相呼應(yīng)。

圖9 不同深度處斷裂帶附近主應(yīng)力演化曲線(xiàn)Fig.9 Evolution of principal stress near faults at different depths

5 結(jié)論與討論

本文通過(guò)對(duì)龍門(mén)山斷裂帶中段及鄰區(qū)構(gòu)造變形以及應(yīng)力場(chǎng)演化進(jìn)行限元數(shù)值模擬，并結(jié)合已有研究成果，得出如下結(jié)論：

(1)龍門(mén)山斷裂帶屬于擠壓逆沖型斷裂，汶川—茂縣、映秀—北川、灌縣—安縣斷裂的逆沖速率分別約0.11，0.40，0.28 mm/a，與實(shí)際觀測(cè)結(jié)果基本一致，中央的映秀—北川斷裂逆沖性質(zhì)最為強(qiáng)烈，對(duì)龍門(mén)山隆升起主要貢獻(xiàn)，龍門(mén)山左側(cè)的巴顏喀拉塊體和斷裂帶區(qū)域呈不斷抬升趨勢(shì)，右側(cè)的四川盆地保持相對(duì)穩(wěn)定。

(2)數(shù)值模擬結(jié)果顯示3條主要斷裂附近6.5～20.4 km深度范圍內(nèi)巖體大小主應(yīng)力比為1.88～3.17，且表現(xiàn)出深度越深主應(yīng)力比越小的規(guī)律，總體上與實(shí)際測(cè)試結(jié)果相符，數(shù)值模擬獲得的主應(yīng)力方向和傾角與利用該地區(qū)震源機(jī)制解反演得到的構(gòu)造應(yīng)力場(chǎng)特征接近。

(3)龍門(mén)山斷裂帶中段區(qū)域構(gòu)造應(yīng)力呈現(xiàn)低態(tài)穩(wěn)定—應(yīng)力累積—臨界穩(wěn)定的階段式演化過(guò)程，且深部要晚于淺部達(dá)到臨界穩(wěn)定狀態(tài)，其臨界穩(wěn)定階段的應(yīng)力值仍表現(xiàn)出降低—升高—降低的波動(dòng)性，這種階段式應(yīng)力演化過(guò)程與該區(qū)域所處的應(yīng)力累積—地震發(fā)生—應(yīng)力再累積的地震活動(dòng)特征相呼應(yīng)。

(4)通過(guò)地表調(diào)查和震前GPS觀測(cè)揭示的斷裂活動(dòng)速率分析評(píng)價(jià)活動(dòng)斷裂區(qū)的地震危險(xiǎn)程度存在一定局限，需要進(jìn)一步結(jié)合其深部動(dòng)力學(xué)過(guò)程和區(qū)域應(yīng)力狀態(tài)進(jìn)行綜合研究。

本文側(cè)重在研究思路與方法上，通過(guò)垂直于龍門(mén)山斷裂中段走向的二維剖面模型，探討構(gòu)造擠壓作用下多斷裂面逆沖運(yùn)動(dòng)過(guò)程中的變形及應(yīng)力演化規(guī)律。數(shù)值模擬尚存在一些待改進(jìn)之處，例如地震的觸發(fā)和破裂過(guò)程未加以考慮，文中二維模型雖然是三維情形的合理簡(jiǎn)化，但該地區(qū)地質(zhì)構(gòu)造十分復(fù)雜，可考慮更多因素建立三維模型。在三維接觸或流變模型的基礎(chǔ)上嘗試引入合理的損傷及愈合機(jī)制，研究震前構(gòu)造應(yīng)力累積—震時(shí)應(yīng)力釋放與調(diào)整—震后應(yīng)力恢復(fù)再累積全過(guò)程的應(yīng)力狀態(tài)變化，對(duì)于科學(xué)認(rèn)識(shí)孕震發(fā)震機(jī)理是很有意義的。

感謝評(píng)審專(zhuān)家提出的寶貴意見(jiàn)和建議，感謝中國(guó)地震局第一監(jiān)測(cè)中心提供的GPS觀測(cè)數(shù)據(jù)。