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

萬永魁, 劉峽, 劉瑞豐, 張揚, 沈小七, 鄭智江

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

0 引言

松潘—甘孜地塊隆升機制主要有側(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的最大垂直位移，對應(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è)殼幔速度差異顯著，相對四川盆地，松潘—甘孜地塊普遍具有低速、高導(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)成上地殼與中下地殼“解耦”運動的滑脫層或剪切帶.龍門山后山、中央、前山和山前隱伏斷裂向下延伸傾角逐漸變緩，最終收斂于該剪切帶(Hubbard and Shaw, 2009; Yu et al., 2010; Wan et al., 2017；許志琴等，2007；張培震，2008).因此，目前對松潘—甘孜地塊隆升之爭，主要集中在僅受控于低阻低速層“解耦”還是同時受殼幔流共同作用這兩種主流觀點上.

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

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

1 模型

1.1 模型構(gòu)建

青藏高原東緣由西向東依次為高海拔強烈變形變質(zhì)作用的松潘—甘孜褶皺帶、地勢起伏劇烈的龍門山?jīng)_斷帶和低海拔弱變形的四川盆地(劉樹根等，2019).晚三疊世以來，受印度板塊與歐亞大陸強烈碰撞及中國陸地主體拼合的影響，松潘—甘孜地塊和龍門山造山帶發(fā)生強烈變形和沖斷隆升(劉樹根等，2009)，在東西寬約30～50 km范圍內(nèi)，西側(cè)松潘—甘孜地塊與東側(cè)四川盆地形成高達約4 km地形差異，成為整個青藏高原地形梯度最大地區(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)，具有強度大、能夠積累高密度彈性應(yīng)變能的特性(張培震等，2008).阿壩—雙流人工地震探測剖面顯示龍門山斷裂帶兩側(cè)地殼厚度變化顯著，松潘—甘孜地塊至四川盆地寬約80 km范圍，地殼厚度由約60 km迅速遞減至約45 km(朱介壽，2008).

依據(jù)阿壩—雙流人工地震探測剖面(圖1黑色粗實線)，參考前人數(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)以強度較高的彭灌雜巖、寶興雜巖及轎子頂雜巖為核心，其彈性模量為四川盆地的1.2倍.松潘—甘孜地塊0～4 km范圍代表高海拔地形，埋深20～25 km范圍為低阻低速層(圖2a黃色區(qū)域)，用于模擬上地殼與中下地殼的“解耦”.模型整體構(gòu)架見圖2a.

圖1 研究區(qū)水準觀測路線、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è)置

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

(1)

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

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

(3)

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

(4)

本文在應(yīng)變率設(shè)定為10-14s-1前提下(石耀霖和曹建玲，2008)，計算得到松潘—甘孜地塊和四川盆地巖石圈各層蠕變系數(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紅色部位).采用三角形單元對模型進行網(wǎng)格劃分，網(wǎng)格尺寸約1 km，網(wǎng)格劃分結(jié)果見圖2b，單元總數(shù)31749，節(jié)點總數(shù)32116個.

圖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 加載計算

2.1 加載

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

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

加載分兩步進行：(1)重力加載；(2)位移(構(gòu)造)加載.重力加載時，模型東、西邊界(x=300及x=0 km)水平向固定、垂向自由，底邊界(y=-100 km)垂向固定、水平向自由，施加重力后維持1 Ma自由變形.位移加載是在重力加載的基礎(chǔ)上完成的，一方面更接近實際，另一方面在構(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ù)計算獲得蠕變參數(shù))，體現(xiàn)低阻低速層“解耦”對高原隆升的影響(Model 2).參考胡幸平等(2012)由地表至深100 km處，加載速率線性增加5倍，所得龍門山地區(qū)理論震源機制解與實際震源機制解最為吻合，劉峽等(2010)在華北地區(qū)現(xiàn)今地殼運動模擬研究中設(shè)定由地表至深100 km處，加載速率線性增加1.2倍，模擬結(jié)果與GPS觀測結(jié)果最為一致，本文通過模型西邊界加載速率隨深度線性增加(1.2、1.5、1.8、2.0、2.5倍或3.0倍)，體現(xiàn)巖石圈水平向運動隨深度增加對高原隆升的影響(Model 3—Model 8).最后，基于Model 5的模擬結(jié)果，通過改變底邊界位移加載條件(類比地幔拖曳力)，體現(xiàn)地幔對流拖曳力對高原隆升的影響(Model 9—Model 10).此外，不考慮Model 5中低阻低速層的作用，模擬高原隆升速率(Model 11)，通過分析兩模型隆升速率差異，進一步研究低阻低速層“解耦”對高原隆升的作用.

