松潘—甘孜地塊地殼和上地幔頂部現(xiàn)今變形模式數(shù)值模擬研究

萬永魁, 劉峽, 劉瑞豐, 張揚(yáng), 沈小七, 鄭智江

1 中國地震局地球物理研究所， 北京 100081 2 中國地震局第一監(jiān)測中心， 天津 300180 3 中國地震局地質(zhì)研究所， 北京 100029

0 引言

松潘—甘孜地塊隆升機(jī)制主要有側(cè)向擠出(Molnar and Tapponnier, 1975; Tapponnier et al., 2001; Hubbard and Shaw, 2009; Wang et al., 2011; Zhang, 2013)和中下地殼流(Royden et al., 1997; Clark and Royden, 2000; Clark et al., 2005a; Burchfiel et al., 2008)兩種模式.汶川MS8.0地震在龍門山斷裂帶中南段和北段分別形成6.2±0.5 m、6.5±0.5 m的最大垂直位移，對(duì)應(yīng)地殼縮短量約7.0 m，直接支持側(cè)向擠出模型(徐錫偉等，2010).然而側(cè)向擠出模式不能解釋汶川地震余震多集中在5～20 km，25～40 km處的中下地殼僅有少量地震發(fā)生(黃媛等，2008；朱艾斕等，2008；陳九輝等，2009；易桂喜等，2012).龍門山前陸盆地新生代沉積分布范圍有限，且發(fā)育不全，缺乏同期顯著的構(gòu)造撓曲，不是正常的盆—山耦合關(guān)系，這也與側(cè)向擠出逆沖推覆模式的預(yù)測相矛盾(Burchfiel et al., 1995; Chen et al., 2000; Zhang et al., 2004; Li et al., 2012; Sun et al., 2018).

龍門山斷裂兩側(cè)殼幔速度差異顯著，相對(duì)四川盆地，松潘—甘孜地塊普遍具有低速、高導(dǎo)的地球物理場特征，指示該區(qū)中下地殼和上地幔頂部可能存在殼幔流(Zhao et al., 2012; Liu et al., 2014；王緒本等，2017；王志等，2017；朱介壽等2017).松潘—甘孜地塊深20～25 km處存在厚約5 km的低阻低速層(Li et al., 2009；朱介壽，2008；劉啟元等，2009)，構(gòu)成上地殼與中下地殼“解耦”運(yùn)動(dòng)的滑脫層或剪切帶.龍門山后山、中央、前山和山前隱伏斷裂向下延伸傾角逐漸變緩，最終收斂于該剪切帶(Hubbard and Shaw, 2009; Yu et al., 2010; Wan et al., 2017；許志琴等，2007；張培震，2008).因此，目前對(duì)松潘—甘孜地塊隆升之爭，主要集中在僅受控于低阻低速層“解耦”還是同時(shí)受殼幔流共同作用這兩種主流觀點(diǎn)上.

姚琪等(2012)將低阻低速層構(gòu)建為軟弱薄層，與上地殼底面構(gòu)成接觸單元，利用與速度相關(guān)的非線性接觸摩擦準(zhǔn)則，模擬了低阻低速層對(duì)高原東緣隆升的影響，結(jié)果顯示，低阻低速層可以控制松潘—甘孜地塊快速隆升，但模型未考慮中下地殼和上地幔頂部變形的影響.尹力和羅綱(2018)、龐亞瑾等(2019)采用連續(xù)模型對(duì)青藏高原東緣垂向變形控制因素展開討論，結(jié)果顯示地形起伏、地殼結(jié)構(gòu)和密度差異、巖石圈流變性質(zhì)均是該區(qū)垂向變形的重要控制因素.汪昌亮(2012)沙箱模擬實(shí)驗(yàn)結(jié)果表明，脆性上地殼與中下地殼流先后對(duì)松潘—甘孜地塊隆升起到作用.Xie等(2020)指出巖石圈均衡和巖石圈撓曲的靜態(tài)支撐，以及下地殼流和地幔對(duì)流的動(dòng)力作用是控制松潘—甘孜地塊高海拔地形的主要機(jī)制.滕吉文等(2008，2014)指出松潘—甘孜地塊中下地殼和上地幔頂部物質(zhì)分別以埋深約20 km處的低阻低速層為第一滑移面，上地幔軟流層頂面(深100±10 km)為第二滑移面，整體向SE方向運(yùn)動(dòng)，受到四川盆地深部“剛性”物質(zhì)阻擋后，向上運(yùn)移，物質(zhì)匯聚、應(yīng)力集中，從而引發(fā)汶川地震.胡幸平等(2012)采用彈性有限元模型對(duì)高原東緣深部構(gòu)造變形模式展開討論，結(jié)果表明，松潘—甘孜地塊一側(cè)，即模型北、西邊界深100 km處水平運(yùn)動(dòng)是地表5倍的前提下，所得理論震源機(jī)制解與實(shí)際震源機(jī)制解最為吻合.此外，已有研究表明地幔對(duì)流拖曳力作為中國陸地巖石圈構(gòu)造運(yùn)動(dòng)的重要驅(qū)動(dòng)力，對(duì)青藏高原內(nèi)部地殼運(yùn)動(dòng)方向有明顯控制作用(黃建平等，2008；祝愛玉等，2019).綜上所述，松潘—甘孜地塊隆升可能與低阻低速層“解耦”，地形起伏、殼幔密度和結(jié)構(gòu)差異、巖石圈流變等因素引起地殼至上地幔頂部水平運(yùn)動(dòng)速率逐漸增加以及地幔流對(duì)巖石圈底部的拖曳作用均有關(guān)聯(lián).

