鄭群凡，石耀霖

(中國科學院大學地球與行星科學學院 中國科學院計算地球動力學重點實驗室， 北京 100049)

印度板塊和歐亞板塊的碰撞是地球上最典型的造山事件之一，形成了喜馬拉雅造山帶和廣闊的青藏高原。其初始碰撞發(fā)生在大約60 Ma前，持續(xù)的碰撞造成了亞洲一側1 000 km以上的地殼縮短量[1]。一直以來，對青藏高原的形成和生長的動力學討論大多側重于碰撞后物質的運移、地殼變形，例如，柔性下地殼流[2-3]、構造逃逸[4]、地殼一致性加厚[5]、亞洲大陸地殼結構橫向不均一性[6]等方面。但對于導致印度—歐亞板塊會聚和青藏高原隆升的動力來源問題始終沒有得到很好的解決。

傳統(tǒng)的板塊構造理論及威爾遜旋回的觀點認為大陸巖石圈在進入海溝之后，地殼不斷加厚，隨著陸陸碰撞持續(xù)進行，俯沖板片發(fā)生斷離或巖石圈地幔拆離，驅動板塊運動的俯沖板片拉力消失，地殼增厚也隨之停止[7-8]。然而，印度—歐亞板塊碰撞后，2個大陸板塊之間快速的會聚已經持續(xù)了數千萬年，地殼持續(xù)增厚形成了青藏高原[9]，但印度板塊仍持續(xù)向北運動?？梢娪《却箨懗掷m(xù)向北運動的驅動力不但克服了長時間尺度上碰撞系統(tǒng)的巨大阻力，并足以使青藏高原持續(xù)地抬升與擴展。

近年來的研究強調由地幔對流造成的巖石圈底部的地幔拖曳力對區(qū)域構造演化的影響。Jolivet等[10]通過比較50 Ma以來的巖石圈板塊的運動軌跡和地震層析成像、各向異性模型，認為軟流圈地幔拖曳印度板塊向北運動，形成喜馬拉雅山脈和青藏高原，并且控制東亞大陸變形。Zou等[11]提出與向北運動的印度板塊耦合的地幔流存在的地震學證據，并認為這種地幔流可能帶動了GPS觀測到的塊體擠出現(xiàn)象。前人也做了有關這一問題的數值模擬研究。例如Liu等[12]的二維模擬結果表明，類似于地殼塊體的側向擠出，印度—歐亞碰撞造成了顯著的大尺度側向地幔流，且這種地幔擠出具有遠程效應，可以影響到中國東部。祝愛玉等[13]對青藏高原東北緣的小尺度地幔對流進行數值模擬，認為東北緣區(qū)域承受了強烈的巖石圈底部水平拖曳力。

那么在實際三維地球中，東亞區(qū)域是否在巖石圈底部存在施加拖曳力的大尺度的地幔流？如果存在，這種拖曳力在印度—歐亞板塊持續(xù)會聚中起到的作用是什么？對青藏高原的構造演化有什么樣的影響？基于這些問題，本文建立了三維全球有限元模型，基于地幔動力學探討印度板塊持續(xù)向北運動的驅動力來源，嘗試給出印度—歐亞會聚的一階控制因素。

1 計算方法

1.1 控制方程及數值實現(xiàn)

作為一階近似，巖石圈和地?？杀灰暈槿S球形幾何中的不可壓縮黏性流體。質量、動量和能量的守恒方程由下式給出：

?iui=0,

(1)

?iσij=-Δρgi,

(2)

(3)

其中：u為速度,σij為全應力張量，Δρ為密度異常，g為重力，T為溫度，κ為熱擴散系數，t為時間。

忽略化學組分的影響，由巖石的熱膨脹性產生的密度異常可以寫為

Δρ=-αρ0(T-T0),

(4)

其中：α,ρ0和T0分別為熱膨脹系數、參考密度和參考溫度。

全應力張量可以表示為

σij=τij-pδij,

(5)

其中：τij,p和δij分別表示偏應力、動壓力和Kroneckerδ函數。

偏應力張量為

(6)

其中μ為黏度。

