青藏高原通道流模型動力環(huán)境的數(shù)值模擬

楊 輝，滕吉文，皮嬌龍

中國科學院地質(zhì)與地球物理研究所，北京 100029

1 引 言

巖石圈具有一定的流變性質(zhì)，這是控制巖石圈構造變形的主導因素之一，也是探討巖石圈動力學的基礎.因此，從20世紀70年代開始，許多學者便開始對巖石圈流變性質(zhì)研究，對大陸構造運動的動力條件及大陸在其作用下的運動和形變場進行研究，以探討地球深部各圈層及地表各地質(zhì)單元塊體的規(guī)律性運動和形變過程.但利用礦物和巖石的實驗規(guī)律來定量或是半定量地外延去定義巖石圈分層結(jié)構的流變參數(shù)，從準確性及精度上尚遠遠不能滿足大陸動力學及板塊運動學的要求，這也使得在對不同區(qū)域基于地質(zhì)、地球物理和地球化學等觀測數(shù)據(jù)而得出的構造演化的動力學過程解釋中出現(xiàn)不同的，甚至相互矛盾的解釋［1］.

當前，對巖石圈內(nèi)部的流變機制認識還相當有限，主要是基于地震波場和地熱場方面的資料，而這些地球物理量的真實測定及解釋還存在很大的不確定因素，且地球介質(zhì)的屬性并不全等同于物理上的流變性，因而準確地說應當為“視流變性”.顯然在這種尚不能給出地球深部介質(zhì)真實黏滯系數(shù)和全基于物理上的流變機制狀況下，采用各種假設條件去試圖解釋巖石圈動力學機制，這會給理解巖石圈流變結(jié)構增加很大的困難，乃至在解釋方面引發(fā)很大爭議.為此，須在這一基點上對地球深部介質(zhì)屬性進行探索.

以青藏高原為例，對這個由印度板塊與歐亞板塊因碰撞而形成的地球上最高最大的高原，國內(nèi)外已有不少人涉足于這一方面的研究，但是對其從表層到深部的三維變形機制，至今仍眾說紛紜.近年來，爭論的焦點主要圍繞著兩個端元模型：一是巖石圈連續(xù)變形模型［2－3］，認為高原在演化發(fā)展的過程中，變形從地殼到地幔是連貫且耦合的；一是通道流（channel flow）模型［4－7］，認為變形是隨深度變化的，是由下方的軟流圈流體驅(qū)動或是由上地殼拖曳運動.

可以說，如何對研究區(qū)動力學模型進行簡化，特別是對邊界條件的正確理解，會給動力學響應過程的合理解釋帶來很多方面的影響.我們已清晰地認識到，由于實驗手段和方法等方面的限制，地球科學中對物理學中流變學的借用，與物理學中的流變學屬性相比并非同一個量級，而同時又因為當前尚缺乏一個對地球深淺部結(jié)構介質(zhì)屬性黏度的真實物理量的直接測定，這不僅限制了對地球動力學和大陸動力學認識的深化，更帶來了對地球深部介質(zhì)流變性量化的判斷方面的困難.

因此，本文基于青藏高原已經(jīng)開展的深部地球物理探測工作所得到的一些規(guī)律性結(jié)果，針對當前地球科學界流行的通道流模型的一些理論或假設基礎，對青藏高原下方物質(zhì)流動所需要的邊界條件開展了數(shù)值模擬工作，以探討該現(xiàn)象存在的一些必要條件.

2 問題的提出

自Kuhn［8］于1979年提出了通道流的原始概念以來，三十年間，通道流的模型得到了迅速的發(fā)展.通道流模型認為以青藏高原為例的一些造山帶區(qū)域，相對于較為剛性的上地殼和巖石圈上地幔，處于中間層的中、下地殼作為強度較低的塑性流變層，具有更強的流動性，因而在高原厚地殼的重力差異驅(qū)動或是上地幔的拖曳作用下易發(fā)生流動而形成通道流［9］.最初，通道流的模型運用于對造山帶腹地存在高階變質(zhì)巖的地帶進行構造解釋，Grujic［10］首先將通道流模型（Poiseuille和Couette流）力學機制用于喜馬拉雅造山帶腹地的高喜馬拉雅結(jié)晶巖系（Greater Himalayan Sequence，GHS）的韌性擠出構造演化模式中.高喜馬拉雅結(jié)晶巖系位于北部的藏南拆 沉 系 （South Tibetan Detachment system，STD）及南部的主中央斷裂帶（Main Central Thrust，MCT）之間，出露大范圍的前震旦紀變質(zhì)巖，一般由高級變質(zhì)巖組成，屬于印度板塊北緣結(jié)晶基底的組分，剪切變形的應變場特征劇烈.藏南拆沉系為北傾的低角度正斷層系，而主中央斷裂帶為北傾的逆沖斷裂帶，兩條斷裂系近于平行，具有相反的剪切特征，活動特征具有明顯的同期性，約發(fā)生于22Ma左右，在中中新世中止［11］.通道流模型作為其中一種可能的演化方式，認為位于藏南拆沉系和主中央斷裂系之間北傾的高喜馬拉雅結(jié)晶巖系是中地殼物質(zhì)向南擠出的條帶區(qū)域，可以說在解釋高喜馬拉雅結(jié)晶巖系的構造演化和應變模式方面得到了很好的應用.