2.2 計算

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

3 模擬結(jié)果

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

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

圖3 重力加載后不同時段節(jié)點平均運動速率(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 重力加載后不同時間地表相對高程(a)和相應(yīng)時間段地表垂向速率(b)；構(gòu)造加載后Model 5在不同摩擦系數(shù)下節(jié)點174(c)、節(jié)點21913(d)相對高程變化；構(gòu)造加載20萬年Model 5在不同摩擦系數(shù)下地表節(jié)點水平壓縮速率(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é)點隆升位移隨位移加載時間變化曲線Fig.5 The graph of uplift displacement of surface nodes with tectonic loading time of Model 5

圖6 Model1-Model11構(gòu)造加載20萬年時地表水平向速率及GPS擬合結(jié)果(a)、(c)和隆升速率及水準結(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)定性分析

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

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

4 討論

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

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

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

4.2 低阻低速層對高原隆升的影響

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

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

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

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)部，水平向運動在低阻低速層發(fā)生了“解耦”，上地殼運動速率小于中下地殼.垂向速率與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ū)域基本擴展至整個松潘—甘孜地塊，強烈隆升區(qū)域范圍顯著擴展；(3)Model 2隆升中心位于模型150 km附近，Model 5向西轉(zhuǎn)移至模型100 km附近，與長時段水準觀測結(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)幅度顯著增強，中下地殼和上地幔水平向運動速率快速減小區(qū)域集中在模型50～150 km處，松潘—甘孜東緣和龍門山斷裂帶(模型150～200 km)，速率減小已不顯著，這與水準觀測隆升中心并非在龍門山斷裂帶，而是在其西側(cè)約100 km處的川西高原，即模型50～150 km范圍相吻合.中下地殼和上地幔物質(zhì)水平向運動速率快速減小與四川盆地的阻擋固然有關(guān)，但與地形起伏、地殼結(jié)構(gòu)和密度差異、巖石圈流變及地幔對流拖曳力可能均有關(guān)聯(liá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對模型底邊界進行了位移加載，用于類比地幔拖曳力，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地表隆升速率與水準觀測結(jié)果均存有了一定差異(圖6d)，垂向速率隨深度變化結(jié)果與Model 5模擬結(jié)果差異也較為顯著(圖7f、7h、7j).上述結(jié)果反映出隆升速率和隆升中心的位置分布對底邊界加載條件較為敏感.

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

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

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

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

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

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

圖11 青藏高原東緣物質(zhì)運移模式示意圖Ⅰ 松潘—甘孜地塊； Ⅱ 四川盆地； Ⅲ 西秦嶺塊體； Ⅳ 川滇“菱形塊體”. ① 龍門山斷裂；② 龍日壩斷裂； ③ 鮮水河斷裂；④ 東昆侖斷裂； ⑤ 岷江斷裂； ⑥ 虎牙斷裂.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倍的前提下，地表隆升速率與長時段水準觀測結(jié)果及現(xiàn)今地貌結(jié)合低溫?zé)崮甏鷮W(xué)獲得的隆升速率最為吻合.在這一動力條件下，中下地殼和上地幔物質(zhì)受四川盆地強烈阻擋，速率快速減小區(qū)域主要分布于松潘—甘孜地塊內(nèi)部(模型50～150 km范圍)，松潘—甘孜地塊東緣和龍門山斷裂帶速率減小已不顯著，故地表隆升中心位于川西高原內(nèi)而非龍門山斷裂帶.

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

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

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