汶川地震后，針對(duì)龍門山斷裂帶及其周緣區(qū)域開展了一系列地球物理探測，新成像技術(shù)的廣泛應(yīng)用，使得獲取更為精細(xì)的內(nèi)部結(jié)構(gòu)圖像成為可能(Zhang et al., 2010；朱介壽，2008；朱介壽等，2017；王緒本等，2017；王志等，2017)，為構(gòu)建二維有限元精細(xì)模型提供了條件.多期GPS觀測、跨龍門山斷裂水準(zhǔn)觀測以及最新的低溫?zé)崮陮W(xué)剝蝕速率相關(guān)研究成果(Kirby et al., 2002; Godard et al., 2009; Li et al., 2012; Tian et al., 2013, 2015; Tan et al., 2017; Wang and Shen, 2020)，為模型邊界加載提供了有效約束，同時(shí)也為模擬結(jié)果合理性檢驗(yàn)提供了有力支持.為獲得松潘—甘孜地塊地殼和上地幔頂部現(xiàn)今變形模式，本文依據(jù)跨龍門山斷裂帶探測剖面，構(gòu)建了長300 km、寬104 km的二維有限元接觸模型(龍門山斷裂帶設(shè)置為接觸單元)，在考慮巖石蠕變特性的前提下，以1991—2016年長期GPS觀測結(jié)果為初始約束(圖1a黑色箭頭，Wang and Shen, 2020)，通過改變低阻低速層蠕變參數(shù)(即是否考慮低阻低速層的存在)以及模型西邊界和底邊界加載條件，共計(jì)開展了11項(xiàng)數(shù)值模擬實(shí)驗(yàn).將模擬結(jié)果與綜合多學(xué)科、不同時(shí)間尺度獲得的隆升速率進(jìn)行對(duì)比，各模型結(jié)果橫向?qū)Ρ龋懻摿怂膳恕首蔚貕K地殼至上地幔頂部變形及物質(zhì)運(yùn)移模式，有助于進(jìn)一步理解松潘—甘孜地塊隆升以及汶川MS8.0地震孕育和發(fā)震機(jī)制.

1 模型

1.1 模型構(gòu)建

青藏高原東緣由西向東依次為高海拔強(qiáng)烈變形變質(zhì)作用的松潘—甘孜褶皺帶、地勢起伏劇烈的龍門山?jīng)_斷帶和低海拔弱變形的四川盆地(劉樹根等，2019).晚三疊世以來，受印度板塊與歐亞大陸強(qiáng)烈碰撞及中國陸地主體拼合的影響，松潘—甘孜地塊和龍門山造山帶發(fā)生強(qiáng)烈變形和沖斷隆升(劉樹根等，2009)，在東西寬約30～50 km范圍內(nèi)，西側(cè)松潘—甘孜地塊與東側(cè)四川盆地形成高達(dá)約4 km地形差異，成為整個(gè)青藏高原地形梯度最大地區(qū)(Kirby et al., 2002; Li et al., 2012；李勇等，2005).龍門山斷裂帶呈北東—南西方向，全長約500 km，寬30～50 km，自北西向南東方向發(fā)育4條主干斷裂，即汶川—茂縣斷裂(后山斷裂)、映秀—北川斷裂(中央斷裂)、灌縣—江油斷裂(前山斷裂)和大邑—廣元斷裂(山前隱伏斷裂).斷裂呈疊瓦狀向四川盆地逆沖推覆，淺部斷層傾角為60°～70°，向下延伸，傾角逐漸減小并收斂合并于埋深約20 km處的低阻低速層(Burchfiel et al., 2008；Yu et al., 2010；許志琴等，2007).低阻低速層厚約5 km，埋深于松潘—甘孜地塊20～25 km處(Yu et al., 2010；劉啟元等，2009).龍門山造山帶南段、中段和北段分別以寶興雜巖、彭灌雜巖和轎子頂雜巖為核心(劉樹根等，2009)，具有強(qiáng)度大、能夠積累高密度彈性應(yīng)變能的特性(張培震等，2008).阿壩—雙流人工地震探測剖面顯示龍門山斷裂帶兩側(cè)地殼厚度變化顯著，松潘—甘孜地塊至四川盆地寬約80 km范圍，地殼厚度由約60 km迅速遞減至約45 km(朱介壽，2008).

依據(jù)阿壩—雙流人工地震探測剖面(圖1黑色粗實(shí)線)，參考前人數(shù)值模擬相關(guān)模型構(gòu)架(Liu et al., 2015；朱守彪和張培震，2009；張竹琪等，2010；柳暢等，2014；陳棋福等，2015；祝愛玉等，2016；萬永魁等，2017)，本文構(gòu)建了長300 km、寬104 km的二維有限元接觸模型.模型中龍門山4條主干斷裂簡化為1條，按照上陡下緩“鏟式”結(jié)構(gòu)，在20 km深度處與低阻低速層頂界相切.松潘—甘孜地塊與四川盆地存在顯著地形高差，中下地殼埋深差異更為顯著，故模型上地殼設(shè)置了寬30 km過渡帶，中下地殼設(shè)置了寬80 km過渡帶.上地殼在過渡帶內(nèi)以強(qiáng)度較高的彭灌雜巖、寶興雜巖及轎子頂雜巖為核心，其彈性模量為四川盆地的1.2倍.松潘—甘孜地塊0～4 km范圍代表高海拔地形，埋深20～25 km范圍為低阻低速層(圖2a黃色區(qū)域)，用于模擬上地殼與中下地殼的“解耦”.模型整體構(gòu)架見圖2a.

圖1 研究區(qū)水準(zhǔn)觀測路線、GPS測站速度矢量(a)及垂直于龍門山斷裂速度剖面(b)Fig.1 Leveling observation route and velocity vector of GPS in the study area(a) and velocity profile perpendicular to Longmenshan Fault (b)