隨后，Nelson［12］在對藏南構造演化研究時首先明確提出了“channel flow”一詞，通道流模型得到了進一步的發(fā)展.期間，青藏高原眾多深部地球物理探測工作，如重力與地表地形、地熱、大地電磁和地震測深等資料解釋對通道流模型的發(fā)展提供了幫助，同時，還開展了相應的數(shù)值模擬工作［7，13－17］.這些工作的展開將通道流模型運用到了更多的構造背景［18］，如：軟流圈倒流（asthenospheric counterflow）；下地殼通道流（lower crustal channels）；造山帶殼內(nèi)通道流（intra－crustal channels）；俯沖帶通道流（subduction channels）；鹽構造通道流（salt tectonics）等，并且，其應用區(qū)域已不僅僅局限于喜馬拉雅造山帶和青藏高原周緣，還應用到如希臘的西奈山（Hellenides）［19］的中地殼韌性變形擠出模式和阿巴拉契亞內(nèi)山麓（Appalachian Inner Piedmont）［20］的通道流模式等這類古老造山帶的構造演化動力學模型中.

由此可見，通道流模型的適用與否，不僅關系到青藏高原隆升的機制問題，還牽涉到全球大陸動力學系統(tǒng)的合理解釋.

在通道流模型的兩種端元模型中，Couette流是由頂部或底部剛性較強層的運動所導致的剪切力系牽連效應使中間的通道流層發(fā)生流動，如重力垮塌引起上部剛性上地殼的橫向運動，地幔流動對下地殼的拖拽運動，以及巖石圈向下俯沖而通過層間耦合驅(qū)動軟弱的中、下地殼的流動，其通道流中速度最大的區(qū)域發(fā)生在頂部或底部；而Poiseuille流的上、下邊界是穩(wěn)定的，由兩側(cè)的壓力梯度差導致通道流層發(fā)生運動，如存在地殼密度差異、地殼厚度差異或是地貌高度差異引起的橫向壓力梯度驅(qū)動殼內(nèi)軟弱層的流動，其最大流動速度位于通道流層的中部.顯然，通道流的存在必須是低黏度的物質(zhì)被高黏度的物質(zhì)所包圍，形成物理意義上的通道.而由于深部地殼真實介質(zhì)屬性難以直接測定，因而低黏度層的特異黏滯系數(shù)主要由數(shù)值模擬工作者來界定，認為只要該層設定為黏滯系數(shù)很低的物質(zhì)組成，便可產(chǎn)生通道流效應.然而，必須清晰地認識到，黏滯系數(shù)的正確給定和理想模型的設定是具有一定的局限性，即必須從實際出發(fā)加以約束，方能提取一個合理的動力學模型.

通道流模型認為低黏度區(qū)域的特異黏滯系數(shù)往往是由部分熔融或含水流體所致，而地球內(nèi)部各圈層的黏滯系數(shù)可以通過多要素約束下的邊界條件進行理論計算，亦可通過高溫高壓實驗、地震學、地磁學、重力學和大地測量學等多種途徑獲得.然而，用不同手段和方法所給出的數(shù)值結(jié)果之間其差異明顯，誤差范圍在最佳情況下亦可達2～3個數(shù)量級，而用固體潮方法得到的黏滯系數(shù)差別更高，甚至可達到1011Pa·s不等（孫和平，個人交流，2010）.

已開展的數(shù)值模擬工作亦存在同樣的問題.Clark和Royden［6］的二維模型很好地重現(xiàn)了通道流的演化，但通道流層內(nèi)物質(zhì)的黏滯系數(shù)設置過低.同樣，Shen等［21］所用三維牛頓黏性流體模擬軟弱的下地殼底部時，對巖石圈地幔的強度設置過高，這些都難以被深部地球物理工作所完全支持.可見，復雜的數(shù)值模型具有很強的特殊性，其中的因果關系亦往往不是輕易就能一目了然辨析而出的.

因此，這樣的模擬所得結(jié)果僅可供討論問題或推斷及猜測認識時的一個參考，若以其作為定量的模型或解決深部物質(zhì)的逼近運移軌跡則尚難以令人信服.

3 青藏高原深部殼、幔物質(zhì)流展的邊界條件和數(shù)值模型的構建

近四十年來，在青藏高原腹地已取得了大量關于地殼和上地幔分層速度、密度和電性結(jié)構的結(jié)果，盡管在高原的北部、南部、東部和中部有著一定的異同，但是在總體上仍存在著一些基本的共同規(guī)律.在提取青藏高原的殼、幔介質(zhì)平均模型時，我們對人工源地震所得的精細速度結(jié)構結(jié)果給予了更大的權重，而以天然地震和大地電磁測深所得結(jié)果為輔助約束，構建了平均模型.從而得到青藏高原地殼與上地幔分層結(jié)構的基本模型為（滕吉文，2011年在中國科學院地質(zhì)與地球物理研究所學術報告）：a）上地殼，厚約20±5km，P波速度為5.9±0.1km／s；b）地殼低速層，厚約8±2km，P波速度為5.7±0.1km／s；c）中、下地殼厚度為40±10km，P波速度為6.7±0.1km／s；d）上地幔蓋層厚度為30±10km，P波速度為8.1±0.1km／s；e）上地幔軟流層頂部埋藏深度為110±10km，P波速度為7.4±0.1km／s.