為建立動力學模型，我們使用開源的大規(guī)模并行有限元代碼ASPECT[14]。ASPECT的核心代碼關注于方程的實現(xiàn)，并與包括DEAL.II、TRILINOS及 p4est在內的庫進行通信。在有限元建模中，對流動變量的空間近似可以有不同的選擇。在這里，對速度和溫度使用二階元，但對壓力使用線性元，以滿足LBB(Ladyzhenkaya、Babouska和Brezzi)條件[15-16]。對于時間離散，由于斯托克斯方程不含時間導數，因此斯托克斯方程可以在每個時間步限制能量方程。這類方程常采用IMPES格式求解，即顯式溫度解和隱式斯托克斯解。能量方程中的時間導數由BDF-2(二階向后差分法)時間步取代[17]。ASPECT使用CFL(Courant-Friedrichs-Lewy)條件[18]來選擇時間步長。有關更多技術細節(jié)，可以參考ASPECT手冊(https:∥aspect.geodynamics.org/)。

1.2 模型設置

為尋求更接近實際地球的初始條件來探究地球演化歷史，數值模擬通常會采用地震層析成像反映地幔不均一性，這種方法已在前人的研究中被廣泛采用[19-23]。假設地震速度異常與密度和溫度變化有關，就可以從地震層析成像中得到地幔密度和溫度的不均一性。我們模型初始溫度場考慮由全球地震層析成像模型S20RTS[24]換算成的溫度異常。在換算中，首先用關系式(δρ/ρ)/(δvs/vs)=0.15將地震波速異常轉換成密度異常[21,25-28]，然后用關系式δT=(δρ/ρ)/α[29]把密度異常轉換為溫度異常。

以往使用二維矩形模型或者三維六面體模型時，在中側邊界及底邊界上施加什么邊界條件是一個難題[12]，往往由于缺乏觀測約束而不得不人為假定。雖然本文只聚焦于青藏高原及鄰區(qū)的地幔運動，但為了減少邊界條件對模擬的影響，以求得更加合理的計算結果，我們的模擬是在三維全球框架下研究區(qū)域問題，計算全球在給定的基于地震層析成像的密度分布下的對流模式。這樣只需要關注地表邊界條件。使用Seton等[30]重構的全球板塊運動模型作為模型的地表邊界條件，利用Gplates[31]把該模型生成為計算所需的1°×1°分辨率格式。

在我們的三維全球模型中，黏度結構為分層常黏度。巖石圈(0～100 km)黏度為1023Pa·s，上地幔(100～660 km)黏度為5×1020Pa·s，下地幔(660～2 890 km)黏度為 3×1022Pa·s[26,32]。模型計算中采用的物性參數如表1所示。地球表面溫度邊界條件取為273 K。核幔邊界溫度約為3 000～4 000 K[33]。上地幔和下地幔的絕熱溫度梯度分別為0.4～0.5 K/km和0.3 K/km[34]。模型中的溫度梯度是絕熱的，溫度邊界條件應考慮絕熱效應。模擬中使用2 600 K作為核幔邊界溫度值，這是真實溫度減去絕熱加熱量的近似值。

表1 模型計算中采用的物性參數Table 1 Material properties used in the numerical experiment

2 計算結果

基于上述計算方法，將研究區(qū)域(圖1)從全球模型中截出，主要展示研究區(qū)域水平流場及4條典型剖面上的結果。這4條剖面位置顯示在圖1。

2.1 軟流圈流場預測

用180 km深度的地幔流場近似代表軟流圈流場的一階特征(圖2(b))。為對比起見，圖2(a)顯示了絕對板塊運動給出的速度[30]，它與文獻中經常給出的以歐亞板塊為參照的GPS速度不同。計算得出的東亞大陸下方軟流圈地幔水平流場(圖2(b))和地表絕對板塊運動速度場方向(圖2(a))大體相符。比如，在板塊內部兩者總體速度基本一致，在板塊邊界也都存在明顯的突變。但在靠近板塊邊緣的局部位置存在一定的差別。如，在青藏高原東南緣及其以南，前者存在一定的轉變過渡，而后者沒有。