1.2 參數(shù)設(shè)置

巖石在長時(shí)間載荷作用下應(yīng)力、應(yīng)變符合冪指數(shù)關(guān)系(Kirby,1983)，滿足

(1)

B=Aexp(-E/RT),(2)

對(duì)于處于柔性變形階段的巖石，其等效黏滯系數(shù)可利用公式(3)計(jì)算：

(3)

中國陸地地殼和上地幔頂部等效黏滯系數(shù)一般在1019～1024Pa·s，相對(duì)穩(wěn)定，因此可以根據(jù)等效黏滯系數(shù)計(jì)算蠕變系數(shù)，即

(4)

本文在應(yīng)變率設(shè)定為10-14s-1前提下(石耀霖和曹建玲，2008)，計(jì)算得到松潘—甘孜地塊和四川盆地巖石圈各層蠕變系數(shù).模型各層相關(guān)參數(shù)見表1和表2.

表1 松潘—甘孜地塊巖石圈物質(zhì)參數(shù)Table 1 Material parameters of rocks in the lithosphere in Songpan-Garzê block

表2 四川盆地巖石圈物質(zhì)參數(shù)Table 2 Material parameters of rocks in the lithosphere in Sichuan basin

1.3 網(wǎng)格劃分

模型由面單元plane182組成，龍門山斷裂由二維接觸單元contact171、targe169組成(圖2a紅色部位).采用三角形單元對(duì)模型進(jìn)行網(wǎng)格劃分，網(wǎng)格尺寸約1 km，網(wǎng)格劃分結(jié)果見圖2b，單元總數(shù)31749，節(jié)點(diǎn)總數(shù)32116個(gè).

圖2 模型構(gòu)架及重力加載過程邊界約束(a)和網(wǎng)格劃分結(jié)果(b)Fig.2 Model structure and boundaries constraint of gravity loading process (a) and meshing results (b)

2 加載計(jì)算

2.1 加載

本文通過改變低阻低速層蠕變參數(shù)(是否考慮低阻低速層的存在)、模型西邊界和底邊界加載條件，共計(jì)開展了11項(xiàng)數(shù)值模擬實(shí)驗(yàn)，詳情見表3.

表3 11項(xiàng)模擬實(shí)驗(yàn)低阻低速層和邊界設(shè)置Table 3 Soft layer and boundary settings of 11 simulated experiments

加載分兩步進(jìn)行：(1)重力加載；(2)位移(構(gòu)造)加載.重力加載時(shí)，模型東、西邊界(x=300及x=0 km)水平向固定、垂向自由，底邊界(y=-100 km)垂向固定、水平向自由，施加重力后維持1 Ma自由變形.位移加載是在重力加載的基礎(chǔ)上完成的，一方面更接近實(shí)際，另一方面在構(gòu)造變形過程中重力的影響不可忽略.依據(jù)1991—2016年GPS觀測結(jié)果，模型300 km范圍內(nèi)，垂直于龍門山斷裂帶水平壓縮速率約3 mm·a-1(圖1b黃色透明區(qū)域)，故位移加載，首先考慮在模型東邊界和底邊界約束不變的前提下，西邊界施加3 mm·a-1水平位移，作為基礎(chǔ)模型(Model 1).基于基礎(chǔ)模型，通過改變低阻低速層蠕變參數(shù)(參考已有研究由低阻低速層黏滯系數(shù)計(jì)算獲得蠕變參數(shù))，體現(xiàn)低阻低速層“解耦”對(duì)高原隆升的影響(Model 2).參考胡幸平等(2012)由地表至深100 km處，加載速率線性增加5倍，所得龍門山地區(qū)理論震源機(jī)制解與實(shí)際震源機(jī)制解最為吻合，劉峽等(2010)在華北地區(qū)現(xiàn)今地殼運(yùn)動(dòng)模擬研究中設(shè)定由地表至深100 km處，加載速率線性增加1.2倍，模擬結(jié)果與GPS觀測結(jié)果最為一致，本文通過模型西邊界加載速率隨深度線性增加(1.2、1.5、1.8、2.0、2.5倍或3.0倍)，體現(xiàn)巖石圈水平向運(yùn)動(dòng)隨深度增加對(duì)高原隆升的影響(Model 3—Model 8).最后，基于Model 5的模擬結(jié)果，通過改變底邊界位移加載條件(類比地幔拖曳力)，體現(xiàn)地幔對(duì)流拖曳力對(duì)高原隆升的影響(Model 9—Model 10).此外，不考慮Model 5中低阻低速層的作用，模擬高原隆升速率(Model 11)，通過分析兩模型隆升速率差異，進(jìn)一步研究低阻低速層“解耦”對(duì)高原隆升的作用.

2.2 計(jì)算

蠕變和接觸涉及材料非線性、幾何非線性并與變形過程密切相關(guān)，接觸狀態(tài)事前通常未知，因此有限元處理蠕變和接觸問題通常采用增量法、自動(dòng)時(shí)間步長和N-R(Newton-Raphson)迭代聯(lián)合求解.增量法首先將載荷分為Q0，Q1，Q2，…，若干步，對(duì)應(yīng)位移分別為a0，a1，a2，…，假設(shè)已知第m步載荷Qm和相應(yīng)的位移am，當(dāng)載荷增加至Qm+1(Qm+1=Qm+ΔQm)，求解位移am+1(am+1=am+Δam)，理論上如果載荷增量ΔQm足夠小，則可以計(jì)算出相應(yīng)的位移增量.選擇自動(dòng)時(shí)間步長求解，ANSYS會(huì)依據(jù)載荷增量和時(shí)間步長自動(dòng)將每一載荷步劃分為若干子步進(jìn)行求解，如果計(jì)算結(jié)果仍難以收斂，會(huì)通過增加子步數(shù)、減小子步時(shí)間，使子步載荷增量ΔQ足夠小，從而實(shí)現(xiàn)收斂.N-R迭代的目的是為了改進(jìn)計(jì)算精度，對(duì)于m+1次增量步的第n+1次迭代可表示為