基于青藏高原深部的分層地震波速度結(jié)構，現(xiàn)設青藏高原巖石圈由三層介質(zhì)組成，即上地殼、下地殼與巖石圈地幔，巖石圈下方為軟流圈.模型取二維，深度為150km，水平向?qū)挾葹?00km.為了便于討論和明晰起見，分別考慮到上地殼與下地殼之間的低速層以及軟流圈頂面的滑移面，構建了三個模型進行比較：

模型一：速度結(jié)構考慮了兩個低速層.a）上地殼厚度為23km，VP＝5.9±0.1km／s；b）下方23～30km深度處為低速層，VP＝5.7±0.1km／s；c）30～70km處為下地殼，VP＝6.7±0.1km／s；d）巖石圈地幔，即上地幔蓋層位于70～97km深處，VP＝8.1±0.1km／s；e）軟流圈頂部滑移面，位于深部97～103km之間，VP＝7.8±0.1km／s；f）103～150km之間為軟流圈，VP＝7.4±0.1km／s.

模型二：與模型一相比，僅考慮上地殼與下地殼之間的低速層，不考慮軟流圈頂部的滑移面，其分層結(jié)構共計5層.a）上地殼厚度為23km，VP＝5.9±0.1km／s；b）下方23～30km 深度處為低速層，VP＝5.7±0.1km／s；c）30～70km 處為下地殼，VP＝6.7±0.1km／s；d）巖石圈地幔，即上地幔蓋層位于70～100km深處，VP＝8.1±0.1km／s；e）巖石圈下方，即100～150km之間為軟流圈，VP＝7.4±0.1km／s.

模型三：與模型一相比，兩個低速層均不予考慮，只考慮上地殼、下地殼、巖石圈地幔和軟流圈四層.a）上地殼厚度為30km，VP＝5.9±0.1km／s；b）30～70km處為下地殼，VP＝6.7±0.1km／s；c）巖石圈地幔，即上地幔蓋層位于70～100km深處，VP＝8.1±0.1km／s；d）100～150km 之間為軟流圈，VP＝7.4±0.1km／s.

為便于計算，三個模型的邊界均為直立或是水平.以上述幾何模型和地震波速度為基礎，本項工作數(shù)值模型的構建及計算處理均采用有限元軟件Ansys，模型采用了8節(jié)點Plane 183的四邊形平面單元，網(wǎng)格密度在上、下之間的低速層及軟流圈頂部的滑移面均加密為0.5km×0.5km，其余層位均設置為1.0km×1.0km，整個有限元模型單元數(shù)為49229，節(jié)點數(shù)為49629.模型二和模型三在同一地層中的網(wǎng)格密度與之相同，但單元網(wǎng)格數(shù)及節(jié)點數(shù)均有所減少.

在地球動力學數(shù)值模擬研究中，對等效黏滯系數(shù)進行合理的估計，是取得可靠結(jié)果的基礎.Royden等［13］采用牛頓流體的方法研究了高原東部上、下地殼流變性質(zhì)差異和地表變形的關系時，假定地殼淺部（h＜55km）的黏度為1021Pa·s，地殼深部黏度隨深度指數(shù)衰減.Shen等［21］利用冪律流體模擬印度的推擠造成40Ma后青藏高原的抬升高度，認為上地殼黏度在0.5×1021～3×1021Pa·s之間.Flesch等［22］利用偏應力和應變速率資料計算了青藏高原及周邊地區(qū)垂向平均的等效黏滯系數(shù)，結(jié)果認為在青藏高原的等效黏滯系數(shù)在0.5×1022～5×1022Pa·s之間.

參考先前的研究結(jié)果，黏度結(jié)構的選取中，模型選擇青藏高原數(shù)值模擬中常用的“三明治”結(jié)構，即具有力學性質(zhì)強的上地殼和巖石圈地幔，黏度取1.0×1022Pa·s，低速層和滑移面作為軟弱層，黏度取1.0×1020Pa·s.而下地殼和軟流圈的力學性質(zhì)比上地殼和巖石圈地幔弱，比低速層和滑移面強，黏度取1.0×1021Pa·s.模型二和模型三中，相同地層的力學參數(shù)選取與模型一相同.

表1 模型一中的力學參數(shù)Table 1 The mechanical parameters of materials in the 1st model

三個模型在計算過程中均選擇平面應變模式，求解器為軟件自主選擇的稀疏矩陣求解器，該求解器的計算以直接消元法為基礎，病態(tài)矩陣的存在不會引起求解的困難，因而計算過程具有良好的收斂性，但對內(nèi)存及硬盤的配置要求較高.

張培震［23］認為高原東側(cè)相對四川盆地其地殼厚度、地貌高差和垂直形變上均存在橫向壓力梯度，但沒發(fā)現(xiàn)上地殼和巖石圈地幔相對于下地殼的運動，因而通道流模型更可能是Poiseuille模式，而非Couette流.Klemperer［9］認為藏南可能是 Poiseuille流，而藏北可能是Poiseuille流與Couette流的混合模式.為此，本文中的初始模型設置為Poiseuille流的動力模式，即在左邊界施加由高原重力加載隨深度變化的靜巖壓力，在模型的計算中考慮了重力項所帶來的影響，平均靜巖壓力簡化取為3倍靜水壓力值.邊界條件的選取中，模型底部施加豎直向約束，使其水平向能自由移動.模型的右邊界施加水平方向的約束，在豎直方向上自由移動.圖1所示為模型一的卡通圖模型，圖中從上到下六種不同顏色的單元層分別代表著上地殼、低速層、下地殼、巖石圈地幔、軟流圈頂部滑移面和軟流圈六個地層.左邊界的箭頭代表了靜巖壓力的方向，箭頭長短表征了其數(shù)值大小隨深度而變化.