從軟流圈地幔水平流場(圖2(b))特征來看，中國大陸下方的大部分地區(qū)軟流圈地幔流動方向大致為東南向。只是在出現(xiàn)類似于巖石圈塊體擠出現(xiàn)象的青藏高原東南緣地區(qū)，軟流圈地幔出現(xiàn)了大尺度順時針的旋轉趨勢，由東南向轉為南向。印度板塊下方北東向的地幔流方向在印度—歐亞板塊邊界附近發(fā)生突然改變。同樣地，西太平洋下北西向的地幔流在洋陸俯沖板塊邊界位置與東南向的來自中國大陸東部的地幔流相遇。據此推定，180 km軟流圈流場與板塊相互作用存在明顯的關系。如果這樣，中新生代以來的印藏碰撞、太平洋板塊俯沖、北側蒙古—鄂霍茨克板塊以及南部緬甸板塊的俯沖都有可能使得中國大陸下方的地幔水平流場發(fā)生系統(tǒng)性偏移。

圖1 研究區(qū)域及剖面位置Fig.1 Study region and locations of cross-sections

圖2 重構的地表板塊運動速度(據Seton等[30])及180 km深度地幔水平流場Fig.2 Reconstructed plate motion from Seton et al.[30] and horizontal velocity field at 180 km depth

與180 km深度的物質流動不同，計算得到的180 km深度以下的上地幔物質流動與地表物質運動存在明顯差異，有關情況從下面的剖面圖來展示。

2.2 剖面結果

剖面AA′主要穿過印度板塊(圖1)，該剖面速度場見圖3(a)。在剖面圖中，箭頭表示平行于剖面方向上的速度分量，顏色則代表垂直于剖面方向上的速度。藍色表示速度方向垂直于剖面向外(近似南南西向)，紅色表示速度方向垂直于剖面向內(近似北北東向)。計算結果顯示在剖面穿過的大部分地區(qū)，巖石圈與軟流圈流場速度都向北北東方向(粉色和紅色)，然而軟流圈速度明顯大于上覆巖石圈的速度，此時軟流圈對上覆巖石圈肯定存在拖曳作用。

在剖面中段和東段，印度板塊向東俯沖至緬甸之下。剖面AA′ (圖3(a))的最東段位于青藏高原東南緣，喜馬拉雅東構造結以東，是亞洲板塊的一部分。垂直于剖面的速度場顯示最東段巖石圈和軟流圈速度方向為南南西方向(藍色)。這顯示亞洲板塊與印度板塊在緬甸附近的板塊邊界東西兩側的深部運動存在很大的差異。

圖3 剖面的速度場Fig.3 Flow profiles of cross-sections

東西向的剖面BB′大部分位于青藏高原，計算得到的速度場如圖3(b)所示，圖中的紅色代表向北的速度，藍色代表向南的速度。剖面西部位于印度板塊的一段仍顯示出與剖面AA′ (圖3(a))相同的特征，即軟流圈向北的流動速度明顯大于其上覆巖石圈。而剖面BB′中位于青藏高原的中段和東段，則顯示出與西段相反的流動方向，即巖石圈與軟流圈均向南流動。此時軟流圈向南的流動速度也大于其上覆巖石圈。剖面AA′與BB′位置相近，但分別位于喜馬拉雅弧的弧南與弧北，印度—歐亞碰撞帶之下復雜的結構使得2個剖面的結果出現(xiàn)了很大的差異。剖面AA′和BB′向下拖曳的地幔流的位置主要是受到印度與歐亞板塊碰撞帶的影響。剖面AA′位于喜馬拉雅弧以南，剖面中段切過主邊界逆沖斷裂帶，由于印度板塊沿歐亞板塊南緣的俯沖作用，帶動了向下的地幔流。而剖面BB′的西段是印度—歐亞碰撞帶的位置，中段位于拉薩地塊內部，因此向下拖曳的地幔流主要出現(xiàn)在剖面西段。

橫穿青藏高原的東西向剖面CC′速度場計算結果為圖3(c)，圖中的紅色代表向北的速度，藍色代表向南的速度。垂直于剖面的速度結果顯示這一區(qū)域軟流圈的速度仍舊大于其上覆巖石圈，軟流圈對巖石圈向南的拖曳作用十分明顯。平行于剖面的速度場在剖面西段出現(xiàn)地幔環(huán)流，可能與印度板塊向歐亞板塊之下的俯沖作用有關。