3 模擬結(jié)果

3.1 重力和摩擦系數(shù)對(duì)模擬結(jié)果的影響

施加重力后模型會(huì)產(chǎn)生一定程度塌陷，將Model 5接觸單元(龍門山斷裂)摩擦系數(shù)設(shè)定為0.6，分別提取僅重力作用下0.5 Ma、1 Ma、2 Ma、4 Ma和6 Ma時(shí)模型各節(jié)點(diǎn)累計(jì)位移，計(jì)算0.5～1 Ma、>1～2 Ma、>2～4 Ma和>4～6 Ma四個(gè)時(shí)段因施加重力節(jié)點(diǎn)運(yùn)動(dòng)速率(圖3)，結(jié)果顯示，0.5～1 Ma和>1～2 Ma模型中、上地殼有一定程度的下沉，量值分別<0.06 mm·a-1和<0.03 mm·a-1，下地殼和上地幔存在由松潘—甘孜地塊向四川盆地流動(dòng)的趨勢，量值分別<0.03 mm·a-1和<0.02 mm·a-1.>2～4 Ma和>4～6 Ma中、上地殼下沉和下地殼和上地幔東向流動(dòng)趨勢進(jìn)一步減小，中、上地殼下沉速率分別<0.015 mm·a-1和<0.01 mm·a-1，反映出重力加載后隨時(shí)間增加模型逐漸趨于穩(wěn)定，這一結(jié)果與前人關(guān)于重力勢作用可能是控制青藏高原邊緣動(dòng)力學(xué)變形的主要因素相矛盾，主要是因?yàn)槟Ｐ瓦吔绮捎昧藦?qiáng)制約束，這與實(shí)際邊界條件存有差異.提取松潘—甘孜地塊和四川盆地地表平均累計(jì)塌陷量，結(jié)果顯示，0.5 Ma松潘—甘孜地塊平均累計(jì)塌陷量為1197.6 m，四川盆地為699.7 m，1.0 Ma松潘—甘孜地塊為1286.1 m，四川盆地為741.7 m(圖4a)，由此得到0.5～1.0 Ma內(nèi)松潘—甘孜地塊和四川盆地因重力加載導(dǎo)致的垂向變形速率分別為0.057 mm·a-1和0.020 mm·a-1(圖4b)，遠(yuǎn)小于兩區(qū)域隆升速率，即加載重力后經(jīng)過1 Ma的變形，因重力作用導(dǎo)致的垂向變形基本可以忽略，所以本文在重力加載1 Ma模型相對(duì)穩(wěn)定后進(jìn)行位移(構(gòu)造)加載，但此時(shí)松潘—甘孜地塊與四川盆地地形高差減小至約3500 m，與4000 m實(shí)際地形高差相差約500 m，因此本文模擬結(jié)果低估了重力勢能的影響.針對(duì)龍門山斷裂摩擦系數(shù)的選取，同樣利用Model 5，嘗試將接觸單元摩擦系數(shù)分別設(shè)置為0.2、0.4、0.6、0.8和1.0，模擬結(jié)果顯示，位移加載后節(jié)點(diǎn)174(x≈150 km)和節(jié)點(diǎn)21913(x≈250 km)在不同摩擦系數(shù)下隆升趨勢幾乎一致(圖4c、4d).地表各節(jié)點(diǎn)垂直于龍門山斷裂水平向壓縮速率(圖4e)和垂向隆升速率(圖4f)也較為一致，反映出摩擦系數(shù)對(duì)模擬結(jié)果的影響并不顯著，據(jù)此本文11項(xiàng)模擬實(shí)驗(yàn)?zāi)Σ料禂?shù)統(tǒng)一設(shè)定為0.6.

圖3 重力加載后不同時(shí)段節(jié)點(diǎn)平均運(yùn)動(dòng)速率(a) 0.5～1 Ma； (b) >1～2 Ma； (c) >2～4 Ma；(d) >4～6 Ma.Fig.3 Average velocity of nodes in different periods after gravity loading

圖4 重力加載后不同時(shí)間地表相對(duì)高程(a)和相應(yīng)時(shí)間段地表垂向速率(b)；構(gòu)造加載后Model 5在不同摩擦系數(shù)下節(jié)點(diǎn)174(c)、節(jié)點(diǎn)21913(d)相對(duì)高程變化；構(gòu)造加載20萬年Model 5在不同摩擦系數(shù)下地表節(jié)點(diǎn)水平壓縮速率(e)和垂向隆升速率(f)Fig.4 Relative elevation of surface nodes at different times after gravity loading (a) and vertical velocity in corresponding time period; relative elevation change of node 174 (c) and node 21913 (d) in Model 5 after displacement loading under different friction coefficients; surface nodes in Model 5 horizontal compression rate (e) and vertical rate (f) under different friction coefficients at loading 0.2 Ma.

圖5 Model 5地表節(jié)點(diǎn)隆升位移隨位移加載時(shí)間變化曲線Fig.5 The graph of uplift displacement of surface nodes with tectonic loading time of Model 5

圖6 Model1-Model11構(gòu)造加載20萬年時(shí)地表水平向速率及GPS擬合結(jié)果(a)、(c)和隆升速率及水準(zhǔn)結(jié)果(b)、(d)Fig.6 Simulation results of 11 models at tectonic loading of 0.2 Ma, surface horizontal rates and GPS (a)、(c) and uplift rates and leveling data (b)、(d)