圖1 模型一結(jié)構示意圖Fig.1 Structure sketch of the 1st numerical model

三個模型總的求解時間均為2.5萬年，遠大于巖石圈和軟流圈介質(zhì)的黏彈性松弛時間（約數(shù)十至數(shù)千年），這就使模型最終計算結(jié)果均能進入黏性效應.求解中選擇的時間步為1500步，平均每一步的運算時間步長為16.7a，小于相應黏度條件下的應力松弛時間，以確保計算結(jié)果的可靠.

4 數(shù)值模擬結(jié)果對比

4.1 滑脫面的效應

圖2為三種模型所計算得到的位移場分布圖，圖2d為模型一在0.2Ma后的位移場分布圖.從圖2a、圖2b和圖2c三個模型的對比可以看出，在左邊界動力源附近位移值最大，向右邊界方向逐漸減小.同一水平位置所對應的深度上，模型一的下地殼和巖石圈地幔位移量均比上層的上地殼和下層的軟流圈位移量大，其變化界面正是地殼間低速層和軟流圈頂部的滑脫面.同時，與沒有低速層和滑移面的模型三相比，僅僅2.5萬年之后，模型一的位移場形態(tài)便有了與模型三差異明顯的變化，而且相同位置處位移量比模型三增量超過10%.模型二只具有地殼間的低速層，因而其變形形態(tài)及大小介于模型一和模型三之間.可見，地殼間低速層和軟流圈頂部滑移面的共同存在，對下地殼和巖石圈地幔相對上下層面的快速有序移動起到了明顯的增強作用.

圖2 位移場分布圖（a）模型一；（b）模型二；（c）模型三；（d）0.2Ma后模型一.Fig.2 Distribution of displacement field（a）in model 1；（b）in model 2；（c）in model 3；（d）in model 1after 0.2Ma.

同 時，從圖2d可見，0.2Ma之后，模型一下地殼及巖石圈地幔相對上下層位的位移變化更為明顯，但巖石圈地幔黏度大，力學強度高，抵抗變形的能力比下地殼強，因而位移量弱于下地殼.而上地殼在下地殼快速移動所導致的拖曳作用之下，其位移情況呈現(xiàn)出與巖石圈地幔近乎同步的變形趨勢，這或許也是導致在這種邊界條件下的深部地震波資料解釋中SKS橫波分裂所表征的各向異性現(xiàn)象不明顯的一個因素.

從三個模型里各自取距離左邊界100km處的下地殼和巖石圈地幔中兩個單元節(jié)點，觀測它們的水平位移量隨時間的變化情況，如圖3所示.圖3a中的節(jié)點為下地殼層中距地表深度50km處，圖3b中的節(jié)點則距地表深度為95km處，位于巖石圈地幔地層中部.

圖3a中，萬年時間尺度之后，比較三個模型的下地殼水平位移量，模型三明顯小于模型一和模型二，而模型一和模型二之間下地殼變形的差別一直較小，表明下地殼頂部低速層對下地殼位移響應的影響效果明顯.至于圖3b巖石圈地幔的水平位移量，三個模型的差別顯著地小于它們之間下地殼水平位移的差別量，這或許是與巖石圈地幔強度較大，較難變形相關.但模型一和模型二的巖石圈地幔水平位移差別則較它們之間的下地殼水平位移差別明顯，表明模型二中由于缺乏軟流圈頂部滑移面的作用，導致其巖石圈地幔位移能力較模型一為弱.

圖4為三個模型右邊界的等效應力值與深度變化的對比，圖上等效應力值的轉(zhuǎn)折處清晰地對應各圈層的界面位置.各圈層內(nèi)部的等效應力值比較穩(wěn)定，而圈層間界面處應力變化明顯.可以看到，軟弱層作為滑脫面，它與上、下圈層界面處是應力差最大的區(qū)域，提供了各界面連接處最大的應力差，易在此區(qū)域形成集中且數(shù)值較大的剪切力差，而軟弱層抗剪切能力差，在有破裂或空隙等情況下，容易發(fā)生剪切破壞現(xiàn)象，增加了通道內(nèi)地層物質(zhì)在壓力下有序運移的機率和可能性.下地殼作為力學性質(zhì)軟弱的圈層，在三個模型中均體現(xiàn)為巖石圈中等效應力強度值最小的區(qū)域，而上地殼和巖石圈地幔均是應力強度值較大的區(qū)域，這與青藏高原地震震源深度基本全在0～20km的上地殼和70～110km的巖石圈地幔相符合.

需要注意的是，在遠離力源作用的遠端邊界處，模型二和模型三在Moho面以下的等效應力值趨勢基本一致，而模型一在巖石圈深度以下的軟流圈深度范圍內(nèi)，其等效應力值亦與模型二和模型三基本一致.可見，通道流效應需要有足夠強的初始力源作為保障，而不像物理學意義上的流體一樣，由于抗剪切能力差，較小的力源便能發(fā)生流動.并且，需要考慮其作用范圍，即在初始力源作用下，通道流效應存在著一定的效應范圍，當水平方向位移尺度增加時，通道流效應逐漸減弱.