剖面DD′由南向北穿越印度、青藏高原、塔里木盆地、天山等區(qū)域。在垂直于剖面DD′的速度場中(圖3(d))，藍色代表近似向東偏南的速度，紅色反之向西偏北?？梢钥吹狡拭鍰D′同樣出現(xiàn)軟流圈對上覆巖石圈向東的拖曳作用。平行剖面的速度場顯示出印度板塊向歐亞板塊之下的深俯沖，并且俯沖板片可能出現(xiàn)向南翻轉的幾何形態(tài)。印度—歐亞碰撞帶之下的俯沖造成了遠場效應，帶動兩側的地幔物質流向俯沖帶位置，驅動小規(guī)模的地幔環(huán)流。該計算結果顯示，地幔環(huán)流有利于印度巖石圈向青藏高原下的俯沖和高原巖石圈的拆沉。

計算結果反映青藏高原下有可能存在印度巖石圈的俯沖和高原巖石圈的拆沉以及拆沉誘發(fā)的鄰側地幔上涌，這與青藏高原北部地區(qū)廣泛發(fā)育新生代的火山巖, 并與普遍具有低Pn波速、缺失Sn波的特點相一致。但地震學中觀察到的俯沖的角度以及板片是否斷離等問題在我們的模型中難以體現(xiàn)。

3 討論

以上模擬結果顯示了青藏高原及周邊地幔對流的特征。印度板塊在靠近喜馬拉雅山脈的板塊邊界處軟流層向北的運動速度大于其上覆巖石圈的運動速度(圖3(a))。這說明向北運動印度軟流層拖曳著印度巖石圈向北運動。除負浮力作用下俯沖，這個持續(xù)的驅動力導致印度大陸與歐亞大陸接觸并發(fā)生碰撞。即使在喜馬拉雅山脈和世界最大的青藏高原形成以后，仍然能夠克服高原巨大位能形成的阻力，持續(xù)拖拽印度板塊向北運動，并導致持續(xù)的碰撞和高原隆升生長。

俯沖的印度板片巖石圈驅動地幔在俯沖位置下面形成了對流環(huán)。從更深部位來看，印度板塊軟流層和巖石圈向北運動的同時，在印度板塊軟流層(400 km深度)以下的深層地幔物質則產生了向南的回流(圖3(a)、3(b)、3(d))。

俯沖的印度地幔也帶動了北方青藏高原之下軟流圈的運動，造成碰撞帶以北的軟流圈向南的運動(圖3(b)、3(c))。人們討論印度板塊和歐亞板塊碰撞運動時，往往相對于歐亞板塊為參考，即假定歐亞板塊靜止不動來進行討論。但利用這種假設將無法考慮印度板塊運動對青藏高原運動的拖拽影響，這在今后的研究中值得引起重視。計算顯示，印度巖石圈向歐亞板塊之下的俯沖(圖3(d))影響歐亞板塊青藏高原遠場的地幔流，并可能在青藏高原下面引起小尺度的地幔擾動和上下流動。Huangfu等[35]小尺度的模擬表明，青藏高原之下發(fā)生巖石圈地幔的拆離和軟流層地幔熱物質的上涌可以影響到高原的溫度結構、地震波速和衰減，以及巖漿活動。但本文由于網格尺度和全球層析成像的分辨能力有限，無法精確反映這些小尺度的擾動。另外值得注意的是，與印度板塊之下不同，青藏高原之下并沒有造成上地幔深部物質的回流(圖3(c)、3(d))，400～670 km的地幔沒有大規(guī)模向北的回流運動。

圖3(a)～3(c)都顯示了青藏高原之下軟流層存在向東的運動，而上地幔深部存在向西的回流，這些可能與太平洋板塊和菲律賓板塊向歐亞板塊下的俯沖有關。由于這些海洋板塊的俯沖規(guī)模更大、俯沖海洋板片負浮力也更大，因此影響的范圍也更大，中國的西部青藏高原深部的地幔運動有可能受到其影響。前人的研究[12]曾提出，印度—歐亞碰撞驅動了大尺度的側向地幔擠出，其遠場效應可能導致中國東部廣泛的軟流圈上涌、裂谷作用和新生代火山作用。我們的三維模擬說明，這種流動的確可能存在，但三維的流動圖形遠比二維復雜。印度板塊下(圖3(a))的上地幔深部也有這種向西的回流，可能與印度板塊在緬甸之下向東俯沖有關。這種俯沖還造成緬甸弧東部軟流層流動方向轉為向東，而不再是川滇軟流層和巖石圈的向南運動(圖2(b))，軟流層運動方向的這種轉變，恰好可以解釋SKS波快波方向從川滇的南北向轉為滇南和緬甸的東西向所體現(xiàn)的各向異性特征。