3.2 模擬結(jié)果穩(wěn)定性分析

為觀察模型表面隆升位移隨時(shí)間變化，在模型表面自西向東每隔約20 km取一個(gè)節(jié)點(diǎn)(node01-node 15).圖5為模型Model 5位移加載1 Ma內(nèi)各節(jié)點(diǎn)隆升位移隨時(shí)間變化曲線(圖5a黑色虛線框放大部分見圖5c，圖5b黑色虛線框放大部分見圖5d).結(jié)果顯示，位移加載開始至約8萬年，各節(jié)點(diǎn)隆升速率(曲線斜率)逐漸增大，加載20萬年后，隆升速率接近線性變化(圖5c、5d)，反映出模擬結(jié)果基本穩(wěn)定.因此，本文選取11項(xiàng)模型位移加載20萬年時(shí)的模擬結(jié)果展開分析.

圖6為11項(xiàng)模擬計(jì)算給出的地表各節(jié)點(diǎn)水平向和垂向運(yùn)動(dòng)速率，垂直于龍門山斷裂GPS速度剖面的擬合曲線(圖1b紅色虛線，扣除整體運(yùn)動(dòng))，以及基于1960—2010年水準(zhǔn)觀測數(shù)據(jù)采用動(dòng)態(tài)平差方法獲得的觀測點(diǎn)垂向速率也顯示在圖6中(中國陸地現(xiàn)代垂直形變圖集的編制與資料整編項(xiàng)目).模擬得到的水平速率自西而東逐步減小，在松潘—甘孜地塊內(nèi)減小緩慢，跨過龍門山斷裂后，四川盆地內(nèi)減小速度加快.從曲線形態(tài)看，Model 1-Model 8由折線逐漸向曲線過渡(圖6a)，Model 9-Model 11與Model 5基本一致(圖6c).GPS擬合曲線在松潘—甘孜地塊內(nèi)與Model 6最為一致，四川盆地內(nèi)與Model 5最為一致.模擬給出的松潘—甘孜地塊垂向隆升速率明顯大于四川盆地，Model 2-Model 8垂向隆升速率逐漸增大，說明模型西邊界加載位移越大，越能夠促使地表隆升加速.水準(zhǔn)觀測結(jié)果與Model 5結(jié)果最為一致，而Model 9-Model 11結(jié)果與Model 5略有差異，反映出底邊界加載和低阻低速層對(duì)松潘—甘孜地塊隆升均有影響.考慮到GPS計(jì)算結(jié)果存在0.1～1.2 mm·a-1的誤差，且大部分測站誤差集中在0.2～0.6 mm·a-1，而Model 5水平向模擬結(jié)果與GPS擬合曲線最大誤差小于0.2 mm，我們認(rèn)為Model 5模擬結(jié)果與實(shí)際的水平向和垂向運(yùn)動(dòng)最為接近.

4 討論

4.1 松潘—甘孜地塊長期穩(wěn)定隆升速率

Hao等(2014)、杜方等(2009)依據(jù)1975—1997年跨龍門山斷裂帶水準(zhǔn)觀測數(shù)據(jù)(圖1白色圓點(diǎn))，計(jì)算得到汶川地震前松潘—甘孜地塊相對(duì)于四川盆地隆升速率高達(dá)4～6 mm·a-1.青藏高原東緣經(jīng)歷了由中生代—早新生代平緩至晚新生代突然加速的剝蝕過程(Arne et al., 1997).基于磷灰石裂變徑跡(AFT)、鋯石裂變徑跡(ZFT)和(U-Th)/He等熱年代學(xué)技術(shù)，前人研究了青藏高原東緣快速剝蝕的時(shí)間和速率(表4)，最新研究成果顯示快速剝蝕始于約10 Ma，最大剝蝕速率約1.0 mm·a-1.如果隆升速率大于剝蝕速率，則地表隆升為正值，反之為負(fù)值.假定青藏高原東緣快速隆升始于距今10 Ma，在不考慮地震破裂引起地表垂向位移的前提下，取最大剝蝕速率1.0 mm·a-1，松潘—甘孜地塊與四川盆地現(xiàn)今地形高差約為4000 m，計(jì)算得到最大隆升速率為1.4 mm·a-1(地表實(shí)際隆升速率=最大隆升速率-最大剝蝕速率，為0.4 mm·a-1)，與汶川地震前跨龍門山斷裂帶水準(zhǔn)觀測結(jié)果差異巨大.

基于中國陸地現(xiàn)代垂直形變圖集的編制與資料整編項(xiàng)目，采用動(dòng)態(tài)平差方法，中國地震局第一監(jiān)測中心完成了中國陸地多期垂直形變圖，其中1960—2010年跨龍門山斷裂帶長時(shí)段垂直形變結(jié)果顯示松潘—甘孜地塊相對(duì)四川盆地最大隆升速率約為1.5 mm·a-1(圖6b、6d黑色虛線)，與低溫?zé)崮甏鷮W(xué)結(jié)合現(xiàn)今地貌計(jì)算得到的約1.4 mm·a-1最大隆升速率較為接近.因此，Hao等(2014)、杜方等(2009)給出的汶川地震前松潘—甘孜地塊相對(duì)四川盆地4～6 mm·a-1的快速隆升，作為長期動(dòng)力地質(zhì)學(xué)數(shù)值模擬的檢驗(yàn)標(biāo)準(zhǔn)有待商榷.綜合多學(xué)科、不同時(shí)間尺度研究成果，松潘—甘孜地塊長期穩(wěn)定最大隆升速率為1.4～1.5 mm·a-1可能更為客觀.

4.2 低阻低速層對(duì)高原隆升的影響