圖3 所選節(jié)點水平位移與時間關系圖（a）下地殼中；（b）巖石圈地幔中.Fig.3 Relation graph between the horizontal displacement and the time of the selected nodes（a）in the lower crust；（b）in the upper mantle.

圖4 模型右邊界等效應力值Fig.4 Sketch of the equivalent effective stress at the right boundary

4.2 重力及黏度的效應

前人的研究工作也表明，通道流中軟弱層的位移量取決于流體通道的形態(tài)、巖石黏度和部分熔融的程度，以及板塊邊界的位移速率.在通道壁不平行的情況下，如果流體黏度足夠低，差應力可以引起通道內(nèi)高速率的物質(zhì)浮力折返流動［24］.顯見，黏度的選取成為通道流形成的一個重要因素［17］.

Ryder等［25］利用Mw7.6級瑪尼地震后4年變形觀測得到下地殼黏度為4×1018Pa·s.Bürgmann和Dreson［26］的綜述文章里面認為根據(jù)震后蠕變變形而估計的黏滯系數(shù)可以比1019Pa·s更小.應該注意到，大地震震后斷層附近百公里范圍內(nèi)的變形應變速率在數(shù)年內(nèi)可維持在10－13s－1，而陸內(nèi)穩(wěn)定板塊地質(zhì)時間尺度范圍內(nèi)長期變形的應變速率通常小于或等于10－15s－1，兩者應變率相差可達兩個數(shù)量級，同時作用時間也存在數(shù)量級方面的差別，所推算出來的等效黏滯系數(shù)亦可以相差20～30倍.

在上文的計算中，模型黏度參數(shù)在可能的取值范圍內(nèi)取值時相互差別較小.為了探討通道流模型在差別更大的黏度范圍內(nèi)的響應情況，我們采用前文的幾何模型和彈性參數(shù)，重新構建了相對應的模型四、模型五和模型六，在上文所述前人工作的黏度范圍內(nèi)取值.與模型一相對應，模型四的力學參數(shù)如表2所示.

模型依然采用常見的“三明治”結(jié)構，具有力學性質(zhì)較強的上地殼和巖石圈地幔蓋層.軟弱的下地殼黏度比用震后變形所得到的黏度約大20倍左右，低速層和滑移面的黏度比下地殼低一階，取值范圍均在前人研究的范圍內(nèi)［17］.單元的選取和網(wǎng)格的劃分方式以及變形時間及時間步的大小不變，底界面及右邊界的位移約束條件同樣不變.但在左邊界動力學邊界條件的施加中，為了減小各層之間由于密度差所帶來的影響，抵消重力效應，模型左端除了施加靜巖壓力外（數(shù)值上為三倍靜水），深部每一圈層均增加與上方所有圈層由于密度差所帶來的隨深度變化的壓力差值.這樣，壓力從數(shù)值上更大地減小了重力對邊界力抵消所帶來的影響，并確保了模型的力源約束屬于Poiseuille流方式.

表2 模型四的力學參數(shù)Table 2 The mechanical parameters of materials in the 4th model

圖5 位移場計算結(jié)果（a）模型四；（b）模型五；（c）模型六.Fig.5 Displacement field（a）in model 4；（b）in model 5；（c）in model 6.

圖5為第二組模型在2.5萬年之后的位移場計算結(jié)果.圖5a中，由于模型四地殼間低速層和巖石圈底部滑脫面的存在，使得下地殼與上地殼之間以及巖石圈地幔與軟流圈之間顯示位移解耦的現(xiàn)象，下地殼和軟流圈位移明顯大于上地殼和巖石圈地幔，而且越靠近左邊界力源差別越顯著，其差異向右邊界方向逐漸減小.圖5b中，模型五由于巖石圈底部與軟流圈頂部沒有了滑脫面，在模型中部開始，巖石圈地幔與軟流圈之間的位移已逐漸一致.而低速層與巖石圈地幔之間的位移差也逐漸向右遞減，在模型中部下地殼、巖石圈地幔和軟流圈位移趨于一致.而殼內(nèi)低速層的存在也確保了從下地殼開始，深部圈層物質(zhì)的運移比上地殼大.圖5c右端界面處，四個圈層的位移趨于相等.模型六上地殼的位移明顯比模型五小，在其牽拉作用下，下地殼頂部和底部位移差別較大，下地殼頂部Moho面位置成為地殼中位移速度最大的地方，而這難以從介質(zhì)屬性和機理上解釋.

軟弱層的存在，從力學性質(zhì)上降低了研究區(qū)的等效力學強度，增大了內(nèi)部位移及應變變形，有利于軟弱層之間物質(zhì)的有序流動.

5 結(jié)論和討論

（1）在青藏高原南、北向擠壓與東、西方向拉張的構造背景下，青藏高原腹地和四周地殼與上地幔介質(zhì)屬性和深部物質(zhì)的運移軌跡清晰地表明：深部物質(zhì)屬性在縱向和橫向存在明顯不均勻性，構造展現(xiàn)出南北分區(qū)、東西分帶的復雜深、淺格局，并且局部存在有構造變異.在高原深、淺部介質(zhì)與結(jié)構的不均勻性和各向異性作用下，形成了一系列大型走滑斷裂，促使深部殼、幔物質(zhì)以這類走滑斷裂為邊界向東（主體）、西（部分）方向流展，故又導致了淺部（尤其是高原西部）一系列張性斷裂的形成.