目前，已有不少地震學研究提出青藏高原下印度俯沖的板片可能發(fā)生撕裂[36-38]。在圖3(d)中，雖然400 km以淺的范圍都顯示出東偏南的速度(藍色)，然而在高原內部軟流圈也存在一定的速度差異，比如在青藏高原中南部的軟流圈速度明顯小于兩側的速度。這種地幔經向運動差異可能為喜馬拉雅俯沖板片的撕裂提供了一定力學條件。

理想的模型應該基于合適的模型參數、精細的橫向和縱向不均勻性等，在更真實的自由邊界下的計算才能夠產生實際觀測到的板塊運動。但人們目前對深部的認識還有限，我們也只能把地震層析成像確定的速度異常換算為初始模型的溫度異常，把地球表面板塊運動作為邊界條件。在這種前提的模擬中，如果計算得到的軟流層運動速度低于地表運動速度，不排除可能是由于人為施加的速度邊界條件的影響，軟流層被地表巖石圈拖曳而動。但是在本文的模型中，很多結果都顯示出軟流層速度快于其上覆巖石圈的運動速度，這種運動特征不可能是由于指定的邊界條件造成的，而是真實地幔對流物理過程的表現(xiàn)。

我們在模擬中也嘗試改變一些介質參數，觀察它們對計算結果的影響。例如上下地幔黏滯系數之比，巖石圈蓋層的存在的影響，但這些因素對結果總體特征沒有大的影響，沒有改變地幔流場的一階特征。在本文中，討論大尺度地幔拖曳驅動力的存在和一階作用，但不排除印度—歐亞板塊碰撞、西太平洋俯沖帶、爪哇海溝俯沖帶或者來自非洲下的超級地幔柱等因素分別對這種東亞地幔拖曳力的影響。

本研究選擇的是加州理工學院發(fā)展的S20RTS面波成像模型，該模型水平采用20階球諧函數，徑向則采用3次樣條曲線[24]。依據的資料不僅包括體波和瑞利波測量數據集，還包括對自由震蕩分裂的分析，是一個優(yōu)秀的模型。但是國際上也還有其他層析成像模型，例如S362ANI[39]、SAW642AN[40]、TX2008[41]以及 HMSL-S[42]等模型。如果更換其他模型，對對流計算結果會有什么影響是需要進一步研究的問題。由于各種模型都反映了太平洋、菲律賓等海洋板塊俯沖形成的高速板片，以及印度板塊巖石圈地幔俯沖形成的高速區(qū)域，因此本文的基本結論應該依然成立。

4 結論

本文通過三維全球有限元數值模擬，避免了局部模型模擬中地幔深部側邊界速度邊界條件不確定性的影響，在全球規(guī)模地幔對流的大背景下，計算東亞區(qū)域地幔流場，從中分析新生代以來的印度板塊持續(xù)向北運動、印度—歐亞持續(xù)會聚的驅動力，并得出以下認識：

1)軟流層在印度巖石圈底部施加的拖曳力可能是印度板塊持續(xù)運動和造成青藏高原隆起的驅動力的主要來源。印度巖石圈向歐亞板塊之下的俯沖帶動了地幔軟流層的運動。

2)印藏碰撞具有遠場效應，造成喜馬拉雅弧南北兩側的地幔物質流向俯沖帶的位置并可能引起大規(guī)模地幔環(huán)流，且有可能在青藏高原下方誘發(fā)小規(guī)模的地幔升降流動。

3)青藏高原東南緣地幔流出現(xiàn)“擠出”現(xiàn)象，川滇地區(qū)向南運動的軟流圈在緬甸弧東部轉為向東。

4)青藏高原中南部的軟流圈運動存在速度差，這可能為印度俯沖板片的撕裂提供了條件。