Model 1與Model 2唯一差別是前者沒有考慮低阻低速層，本研究通過對(duì)比這兩項(xiàng)模擬結(jié)果，討論低阻低速層的作用.從隆升速率看，Model 1和Model 2在松潘—甘孜東緣均出現(xiàn)了局部“鼓包”，最大隆升速率約為1.0 mm·a-1，與1.4～1.5 mm·a-1長期穩(wěn)定最大隆升速率存有較大差異.而橫坐標(biāo)0～75 km模型范圍，Model 1的地表垂向隆升速率比Model 2小約0.2 mm·a-1(圖6b).

表4 青藏高原東緣快速剝蝕的時(shí)間和速率Table 4 Time and rate of rapid denudation in the eastern Margin of the Tibet Plateau

根據(jù)圖7a—7d顯示的水平向和垂向速率隨深度變化，兩模型運(yùn)動(dòng)速率的整體趨勢基本一致.如松潘—甘孜內(nèi)部隆升速率為0.4～0.7 mm·a-1，松潘—甘孜東緣至龍門山斷裂帶隆升速率最強(qiáng)為0.7～1.0 mm·a-1，四川盆地相對(duì)微弱為0.1～0.4 mm·a-1.所不同的是Model 2垂向隆升的范圍大于Model 1，如隆升大于0.6 mm·a-1的區(qū)域，Model 1僅分布在地表橫坐標(biāo)50～200 km范圍內(nèi)，Model 2則覆蓋了整個(gè)松潘—甘孜地塊；又如隆升速率大于0.8 mm·a-1的區(qū)域，Model 1分布于橫坐標(biāo)100～175 km范圍內(nèi)，而Model 2擴(kuò)展至模型100 km以西區(qū)域.上述結(jié)果揭示低阻低速層容易變形，有利于松潘—甘孜內(nèi)部整體隆升，但不足以形成松潘—甘孜地塊如此規(guī)模的強(qiáng)烈的地表抬升.

4.3 松潘—甘孜地塊地殼和上地幔頂部變形模式

Model 2與Model 5唯一差別是后者西邊界加載位移隨深度線性增加1.8倍，即由地表3 mm·a-1至深100 km處增至5.4 mm·a-1.根據(jù)圖7e—7f顯示的Model 5水平向和垂向速率隨深度變化，松潘—甘孜內(nèi)部，水平向運(yùn)動(dòng)在低阻低速層發(fā)生了“解耦”，上地殼運(yùn)動(dòng)速率小于中下地殼.垂向速率與Model 2相比存有較大差異：(1)Model 2地表最大隆升速率約為1.0 mm·a-1，而Model 5約為1.44 mm·a-1，隆升速率顯著增加；(2)Model 2僅在上地殼150 km附近隆升速率接近1.0 mm·a-1，Model 5隆升速率大于1.0 mm·a-1的區(qū)域基本擴(kuò)展至整個(gè)松潘—甘孜地塊，強(qiáng)烈隆升區(qū)域范圍顯著擴(kuò)展；(3)Model 2隆升中心位于模型150 km附近，Model 5向西轉(zhuǎn)移至模型100 km附近，與長時(shí)段水準(zhǔn)觀測結(jié)果更為一致(圖6b).

圖7 Model 1、Model 2、Model 5、Model 9-Model 11位移加載20萬年水平向和垂向速率隨深度變化(a) Model 1水平向速率； (b) Model 1垂向速率； (c) Model 2水平向速率； (d) Model 2垂向速率； (e) Model 5水平向速率； (f) Model 5垂向速率； (g) Model 9水平向速率； (h) Model 9垂向速率； (i) Model 10水平向速率； (j) Model 10垂向速率； (k) Model 11水平向速率；(l) Model 11垂向速率.Fig.7 Model 1, Model 2, Model 5 and Model 9-Model 11 simulation results of horizontal and uplift rate vary with depth at tectonic loading of 0.2 Ma(a) Model 1 horizontal rate; (b) Model 1 uplift rate; (c) Model 2 horizontal rate; (d) Model 2 uplift rate; (e) Model 5 horizontal rate; (f) Model 5 uplift rate; (g) Model 9 horizontal rate; (h) Model 9 uplift rate; (i) Model 10 horizontal rate; (j) Model 10 uplift rate; (k) Model 11 horizontal rate; (l) Model 11 uplift rate.

Model 1、Model 2和Model 5速度矢量見圖8.Model 1和Model 2速率減小區(qū)域主要集中在松潘—甘孜東緣及龍門山斷裂帶中下地殼和上地幔頂部(圖8a、8b)，故松潘—甘孜東緣地表隆升速率出現(xiàn)局部“鼓包”.與Model 1、Model 2相比，Model 5最顯著的特征是松潘—甘孜內(nèi)速度矢量旋轉(zhuǎn)幅度顯著增強(qiáng)，中下地殼和上地幔水平向運(yùn)動(dòng)速率快速減小區(qū)域集中在模型50～150 km處，松潘—甘孜東緣和龍門山斷裂帶(模型150～200 km)，速率減小已不顯著，這與水準(zhǔn)觀測隆升中心并非在龍門山斷裂帶，而是在其西側(cè)約100 km處的川西高原，即模型50～150 km范圍相吻合.中下地殼和上地幔物質(zhì)水平向運(yùn)動(dòng)速率快速減小與四川盆地的阻擋固然有關(guān)，但與地形起伏、地殼結(jié)構(gòu)和密度差異、巖石圈流變及地幔對(duì)流拖曳力可能均有關(guān)聯(lián)，至于每一項(xiàng)因素的影響程度還需開展進(jìn)一步研究.