青藏高原的通道流模型中［4－7，9，13－18，27］，通常認為通道流的發(fā)生層位，在拉薩塊體南部是中地殼，北部是中、下地殼，而在青藏高原北部的羌塘和松潘—甘孜塊體則發(fā)生在中地殼或者是下地殼.雖然中、下地殼在長時間地質(zhì)年代下呈黏彈性，低Q值，物理學意義上存在流動性，但青藏高原地殼底界的反射波能量強，體現(xiàn)出其Moho面上部的下地殼屬于高地震波速的物質(zhì).這種高地震波速的物質(zhì)，在地震波場的短時間內(nèi)，力學性質(zhì)可體現(xiàn)為彈性.在固定外力的作用下，隨著時間的增加，巖石的應變不斷增加，但應變速率逐漸減慢，最終達到應變以恒定的速率增長，呈現(xiàn)穩(wěn)態(tài)蠕變現(xiàn)象.這時候，在應變率特別低的情況下，應力可以相當?shù)停粫榱?而當深部巖石由于外作用力突然降低為零而保持恒定的形變時，巖石圈內(nèi)的應力會隨時間的增加而逐漸減小，呈現(xiàn)松弛現(xiàn)象.這類呈黏性流動的蠕變和松弛現(xiàn)象，說明深部高溫巖石在一些方面具有物理學上高黏度流體的性質(zhì)，但另一方面，在低于屈服應力的情況下，巖石的形變達到一定程度后將保持不變.因而巖石圈的深部巖石還是呈現(xiàn)為固體.

在下地殼和上地幔蓋層均為地震波速度隨深度增加而遞增的高速固態(tài)物質(zhì)的情況下，即下地殼物質(zhì)和上地幔蓋層之間為固態(tài)物質(zhì)的接觸方式［9，28－33］，這種邊界限定了高原下地殼不可能像物理學意義上一樣任意流動.

同時，高原深部整個地殼屬于同步增厚，Moho界面被諸多深大斷裂所切割，分別構成了喜馬拉雅地體、拉薩地體、羌塘地體、巴顏喀拉地體和昆侖地體，而且Moho面起伏變化劇烈，極為凹凸不平，它不僅不可能構成一個下地殼流動運移的滑移面，而應體現(xiàn)為阻止下地殼物質(zhì)沿其界面運動的阻礙面.另外，若中、下地殼相對于上地殼和上地幔蓋層強行運動，則施加于中地殼或下地殼介質(zhì)上的力必將足夠大，而在其強迫運移的進程中必當會產(chǎn)生一系列新的物理－力學響應現(xiàn)象，如高強度的摩擦生熱效應，沿上、下地殼之間的界面和下地殼與上地幔的交界面處亦會有一系列的同源、同深和一定強度的破裂效應以及一系列小地震的發(fā)生.這種效應會影響這兩個交界面上、下及周邊介質(zhì)發(fā)生變形和地震波場及速度的變化.然而，至今均未見這些必然現(xiàn)象的產(chǎn)生［1，34－35］.同時，通道流模型中，中、下地殼以塑性變形的方式發(fā)生流動，流動速率隨深度增大，在下地殼處達到最大，而Moho界面處流動速率最大或相對很大的分布情況，難以從介質(zhì)屬性和機理上（即高速介質(zhì)＝低黏度介質(zhì)）得到解釋.

由于不具備中、下地殼能夠運移的上、下邊界條件，故中、下地殼難以有序運移，這是任何在給定具有相當誤差范圍黏滯系數(shù)條件下的數(shù)學模擬和推斷所不能解釋的.顯然，青藏高原深部下地殼物質(zhì)流動，必須具備一定的物理－力學邊界條件和可行的運移軌跡，否則難以構成該區(qū)殼、幔物質(zhì)的有序流展.

（2）在青藏高原的地殼中，上地殼底部普遍存在一低速層（低阻層），它與上地殼為軟－硬銜接，故可解耦而構成殼、幔物質(zhì)運移的上滑移面，即提供了一個物質(zhì)可以移動的自由邊界.

下地殼與上地幔蓋層物質(zhì)亦為兩堅硬的固態(tài)物質(zhì)相接，且Moho界面被諸多深大斷裂所切割、形成了一個犬牙交錯的不規(guī)則界面（帶），故不具備可運動的邊界條件，難以發(fā)生相對運動.而青藏高原巖石圈底界只能是上地幔低速層的頂部，在板塊構造中一個十分核心的問題就是巖石圈板塊漂曳在軟流圈上運動，這亦與巖石圈漂曳在其上運動的邊界層相一致，因而下地殼與上地幔蓋層物質(zhì)只能以上地幔軟流層頂部為下滑移面，才有可能在力系作用下促使其運動，下滑移面即為軟流圈的頂部.

青藏高原上、中地殼之間存在一低速層（深度為20±5km），而上地幔軟流層埋深亦相對較淺（深度為110±10km），并向高原邊緣地域減薄、抬升，這兩個滑脫層面是構成高原內(nèi)部殼、幔物質(zhì)能有序運移的唯一邊界條件.這種條件下，不僅下地殼及巖石圈地幔的固態(tài)物質(zhì)較易運移，而且在上、下邊界上不會產(chǎn)生強烈的破裂、變形和熱事件.在研究巖石圈深部流變性質(zhì)時，多種方法的比較研究是不可缺少的.但在比較中要特別注意的是要了解各種方法估算的前提條件，明確各種方法的特點和弱點，以及給出的等效黏滯系數(shù)所隱含的使用條件.因此，這種邊界條件的設定滿足于青藏高原具體地區(qū)的實際觀測結(jié)果和該邊界條件在運動中所導致的效應，同時不必拘泥于在相當誤差條件下對黏度的估計.

而上述的數(shù)值模型，從位移場和等效應力場等的分布和隨時間演化情況等幾個方面的數(shù)值模擬工作證實了軟弱層的存在很大程度地增強了高原內(nèi)部殼、幔物質(zhì)有序運移的動力學機制.

由此，青藏高原的巖石圈結(jié)構，在這樣的殼、幔分層速度結(jié)構為邊界條件下，在兩陸－陸板塊長期碰撞－擠壓作用下，且高原下地殼與上地幔蓋層在其北部和東部均被周邊高速、高電阻率和高密度物質(zhì)所圈閉，阻隔它們“溢出”，故使地殼短縮增厚.待高原在垂直方向增厚達一定高度后，深部物質(zhì)在兩陸－陸板塊強烈擠壓的力系以及高大地形重力勢能差的共同作用下，殼、幔物質(zhì)運移方式在平面為向東、西方向上伸展，即在張應力作用下“突破重圍”，沿大型走滑斷裂系所圍通道進行側(cè)向運移.這種運移方式在深部體現(xiàn)為以上、下地殼間的低速層（低阻層）底界為上滑移面，以上地幔軟流層頂部為下滑移面，巖石圈在地殼低速層與上地殼解耦條件下，并在驅(qū)動力系作用下，實現(xiàn)“下地殼＋上地幔蓋層”高速物質(zhì)在這上、下兩個滑移面之間，沿巖石圈底部，即軟流圈頂部呈同步有序地運動，而不是下地殼流的力學機制［1，34－35］.當然，并非在任何地域的深部都會任意產(chǎn)生這樣的物質(zhì)運移，只能是局部的和在特異環(huán)境下才能實現(xiàn).因為它必須滿足于：1）殼、幔介質(zhì)中存在可供高速物質(zhì)運動的邊界條件，即存在軟弱層，使上、中地殼中的低速層與上地殼解耦，以及上地幔蓋層與軟流圈解耦；2）被移動的介質(zhì)必須有足夠強的力系作用，即南—北向碰撞、擠壓和東—西向拉張；3）存在主體向東的合理的運移軌跡，即高原東部大型走滑斷層為水平面上深部物質(zhì)運移的通道邊界，深部存在“下地殼＋上地幔蓋層”同步運移.

（References）

［1］滕吉文，楊輝，張雪梅.中國地球動力學研究的方向和任務.巖石學報，2010，26（11）：3159－3176.Teng J W，Yang H，Zhang X M.Development direction and task of the geodynamical research in China.Acta Petrologica Sinica（in Chinese），2010，26（11）：3159－3176.

［2］England P C，Houseman G.Finite strain calculations of continental deformation：2.Comparison with the India－Asia collision.J.Geophys.Res.，1986，91（B3）：3664－3676.

［3］Flesch L M，Holt W E，Silver P G，et al.Constraining the extent of crust－mantle coupling in central Asia using GPS，geologic，and shear wave splitting data.Earth Planet.Sci.Lett.，2005，238（1－2）：248－268.

［4］Zhao W L，Morgan W J P.Injection of Indian crust into Tibetan lower crust：a two－dimensional finite element model study.Tectonics，1987，6（4）：489－504.

［5］Bird P.Lateral extrusion of lower crust from under high topography in the isostatic limit.J.Geophys.Res.，1991，96（B6）：10275－10286.

［6］Clark M K，Royden L H.Topographic ooze：building the eastern margin of Tibet by lower crustal flow.Geology，2000，28（8）：703－706.

［7］Beaumont C，Jamieson R A，Nguyen M H，et al.Crustal channel flows：1.Numerical models with applications to the tectonics of the Himalayan－Tibetan orogen.J.Geophys.Res.，2004，109，B06406，doi：10.1029／2003JB002809.

［8］Kuhn T S.The Structure of Scientific Revolutions.Chicago：University of Chicago Press，1979.

［9］Klemperer S L.Crustal flow in Tibet：geophysical evidence for the physical state of Tibetan lithosphere，and inferred patterns of active flow.∥ Law R D，Searle P，Godin L eds.Channel Flow， Ductile Extrusion and Exhumation in Continental Collision Zones.Geological Society，London，Special Publications，2006，268：39－70.

［10］Grujic D，Casey M，Davidson C，et al.Ductile extrusion of the Higher Himalayan Crystalline in Bhutan：evidence from quartz microfabrics.Tectonophysics，1996，260（1－3）：21－43.

［11］Searle M P，Rex A J.Thermal model for the Zanskar Himalaya.J.Metamorphic Geol.，1989，7（1）：127－134.

［12］Nelson K D，Zhao W J，Brown L D，et al.Partially molten middle crust beneath southern Tibet：synthesis of project INDEPTH results.Science，1996，274（5293）：1684－1688.

［13］Royden L H，Burchfiel B C，King R W，et al.Surface deformation and lower crustal flow in eastern Tibet.Science，1997，276（5313）：788－790.

［14］Royden L H.Coupling and decoupling of crust and mantle in convergent orogens：implications for strain partitioning in the crust.J.Geophys.Res.，1996，101（B8）：17679－17705.