圖8 Model 1、Model 2和Model 5構(gòu)造加載20萬年速度矢量模擬結(jié)果(a) Model 1速度矢量； (b) Model 2速度矢量； (c) Model 5速度矢量Fig.8 Model 1, Model 2 and Model 5 simulation results of velocity vector at tectonic loading of 0.2 Ma(a) Model 1 velocity vector； (b) Model 2 velocity vector； (c) Model 5 velocity vector.

與Model 5相比，Model 9和Model 10對(duì)模型底邊界進(jìn)行了位移加載，用于類比地幔拖曳力，Model 9底邊界加載速率采用線性遞減方式，由西邊界5.4 mm·a-1至東邊界遞減為0，Model 10采用非線性遞減方式，速率減小主要集中在松潘—甘孜東緣和龍門山斷裂帶(模型100～200 km范圍)，兩模型底邊界加載速率變化曲線見圖9.模擬結(jié)果顯示，Model 9松潘—甘孜內(nèi)地表隆升速率約1.0 mm·a-1，四川盆地內(nèi)約0.50 mm·a-1，Model 10松潘—甘孜內(nèi)地表隆升速率由0逐漸增至約1.51 mm·a-1，而后快速減小，四川盆地內(nèi)約0.35 mm·a-1，但最大隆升速率位于松潘—甘孜東緣(模型約150 km處).整體而言，Model 9和Model 10地表隆升速率與水準(zhǔn)觀測結(jié)果均存有了一定差異(圖6d)，垂向速率隨深度變化結(jié)果與Model 5模擬結(jié)果差異也較為顯著(圖7f、7h、7j).上述結(jié)果反映出隆升速率和隆升中心的位置分布對(duì)底邊界加載條件較為敏感.

圖9 Model 9、Model 10底邊界加載速率變化曲線Fig.9 Bottom boundary loading rate curve of Model 9 and Model 10

與Model 5相比，Model 11不考慮低阻低速層的作用，模擬的水平向和垂向運(yùn)動(dòng)速率見圖7k、7l，與Model 5相比，Model 11中上地殼與中下地殼水平方向未發(fā)生“解耦”，隆升中心向龍門山斷裂帶一側(cè)偏移，隆升速率大于1.2 mm·a-1的區(qū)域向深部延伸.模型西側(cè)(約0～75 km)垂向隆升速率顯著減小，地表隆升速率與水準(zhǔn)觀測結(jié)果產(chǎn)生了一定差異(圖6d)，這一差異與Model 1和Model 2在模型西側(cè)地表隆升速率存有的差異是類似的(圖6b)，再次揭示低阻低速層易于變形，可以作為中下地殼和上地幔物質(zhì)運(yùn)移的一個(gè)滑移面，構(gòu)成上地殼與中下地殼的“解耦”帶，有利于松潘—甘孜地塊整體隆升.

4.4 青藏高原東緣地殼和上地幔頂部物質(zhì)運(yùn)移模式

Model 5在構(gòu)造加載50萬年和100萬年時(shí)水平向和垂向速率隨深度變化見圖10.與構(gòu)造加載20萬年結(jié)果相比，構(gòu)造加載50萬年和100萬年時(shí)松潘—甘孜地塊整體運(yùn)動(dòng)趨勢基本一致，即水平方向上地殼與中下地殼“解耦”，垂向方向整體隆升，隆升中心位于模型約100 km處.但值得注意的是最大隆升速率卻不斷增加，構(gòu)造加載20萬年最大隆升速率約1.44 mm·a-1，50萬年增至約1.8 mm·a-1，100萬年增至約2.0 mm·a-1.Clark and Royden(2000)指出青藏高原東流物質(zhì)受堅(jiān)硬四川盆地阻擋后向兩方向流出，一部分轉(zhuǎn)至北東向，另一部分轉(zhuǎn)至南東向.鄒鎮(zhèn)宇等(2015)計(jì)算的多期GPS結(jié)果顯示，相對(duì)于穩(wěn)定華南塊體GPS速度場在青藏高原東緣逐漸衰減，并發(fā)生明顯轉(zhuǎn)向，在龍門山斷裂帶中北段附近向北東方向偏轉(zhuǎn)，在鮮水河斷裂帶附近向南東方向偏轉(zhuǎn).快剪切波偏振優(yōu)勢方向在龍門山斷裂北段為NE向，南段為NW向(石玉濤等，2009).汶川MS8.0強(qiáng)余震震源機(jī)制解P軸分布也表現(xiàn)出類似的特征(胡幸平等，2008；鄭勇等，2009).上述研究說明青藏高原東緣物質(zhì)匯聚到一定程度后，可以向兩側(cè)運(yùn)移，這與本文數(shù)值模擬獲取的隆升速率隨時(shí)間逐漸增加的結(jié)果有所不同，考慮到本文2D模型具有一定封閉性，匯聚物質(zhì)并不能向周緣擴(kuò)散，導(dǎo)致應(yīng)力不斷增加，隆升不斷加快，因此本文選定穩(wěn)定初期(20萬年)模擬結(jié)果展開分析.另外，模擬過程也未考慮地殼增厚導(dǎo)致的拆沉和地幔物質(zhì)的底侵作用，同樣會(huì)對(duì)模擬結(jié)果造成影響.實(shí)際松潘—甘孜地塊地殼和上地幔頂部物質(zhì)，受四川盆地阻擋，當(dāng)物質(zhì)匯聚到一定程度后會(huì)被迫向兩側(cè)運(yùn)動(dòng)，同時(shí)也會(huì)發(fā)生地殼拆沉和地幔物質(zhì)底侵作用，實(shí)現(xiàn)物質(zhì)運(yùn)移的動(dòng)態(tài)平衡(圖11).本文模擬結(jié)果隨構(gòu)造加載不斷進(jìn)行，隆升速率逐漸增大，也是對(duì)青藏高原東緣物質(zhì)匯聚、擴(kuò)展，存在穩(wěn)定開放物質(zhì)流動(dòng)系統(tǒng)的側(cè)面印證.此外，構(gòu)造加載50萬年和100萬年時(shí)，模型西邊界中下地殼和上地幔出現(xiàn)了不同程度的塌陷(圖10b、10d黑色區(qū)域)，也是因?yàn)槲镔|(zhì)東移，在本文模型中缺少相應(yīng)的物質(zhì)補(bǔ)充導(dǎo)致.