［15］Beaumont C，Jamieson R A，Nguyen M H，et al.Himalayan tectonics explained by extrusion of a low－viscosity crustal channel coupled to focused surface denudation.Nature，2001，414（6865）：738－742.

［16］Clark M K，Bush J W M，Royden L H.Dynamic topography produced by lower crustal flow against rheological strength heterogeneities bordering the Tibetan Plateau.Geophys.J.Int.，2005，162（2）：575－590.

［17］王曉芳，何建坤.下地殼管道流動與青藏高原東緣大尺度構造地貌關系.中國科學（D輯），2012，42（4）：505－512.Wang X F，He J K.Channel flow of the lower crust and its relation to large－scale tectonic geomorphology of the eastern Tibetan Plateau.Science in China （Series D），2012，55（8）：1338－1390.

［18］Godin L，Grujic D，Law R D，et al.Channel flow，ductile extrusion and exhumation in continental collision zones：an introduction.∥Law R D，Searle P，Godin L eds.Channel Flow，Ductile Extrusion and Exhumation in Continental Collision Zones.Geological Society，London，Special Publications，2006，268：1－23.

［19］Xypolias P，Kokkalas S.Heterogeneous ductile deformation along a mid－crustal extruding shear zone：an example from the External Hellenides（Greece）.∥Law R D，Searle P，Godin L eds.Channel Flow，Ductile Extrusion and Exhumation in Continental Collision Zones.Geological Society，London，Special Publications，2006，268：497－516.

［20］Hatcher R D Jr，Merschat A J.The Appalachian Inner Piedmont：an exhumed strike－parallel，tectonically forced orogenic channel.∥Law R D，Searle P，Godin L eds.Channel Flow， Ductile Extrusion and Exhumation in Continental Collision Zones.Geological Society，London，Special Publications，2006，268：517－541.

［21］Shen F，Royden L H，Burchfiel B C.Large－scale crustal deformation of the Tibetan Plateau.J.Geophys.Res.，2001，106（B4）：6793－6816.

［22］Flesch L M，Haines A J，Holt W E.Dynamics of the India－Eurasia collision zone.J.Geophys.Res.，2001，106（B8）：16435－16460.

［23］張培震.青藏高原東緣川西地區(qū)的現(xiàn)今構造變形、應變分配與深部動力過程.中國科學（D輯），2008，38（9）：1041－1056.Zhang P Z.The present day′s deformation，strain distribution and deep dynamic process on the Western Sichuan，Eastern Tibetan Plateau.Science in China （Series D）（in Chinese），2008，38（9）：1041－1056.

［24］Gerya T V，St?ckhert B，Perchuk A L.Exhumation of highpressure metamorphic rocks in a subduction channel：a numerical simulation.Tectonics，2002，21（6），1056，doi：10.1029／2002TC001406.

［25］Ryder I，Parsons B，Wright T J，et al.Post－seismic motion following the 1997Manyi （Tibet）earthquake：InSAR observations and modelling.Geophys.J.Int.，2007，169（3）：1009－1027.

［26］Bürgmann R，Dreson G.Rheology of the lower crust and upper mantle：evidence from rock mechanics，geodesy，and field observations.Annu.Rev.Earth Planet.Sci.，2008，36（1）：531－567.

［27］Grujic D.Channel flow and continental collision tectonics：an overview.∥Law R D，Searle P，Godin L eds.Channel Flow，Ductile Extrusion and Exhumation in Continental Collision Zones.Geological Society，London，Special Publications，2006，268：25－37.

［28］Chen W P，Yang Z H.Earthquakes beneath the Himalayas and Tibet：evidence for strong lithospheric mantle.Science，2004，304（5679）：1949－1952.

［29］Fan G W，Lay T.Characteristics of Lg attenuation in the Tibetan Plateau.J.Geophys.Res.，2002，107，B10，2256，doi：10.1029／2001JB000804.

［30］Hacker B R，Gnos E，Ratschbacher L，et al.Hot and dry deep crustal xenoliths from Tibet.Science，2000，287（5462）：2463－2466.

［31］Kusznir N J，Matthews D H.Deep seismic reflections and the deformational mechanics of the continental lithosphere.J.Petrol.，1988，29（1）：63－87.

［32］Ratschbacher L，F(xiàn)risch W，Liu G H，et al.Distributed deformation in southern and western Tibet during and after the India－Asia collision.J.Geophys.Res.，1994，99（B10）：19917－19945.

［33］Wernicke B.The fluid crustal layer and its implications for continental dynamics.∥Salisbury M H，F(xiàn)ountain D M eds.Exposed Cross－Sections of the Continental Crust，NATO ASI Series，1990，317：509－544.

［34］滕吉文，白登海，楊輝等.2008汶川Ms8.0地震發(fā)生的深層過程和動力學響應.地球物理學報，2008，51（5）：1385－1402.Teng J W，Bai D H，Yang H，et al.Deep processes and dynamic responses associated with the Wenchuan Ms8.0 earthquake of 2008.Chinese J.Geophys.（in Chinese），2008，51（5）：1385－1402.

［35］滕吉文.地球內(nèi)部物質(zhì)與能量交換和動力過程.∥21世紀100個交叉科學難題.北京：科學出版社，2005：327－338.Teng J W.The exchange of substance and energy and dynamic process in the earth interior.∥ The 100Crossing－Subject Problems in 21st Century （in Chinese）.Beijing：Science Press，2005：327－338.