圖10 Model 5構(gòu)造加載50萬年和100萬年水平向和垂向隆升速率結(jié)果(a) 50萬年水平向速率； (b) 50萬年垂向速率； (c) 100萬年水平向速率； (d) 100萬年垂向速率.Fig.10 Model 5 simulation results of horizontal and uplift rate at tectonic loading of 0.5 Ma and 1 Ma(a) Horizontal rate at 0.5 Ma; (b) Uplift rate at 0.5 Ma; (c) Horizontal rate at 1 Ma; (d) Uplift rate at 1 Ma.

5 結(jié)論

本文基于二維有限元接觸模型，在考慮地形起伏、地殼結(jié)構(gòu)和密度差異、重力及巖石長期載荷蠕變作用的前提下，依據(jù)1991—2016年GPS觀測結(jié)果，通過改變低阻低速層蠕變參數(shù)(是否考慮低阻低速層的存在)、模型西邊界和底邊界加載條件，共計(jì)開展了11項(xiàng)數(shù)值模擬實(shí)驗(yàn)，將模擬結(jié)果與1960—2010年長時(shí)段水準(zhǔn)觀測結(jié)果進(jìn)行對(duì)比，結(jié)合低溫?zé)崮甏鷮W(xué)隆升及剝蝕速率相關(guān)研究成果，討論了松潘—甘孜地塊地殼至上地幔頂部現(xiàn)今變形及物質(zhì)運(yùn)移模式.主要得到以下結(jié)論：

(1) 不考慮剝蝕作用的前提下，松潘—甘孜地塊長期穩(wěn)定最大隆升速率為1.4～1.5 mm·a-1可能更為客觀.

(2) 松潘—甘孜塊體內(nèi)低阻低速層可以作為中下地殼和上地幔物質(zhì)運(yùn)移的一個(gè)滑移面，構(gòu)成上地殼與中下地殼的“解耦”帶，促進(jìn)松潘—甘孜塊體整體隆升，但低阻低速層對(duì)隆升影響有限，不足以形成如此規(guī)模的強(qiáng)烈的地表抬升.

圖11 青藏高原東緣物質(zhì)運(yùn)移模式示意圖Ⅰ 松潘—甘孜地塊； Ⅱ 四川盆地； Ⅲ 西秦嶺塊體； Ⅳ 川滇“菱形塊體”. ① 龍門山斷裂；② 龍日壩斷裂； ③ 鮮水河斷裂；④ 東昆侖斷裂； ⑤ 岷江斷裂； ⑥ 虎牙斷裂.Fig.11 The pattern of material transport mode in the eastern margin of the Tibet Plateau Ⅰ Songpan-Garzê block; Ⅱ Sichuan basin; Ⅲ West Qinling block; Ⅳ Sichuan-Yunnan rhombic block； ① Longmenshan fault; ② Longriba fault; ③ Xianshuihe fault; ④ East Kunlun fault; ⑤Minjiang fault; ⑥Huya fault.

(3) 考慮低阻低速層作用，模型西邊界加載速率由3 mm·a-1隨深度線性增加1.8倍的前提下，地表隆升速率與長時(shí)段水準(zhǔn)觀測結(jié)果及現(xiàn)今地貌結(jié)合低溫?zé)崮甏鷮W(xué)獲得的隆升速率最為吻合.在這一動(dòng)力條件下，中下地殼和上地幔物質(zhì)受四川盆地強(qiáng)烈阻擋，速率快速減小區(qū)域主要分布于松潘—甘孜地塊內(nèi)部(模型50～150 km范圍)，松潘—甘孜地塊東緣和龍門山斷裂帶速率減小已不顯著，故地表隆升中心位于川西高原內(nèi)而非龍門山斷裂帶.

(4) 青藏高原東緣物質(zhì)匯聚到一定程度后，會(huì)被迫向兩側(cè)運(yùn)動(dòng)，從而實(shí)現(xiàn)物質(zhì)運(yùn)移的動(dòng)態(tài)平衡.

本文模擬過程存在以下不足：(1)重力加載模型穩(wěn)定后，松潘—甘孜地塊相對(duì)四川盆地地形高差減至約3500 m，與4000 m實(shí)際地形高差相差約500 m，低估了重力勢的影響；(2)模擬過程未考慮地殼增厚導(dǎo)致的拆沉和地幔物質(zhì)的底侵作用.基于獲得的松潘—甘孜地塊至四川盆地地殼和上地幔頂部現(xiàn)今變形及物質(zhì)運(yùn)移模式，結(jié)合龍門山斷裂帶周緣精細(xì)速度圖像結(jié)果，構(gòu)建三維有限元模型，研究汶川MS8.0地震NW走滑型余震帶(米亞羅斷裂)動(dòng)力學(xué)機(jī)制，揭示青藏高原東緣擠壓動(dòng)力學(xué)背景下，走滑型強(qiáng)震成因，可以作為今后工作的一個(gè)研究方向.

致謝圖件由GMT繪制而成，GPS數(shù)據(jù)采用了中國地震局地質(zhì)研究所王敏研究員計(jì)算結(jié)果，水準(zhǔn)觀測結(jié)果采用了中國陸地現(xiàn)代垂直形變圖集的編制與資料整編項(xiàng)目組計(jì)算結(jié)果，在此一并表示感謝！