底辟熱流體上涌的數(shù)值模擬及其對早古生代青藏高原柴北緣祁連弧形成的啟示

原一哲， 胡才博* ， 魏東平

1 中國科學(xué)院大學(xué)地球與行星科學(xué)學(xué)院， 北京 100049 2 中國科學(xué)院計(jì)算地球動(dòng)力學(xué)重點(diǎn)實(shí)驗(yàn)室， 北京 100049

0 引言

超高壓變質(zhì)巖是研究固體地球不同圈層地球動(dòng)力學(xué)過程的重要依據(jù)，歷來受到地質(zhì)科學(xué)家的高度重視.中央超高壓變質(zhì)帶規(guī)模是世界之最，位于我國中部，東起蘇魯?shù)貐^(qū)，西至阿爾金，全長4000多千米，是研究超高壓變質(zhì)巖形成過程的天然實(shí)驗(yàn)室.超高壓變質(zhì)巖石是如何從地下90 km左右的深度折返至地殼淺層并最終出露地表，在科學(xué)界至今還沒有統(tǒng)一的認(rèn)識(shí).柴達(dá)木北緣超高壓變質(zhì)帶和南祁連造山帶是距今520～400 Ma由古特提斯洋持續(xù)南向俯沖引起的.有學(xué)者認(rèn)為青藏高原北部柴達(dá)木北緣超高壓變質(zhì)帶是由于海洋俯沖過程中底辟流和角流上升過程中的變質(zhì)巖形成的(Lardeaux et al., 2001;Yin et al., 2012).現(xiàn)今地質(zhì)勘察表明，柴達(dá)木北緣超高壓變質(zhì)體在時(shí)間和空間上重疊了早期古生代祁連巖漿弧變質(zhì)帶.

目前，主要有兩種俯沖構(gòu)造模型描述超高壓變質(zhì)巖的形成過程：(1)喜馬拉雅陸陸碰撞模型(Liou et al., 2004)，強(qiáng)調(diào)在陸-陸碰撞中引起的超高壓變質(zhì)作用，超高壓變質(zhì)地體應(yīng)該位于俯沖板塊的被動(dòng)大陸邊緣(O′Brien et al., 2001)，碰撞下行的大陸巖石層地質(zhì)年代學(xué)和地球化學(xué)特征留存于向上運(yùn)移的超高壓變質(zhì)巖石(Jahn, 1999)；(2)安第斯洋陸俯沖模型(Mancktelow, 1995; Dobretsov, 2000)，即由于大洋俯沖和弧形巖漿形成的超高壓變質(zhì)作用.根據(jù)前人的研究，祁連弧超高壓變質(zhì)巖的形成可能由安第斯洋陸俯沖模型來解釋，但巖漿上涌時(shí)空演化的細(xì)節(jié)仍需進(jìn)一步研究.

超高壓變質(zhì)巖石往往隨著海洋板塊俯沖被帶入到大陸90 km以下深處，深俯沖的混雜巖可能轉(zhuǎn)變?yōu)楦呙芏鹊膸r石，然后迅速折返至地殼淺部并部分出露地表.變質(zhì)巖在上升過程中釋放出來的水在浮力的作用下向上進(jìn)入上覆熱的地幔橄欖巖中，導(dǎo)致部分熔融.海洋板塊俯沖導(dǎo)致熱流體上涌存在兩種模型，即管道流模型和底辟流模型，兩者最主要的區(qū)別在于是否形成穩(wěn)定的巖漿通道.圖1a為管道流模型(Lardeaux et al., 2001 )，圖1b為底辟流模型(Yin et al., 2012).本文僅研究底辟流模型，關(guān)于經(jīng)典管道流模型詳見(李忠海等，2015).

圖1 兩種海洋俯沖上涌模型(a) 海洋俯沖管道流模型(修改自Lardeaux et al., 2001)； (b) 海洋俯沖底辟流模型(修改自Yin et al., 2012).Fig.1 Two types of ascent flow during oceanic subduction(a) Channel flow during oceanic subduction (Lardeaux et al., 2001); (b) Diapiric ascent and final emplacement during oceanic subduction (Yin et al., 2012).

很多的研究表明，柴達(dá)木北緣超高壓變質(zhì)帶是早古生代祁連地體的一部分(Zhang et al., 1984; Hsü et al., 1995; Yin and Nie, 1996; Sengoer and Natal′in, 1996; Sobel and Arnaud, 1999; Gehrels et al., 2003a, 2003b).研究柴達(dá)木北緣超高壓變質(zhì)帶可以深入認(rèn)識(shí)古生代祁連地體的演化歷史.總的來說，祁連山—南山地區(qū)包含三個(gè)元古宙被動(dòng)大陸邊緣地質(zhì)構(gòu)造帶，即中央的大雪山構(gòu)造帶，北方的北祁連弧變質(zhì)帶，南方的北柴達(dá)木變質(zhì)帶(Gehrels et al., 2003a, 2003b; Song et al., 2005; Zhang et al., 2005).從空間分布來看，柴達(dá)木北緣早古生代巖漿弧變質(zhì)體和北部柴達(dá)木地區(qū)超高壓變質(zhì)體有很好的一致性(Cowgill et al., 2003; Gehrels et al., 2003a, 2003b)，暗示柴北緣超高壓變質(zhì)體的深部成因.

針對祁連地體發(fā)育過程，Yin等(2012)提出六階段演化觀點(diǎn)： (1) 900～600 Ma，昆侖—柴達(dá)木微板塊邊界被動(dòng)大陸邊緣發(fā)育； (2) 520～400 Ma，祁連海洋板塊的南向俯沖和祁連弧形成，一個(gè)巨大的增生楔沿著昆侖—柴達(dá)木微板塊的北緣不斷向北發(fā)育； (3) 400～375 Ma，祁連弧覆蓋在華北被動(dòng)大陸邊緣，以逆沖方式和華北板塊碰撞；(4)230～100 Ma，從晚三疊紀(jì)到白堊紀(jì)產(chǎn)生弧后擴(kuò)張盆地； (5) 50～40 Ma，印—?dú)W板塊碰撞，祁連山—南山地區(qū)構(gòu)造縮短； (6) 40～0 Ma，華北克拉通的被動(dòng)大陸邊緣局部隆起造成新生祁連弧，與距今520～400 Ma以前早古生代祁連巖漿弧在空間上有重置現(xiàn)象.我們的研究主要關(guān)注祁連地體前三個(gè)階段的地質(zhì)演化過程，如圖2所示.

圖2 早古生代祁連地體和巖漿弧地質(zhì)演化(Yin et al., 2012)(a) 階段一(900～600 Ma)，表示昆侖柴達(dá)木微板塊兩側(cè)被動(dòng)大陸邊緣的發(fā)育； (b) 階段二(520～400 Ma)，表示祁連海洋板塊的南向俯沖和祁連弧的發(fā)育； (c) 階段三(400～375 Ma)，表示祁連弧和華北大陸板塊碰撞.Fig.2 Tectonic evolution of the early Paleozoic Qilian orogen and the magma arc(Yin et al., 2012)(a) Stage 1 (900～600 Ma): Development of passive continental margins on both sides of the Kunlun-Qaidam microcontinent; (b) Stage 2 (520～400 Ma): South-dipping subduction of the Qilian oceanic plate and development of the Qilian arc; (c) Stage 3 (400～375 Ma): Collision between the Qilian arc and North China.

過去幾十年，隨著計(jì)算機(jī)硬件的快速發(fā)展和高性能并行計(jì)算技術(shù)的不斷提升，數(shù)值模擬方法成為理論解析解和物理模擬之外的地球動(dòng)力學(xué)的重要方法，能夠研究在復(fù)雜的幾何、材料、邊界條件之下地球介質(zhì)在相對真實(shí)時(shí)空尺度的變形和流動(dòng).數(shù)值模擬研究主要包括兩大類：簡化模型和綜合模型.簡化模型只考慮簡化的邊界條件和幾何模型，忽略板塊運(yùn)動(dòng)中的很多關(guān)鍵物理參數(shù)，如不考慮相變問題，僅設(shè)定最簡單的溫度分布和邊界速率等；綜合模型能夠更加真實(shí)地反映地球動(dòng)力學(xué)過程，加入了很多物理、化學(xué)參數(shù)，比如復(fù)雜的流變學(xué)性質(zhì)，溫度和壓力變化導(dǎo)致的巖石相變過程，可以模擬實(shí)時(shí)變化的溫度、密度、黏度場等.底辟流上涌數(shù)值模擬難點(diǎn)主要體現(xiàn)在以下幾個(gè)方面： (1) 要合理地模擬再現(xiàn)底辟上涌過程的各個(gè)階段； (2) 各圈層之間巨大的黏滯度差異； (3) 自由表面的邊界條件； (4) 復(fù)雜的巖體黏-彈-塑性本構(gòu)關(guān)系(Hall et al., 2003)； (5) 復(fù)雜的相變過程(Hebert et al., 2009)； (6) 地幔內(nèi)流體/熔體的遷移(Davies and Stevenson, 1992).

目前較為流行研究底辟流上涌過程的數(shù)值模擬算法有以下幾種： (1) FEM-BEM：耦合了固體和流體兩種算法(Morra et al., 2006)，地幔黏性流體區(qū)域使用了邊界元算法(BEM)中的格林函數(shù)；將巖石圈和平板變形視為固體材料，使用有限元法(FEM)進(jìn)行計(jì)算.(2)Underworld：一種流體動(dòng)力學(xué)方法，使用標(biāo)準(zhǔn)Euler有限元網(wǎng)格離散化模型，并采用拉格朗日粒子表示不同的材料(Moresi et al., 2007).(3)E-Script：運(yùn)用“水平集”的數(shù)值方法(Gross et al., 2007)，該方法常常用于追蹤物體界面及形狀，離散模型則使用有限單元網(wǎng)格.(4)FDCON：一種有限差分算法.將動(dòng)量方程和質(zhì)量守恒方程改寫為與黏滯度和流線方程相關(guān)的雙調(diào)和方程(Schmeling and Marquart, 1991)，用Kolesky分解法直接求解.(5)LAPEX-2D：該算法使用網(wǎng)格-粒子方法(Sulsky et al., 1995)和顯式拉格朗日有限差分技術(shù)(Babeyko et al., 2002)求解質(zhì)量、動(dòng)量和能量平衡方程，適用于求解黏彈性問題.(6)CITCOM：一種有限元算法(Moresi and Solomatov, 1995).有限單元是以恒定壓力和線性速度為特征的雙線型四邊形.每個(gè)單元都進(jìn)行屬性內(nèi)插并且直接應(yīng)用于積分點(diǎn).(7)I2VIS：基于有限差分方法和粒子法的混合方法(Gerya and Yuen, 2003a, 2007)，目前應(yīng)用廣泛.

本文所采用的I2VIS數(shù)值模擬是基于單元粒子和有限差分混合的計(jì)算方法(Gerya and Yuen, 2003a, 2007).模型參數(shù)的敏感性實(shí)驗(yàn)通過多因素多水平正交試驗(yàn)設(shè)計(jì)方法實(shí)現(xiàn)，“因素”的多少根據(jù)考察參數(shù)的個(gè)數(shù)而定，而“水平”的多少則主要依據(jù)實(shí)際地質(zhì)過程中的取值范圍而定.在模型實(shí)驗(yàn)結(jié)果的基礎(chǔ)上，為了得出更為直觀而精確的定量化結(jié)果采用數(shù)理方法進(jìn)行量化分析.

柴達(dá)木北緣弧巖漿變質(zhì)帶和超高壓變質(zhì)地體時(shí)空上具有良好一致性，因此需要考慮不同圈層的耦合作用.本文分別采用簡單模擬和綜合模擬的方法，根據(jù)有無溫度耦合作用和是否存在上涌通道，分別建立三組數(shù)值模型： (1) 未考慮溫度影響的底辟流模型； (2) 考慮溫度耦合的無上涌通道底辟流模型； (3) 考慮溫度耦合且具有上涌通道的底辟流模型.通過有限差分法和粒子追蹤技術(shù)，研究巖石圈、軟流圈尺度內(nèi)地球流體介質(zhì)在熱、重力及構(gòu)造作用下的流動(dòng)、變形及應(yīng)力狀態(tài)，定量模擬底辟流上涌動(dòng)力學(xué)過程.通過這三個(gè)模型我們可以模擬出底辟流上涌過程中速度、溫度以及黏度的時(shí)空變化過程，并與流體冪次定律速度公式(Weinberg and Podladchikov, 1994)以及Hall和Kincaid(2001)提供的針對海洋俯沖板塊流體上涌機(jī)理進(jìn)行對比研究，對柴達(dá)木北緣古生代祁連弧的形成過程具有重要的啟示意義.

1 數(shù)值模擬方法

1.1 控制方程

本文模型采用二維的地球動(dòng)力學(xué)算法“I2VIS”進(jìn)行數(shù)值計(jì)算(Gerya and Yuen, 2003a).該算法程序考慮了熱-流體動(dòng)力學(xué)的耦合，有三組控制方程：斯托克斯流體動(dòng)力學(xué)方程、不可壓縮質(zhì)量守恒方程、溫度演化的能量守恒方程.

(1)斯托克斯流體動(dòng)力學(xué)方程：

(1)

(2)

式中，g表示重力加速度，密度ρ是巖石類型C、部分熔融比例M、壓力P及溫度T的函數(shù)，σ′xx、σ′zz和σ′xz是偏應(yīng)力張量分量.

(2)不可壓縮質(zhì)量守恒方程：

(3)

其中，vx和vz表示水平速度和垂向速度.

(3)溫度演化的能量守恒方程：

(4)

(5)

(6)

1.2 巖體流變性質(zhì)

本文采用黏-塑性本構(gòu)方程表示偏應(yīng)力和應(yīng)變率張量，對不可壓縮的黏滯性變形則采用下式表示：

(7)

(8)

(9)

式中，ηeff表示有效黏滯系數(shù)，將熱-流體動(dòng)力學(xué)耦合在一起，可表示為

(10)

實(shí)際巖體的黏-塑性流變結(jié)合了韌性與塑/脆性流變.因此，我們在數(shù)值計(jì)算中采用改進(jìn)的Drucker-Prager屈服準(zhǔn)則(Ranalli, 1995)：

σyield=C0+Psin(φeff),

(11)

sin(φeff)=sin(φ)(1-λ),

(12)

(13)

地球動(dòng)力學(xué)算法“I2VIS”(Gerya and Yuen, 2003a)，選取ηductile和ηplastic最小者作為實(shí)際的黏度參與計(jì)算(Ranalli, 1995):

ηcreep=min(ηductile,ηplastic).

(14)

1.3 部分熔融

模型計(jì)算中考慮了地殼巖石的部分熔融行為(Gerya and Yuen, 2003b; Burg and Gerya, 2005)，根據(jù)實(shí)驗(yàn)巖石學(xué)的約束條件，將溫度與巖石部分熔融體積比例之間的關(guān)系看做近似的線性關(guān)系：

M=0，當(dāng)T≤Tsolidus,

(15a)

(15b)

M=1，當(dāng)T≥Tliquidus,

(15c)

式中,Tsolidus和Tliquidus分別表示巖石的濕固相線溫度和干液相線溫度.部分熔融巖石的有效密度與熔融比例之間的關(guān)系如下所示：

ρeff=ρsolid-M(ρsolid-ρmolten),

(16)

式中，ρsolid和ρmolten分別表示固相巖石的密度和熔融巖石的密度，兩者均為溫度和壓力的函數(shù)：

ρ=ρ0(1-α(T-T0))(1+β(P-P0)),

(17)

式中，ρ0代表巖石在P0=0.1 MPa壓強(qiáng)條件，T0=298 K溫度條件下測得的標(biāo)準(zhǔn)密度；α和β分別代表熱膨脹系數(shù)和可壓縮系數(shù).

以上公式表明，地球介質(zhì)的黏度、密度變化受溫度和壓力的影響非常大，而隨溫度、壓力變化的黏度、密度對固體地球不同圈層的俯沖、上涌等動(dòng)力學(xué)過程也有直接的影響，因此，以上公式構(gòu)成了地球動(dòng)力學(xué)熱-流體耦合數(shù)值模擬的理論基礎(chǔ).

2 柴北緣祁連弧底辟流上涌過程的數(shù)值模擬

2.1 未考慮溫度影響的底辟流上涌模型

未考慮溫度影響的底辟流上涌模型主要研究底辟流的初始半徑對底辟流上涌速度的影響，為論述的方便，下文稱為模型一.考慮到柴達(dá)木北緣超高壓變質(zhì)地體在綠良山、錫鐵山、都蘭區(qū)域南北方向?qū)?5～50 km(Zhang et al., 2005)，并且熱流體或熔體上涌時(shí)伴隨著自身物質(zhì)的擴(kuò)散和體積不斷增大，最終在地殼中的侵蝕范圍大于上涌初始半徑，本模型設(shè)置三組不同半徑(5 km, 10 km, 15 km)的底辟流上涌模型進(jìn)行對比實(shí)驗(yàn).模型一未考慮溫度對底辟上涌的影響，黏度和密度均采用為常數(shù).初始模型寬度50 km，深度120 km.幾何模型分為三個(gè)不同的區(qū)域(如圖3)，不同區(qū)域材料賦值不同：1是地殼模塊，厚度35 km，黏度1024Pa·s，密度2700 kg·m-3；2是地幔模塊，厚度95 km，黏度1020Pa·s，密度3300 kg·m-3；3是初始底辟流，初始黏度1012Pa·s，密度2500 kg·m-3.底辟流初始中心位置位于地表下90 km處開始上涌，模型離散化采用規(guī)則的有限差分網(wǎng)格，網(wǎng)格劃分為兩個(gè)部分：第一，網(wǎng)格的劃分分辨率50×120，相當(dāng)于1 km一個(gè)格網(wǎng)點(diǎn)；第二，網(wǎng)格內(nèi)部標(biāo)記點(diǎn)的劃分，這里采用10×10的標(biāo)記，相當(dāng)于每100 m一個(gè)標(biāo)記點(diǎn)粒子.模型邊界均采用自由滑動(dòng)邊界條件.

圖3 半徑為10 km底辟流初始密度示意圖上圖中，1表示地殼，2表示上地幔，3表示初始底辟流.Fig.3 Schematic diagram of initial density of diapir flow with a radius of 10 kmIn the above figure, 1 denotes the earthcrust, 2 denotes the upper mantle, 3 denotes initial diapir flow.

底辟流上涌過程主要分為四個(gè)階段，這里以半徑為10 km的底辟流上涌為例說明(如圖4)：(1)第一階段是快速上涌階段：0～0.01 Ma，底辟流上涌過程中由于密度差形成的正浮力作用值較大，因此在底辟上涌的初期，初始速度值在短時(shí)間內(nèi)從0迅速增大至一個(gè)峰值，由于加速過程非常短暫，這段時(shí)間相對于模擬研究的時(shí)間尺度而言可以忽略不計(jì)，這里為了研究方便統(tǒng)一取上涌持續(xù)0.01 Ma后的速度為起點(diǎn)進(jìn)行分析.(2)第二階段是勻減速上涌階段：0.01～0.5 Ma，底辟流上涌速度均勻減小，流體和周圍地幔進(jìn)行交代作用，底辟流初始上涌速度為2.1 cm·a-1，上涌速度均勻下降至0.6 cm·a-1以后進(jìn)入非均勻減速階段；(3)第三階段是非均勻減速上涌階段：0.5～1.8 Ma，底辟流上涌速度隨著時(shí)間的推移非均勻放緩，在這個(gè)時(shí)間段內(nèi)上涌速度逐漸減緩為零，同時(shí)由于底辟流屬于高塑性流體，流體組分逐漸向水平方向擴(kuò)散；(4)第四個(gè)階段是穩(wěn)定階段：1.8～3.2 Ma，該階段底辟流上涌速度幾乎為零，底辟流雙向擴(kuò)散均達(dá)到穩(wěn)定狀態(tài).

圖4 半徑為10 km的底辟流上涌速度圖(從0.01 Ma作為起點(diǎn)，不包括階段1)Fig.4 Upwelling velocity of diapir with a radius of 10 km (Starting from 0.01 Ma, excluding stage 1)

以底辟流中心位置為觀測點(diǎn)，得出半徑為10 km底辟流上涌過程中上涌速度變化的實(shí)時(shí)曲線，如圖4所示.

圖5為半徑為r=10 km的底辟流半徑上升速率模擬示意圖.

圖5 半徑為r=10 km的底辟上涌模擬示意圖(水平速度Vx，垂直速度Vz，密度及黏度)(a) 0.01 Ma； (b) 0.7 Ma； (c) 1.8 Ma.Fig.5 Simulation diagram of diapir upwelling with radius r=10 km (Horizontal speed Vx, vertical speed Vz, density and viscosity)

半徑為5 km, 10 km, 15 km的底辟流上涌速度如圖6所示.

通過對比分析三個(gè)不同半徑(r=5 km,10 km,15 km)的底辟流垂直上涌速度(如圖6)，我們得到以下實(shí)驗(yàn)結(jié)論：

圖6 三個(gè)不同半徑(r=5 km,10 km,15 km)的底辟流上升速度圖Fig.6 Ascent speed graph of diapir flow with three different radius (r=5 km,10 km,15 km)

(1)半徑大小對底辟流上涌速度的變化具有明顯的影響，底辟流半徑越大所受到的正浮力作用越大，上涌速度越快，達(dá)到穩(wěn)定階段的時(shí)間也越短，半徑為5 km的底辟流在開始上涌2.2 Ma以后速度為零，達(dá)到穩(wěn)定狀態(tài)；而半徑為10 km、15 km的底辟流分別在1.8 Ma、1.5 Ma以后達(dá)到穩(wěn)定狀態(tài).

(2)底辟流勻減速上涌階段，半徑越大單位時(shí)間內(nèi)降低的速度值越少，半徑5 km的底辟流每百萬年下降速率為4.3 cm·a-1，而半徑10 km、15 km的底辟流每百萬年下降速率分別為2.6 cm·a-1、1.9 cm·a-1.

(3)在不考慮熔融和溫度的影響下，底辟流上涌過程非常迅速，半徑大小對上涌速度的影響至關(guān)重要，三個(gè)不同半徑(5 km、10 km、15 km)的底辟流從開始上涌至速度為零的穩(wěn)定時(shí)間均不超過3 Ma.

綜上所述，通過實(shí)驗(yàn)?zāi)M三組不同半徑的底辟流上涌過程得到相關(guān)的實(shí)驗(yàn)結(jié)果為：半徑5 km、10 km、15 km的底辟流平均上涌速度分別為2.7 cm·a-1、3.3 cm·a-1、4.0 cm·a-1.由上圖可知，在同等地質(zhì)條件下，隨著半徑的增大底辟流上升至下地殼的速度越快，所消耗的時(shí)間越短.三個(gè)不同半徑的底辟上涌穩(wěn)定時(shí)間都不超過3 Ma，這和目前的觀測成果底辟流從90 km深的地幔深度運(yùn)移至地殼30 km處耗時(shí)小于10 Ma在時(shí)間上具有較好的一致性(Yin et al., 2012).

2.2 考慮溫度耦合的無上涌通道底辟流上涌模型

在熱的底辟流上涌過程中，不可避免地伴隨著與周圍地幔的熱交換作用，自身溫度黏度和密度都會(huì)實(shí)時(shí)發(fā)生變化.而模型一未考慮溫度變化的影響，黏度、密度均為常數(shù).為了更加真實(shí)的模擬現(xiàn)實(shí)的底辟流上涌過程，我們考慮巖石的熔融效應(yīng)，在流體速度場的基礎(chǔ)上耦合溫度場，建立溫度耦合的無上涌通道底辟流上涌模型，為論述的方便，將之稱為模型二.與模型一進(jìn)行對比分析，研究溫度變化對底辟流上涌過程的影響.同時(shí)為了對比分析多個(gè)底辟流上涌的影響，將模型二分為單底辟流上涌模型和雙底辟流上涌模型.

如前文所述，底辟流上涌過程中溫度的變化會(huì)導(dǎo)致自身黏度以及密度隨時(shí)間變化.模型二中，我們假設(shè)巖漿來自巖石圈下面的巖漿源區(qū)(SMSR)(Gerya and Burg, 2007)，初始的熔融體為20%的熔體分?jǐn)?shù)，隨著底辟流的上涌，熔融的巖漿溫度逐漸下降，熔體分?jǐn)?shù)也隨之改變，巖漿上涌至地表時(shí)與巖石圈發(fā)生熔融反應(yīng)并逐漸冷卻.模型二的主要參數(shù)見表1.

表1 主要模型參數(shù)(參考Ranalli, 1995)Table 1 Main model parameters (Refer to Ranalli, 1995)

模型二寬度100 km，深度120 km，底辟流初始中心位于地表下90 km處(圖7)，模型頂部設(shè)置5 km沉積層，沉積層下地殼厚度為35 km，其中上地殼20 km，下地殼15 km.上地殼、下地殼和底辟流分別近似采用石英巖，斜長石和基性熔融體的流變參數(shù)，見表1.關(guān)于熱邊界條件，模型頂部保持溫度0 ℃，大陸巖石圈底部邊界溫度1300 ℃(Turcotte and Schubert, 2002)，軟流圈地幔的溫度梯度為 0.5 ℃/km，兩側(cè)邊界的水平方向溫度梯度為零(即零熱流).模型底部采用的是外部邊界固定溫度條件(Li et al., 2010)，即在模型底部120 km處的滲透邊界上，溫度和熱流均可以隨著模型計(jì)算而動(dòng)態(tài)調(diào)整.模型二所采用的不同巖層初始黏度和密度值參考模型一，模型離散化采用規(guī)則的有限差分網(wǎng)格，網(wǎng)格劃分為兩個(gè)部分：第一，格網(wǎng)點(diǎn)的劃分分辨率100×120，相當(dāng)于1 km一個(gè)格網(wǎng)點(diǎn)；第二，網(wǎng)格內(nèi)部標(biāo)記點(diǎn)的劃分，這里采用10×10的標(biāo)記，相當(dāng)于每100 m一個(gè)標(biāo)記點(diǎn)粒子.模型頂部設(shè)置黏滯層，地殼表面之上，與自由滑動(dòng)的模型頂界面之間，設(shè)計(jì)一層相對高黏滯度的偽空氣層，其與上地殼的接觸面用以模擬地形起伏面，該地形起伏面包含了近似的剝蝕和沉積作用(如Li and Gerya, 2009).模型中除了底部邊界，其他邊界均采用自由滑動(dòng)速度邊界條件.

圖7 單個(gè)底辟流上涌示意圖(a) 初始溫度示意圖； (b) 初始黏度示意圖； (c) 初始密度示意圖.Fig.7 Schematic diagram of a single diapir flow upwelling(a) Initial temperature; (b) Initial viscosity; (c) Initial density.

單個(gè)底辟流上涌實(shí)驗(yàn)結(jié)果表明，底辟流上涌過程中由于自身溫度較高和周圍地幔的溫差大，上涌初期進(jìn)行劇烈熱交換，隨著上涌時(shí)間的延續(xù)溫度迅速減小，并且伴隨著自身物質(zhì)擴(kuò)散和體積不斷增大，最終形成以底辟流中心為分界點(diǎn)的雙拱形地幔對流，在上涌2.5 Ma以后，在地表下75 km左右停止上涌(圖8).這是由于底辟流在上涌過程中，不僅僅是流體本身的運(yùn)動(dòng)，還存在由于溫度差異而引起的密度、黏度等物理參數(shù)的實(shí)時(shí)變化.伴隨著上涌過程中不斷地與周圍環(huán)境進(jìn)行熱交換導(dǎo)致自身熱量的流失，底辟流的黏度和密度隨之增大，導(dǎo)致底辟流上涌至地表下75 km左右時(shí)溫度和黏度與周圍地幔相近，正浮力為零，因此相比常黏度模型相同半徑的底辟流上涌時(shí)不能完全抵達(dá)殼幔交界處.溫度場是底辟流上涌過程重要影響因素.單個(gè)底辟流無法上涌至殼幔邊界處，難以穿越地幔楔直接出露地表.

圖8 半徑為10 km的單個(gè)底辟上涌模型圖 (1.0,2.0,3.0 Ma三個(gè)不同時(shí)刻溫度、密度和黏度的分布)Fig.8 Single dip upwelling model with a radius of 10 km (The distribution of temperature, density and viscosity at three different moments of 1.0, 2.0, 3.0 Ma)

根據(jù)柴達(dá)木北緣祁連弧地區(qū)地表超高壓變質(zhì)巖的分布情況可知底辟流上涌過程并非是單個(gè)上涌.伴隨著海洋板塊的俯沖作用，底辟流上涌實(shí)際上是一個(gè)持續(xù)增量的過程，為了進(jìn)一步探究底辟流上涌的情形，我們設(shè)計(jì)雙底辟流上涌模型(如圖9所示)，分別在地表下90 km，70 km設(shè)置半徑10 km的底辟流，模型參數(shù)與單個(gè)底辟流一致，模型運(yùn)行結(jié)果如圖10所示.

圖9 雙底辟流上涌示意圖(a) 初始溫度示意圖； (b) 初始黏度示意圖； (c) 初始密度示意圖.Fig.9 Schematic diagram of double diapir flow upwelling(a) Initial temperature; (b) Initial viscosity; (c) Initial density.

圖10 半徑為10 km的雙底辟上涌模型圖 (1.0,2.0,3.0 Ma三個(gè)不同時(shí)刻溫度、密度和黏度的分布)Fig.10 Double diapir upwelling model with a radius of 10 km (The distribution of temperature, density and viscosity at three different moments of 1.0, 2.0, 3.0 Ma)

雙底辟流上涌實(shí)驗(yàn)結(jié)果表明，底辟流上涌過程中伴隨著與周圍巖體劇烈的熱交換和自身物質(zhì)擴(kuò)散，從最初上涌時(shí)的兩個(gè)半徑為10 km的球形迅速合成一股較大的底辟流，同時(shí)以新形成的大的底辟流中心為分界點(diǎn)雙向侵蝕，在上涌至3 Ma以后，持續(xù)上涌至地表下60 km左右(圖10)，而單個(gè)底辟流則上涌至2.5 Ma時(shí)則停止上涌.在相同的時(shí)間內(nèi)(3 Ma)，雙底辟流上涌比單底辟流在垂直方向侵蝕的距離更長，并且水平方向也沒有發(fā)生彌散現(xiàn)象.因此，我們認(rèn)為當(dāng)?shù)妆倭魃嫌康牧孔銐虼髸r(shí)，會(huì)在原來底辟流上涌的位置上方形成穩(wěn)定的巖漿通道，在此基礎(chǔ)上，底辟流會(huì)在正浮力的推動(dòng)下沿著巖漿通道上涌，不斷地侵蝕殼幔邊界，在短時(shí)間內(nèi)(小于1 Ma)上涌融體不斷地侵蝕減薄上覆地殼，最終出露地表形成巖漿弧.當(dāng)巖漿通道足夠?qū)挄r(shí)，隨著俯沖的進(jìn)行，在熔融過程中密度和黏度相對較小的變質(zhì)巖沿著巖漿通道快速上涌至地表，形成地表可觀測到的超高壓變質(zhì)巖石.

由于巖漿通道形成后的上涌時(shí)間很短，北柴達(dá)木超高壓變質(zhì)帶和巖漿弧的形成在時(shí)間和空間上一致.

2.3 考慮溫度耦合和上涌通道的底辟流上涌模型

為了進(jìn)一步說明柴達(dá)木北緣祁連弧的形成機(jī)理，針對海洋俯沖底辟流模式中存在穩(wěn)定巖漿上涌通道(Yin et al., 2012)，我們建立溫度耦合和上涌通道的底辟流上涌模型，為論述的方便，下文稱為模型三(如圖11).模型三設(shè)置了寬為5 km的軟弱層模擬巖漿上涌通道.模型三所采用的不同巖層初始黏度和密度值參考模型一，模型離散化采用規(guī)則的有限差分網(wǎng)格，模型尺寸為100 km×120 km，網(wǎng)格分辨率為100×120，相當(dāng)于1 km一個(gè)網(wǎng)格；第二，網(wǎng)格內(nèi)部標(biāo)記點(diǎn)的劃分，這里采用10×10的標(biāo)記，相當(dāng)于每100 m一個(gè)標(biāo)記點(diǎn)粒子.地層結(jié)構(gòu)、熱邊界條件參考模型二.上地殼、下地殼和底辟流采用的流變參數(shù)見表1.模型頂部設(shè)置黏滯層，即地殼表面之上，與自由滑動(dòng)的模型頂界面之間，設(shè)計(jì)一層相對高黏滯度的偽空氣層，其與上地殼的接觸面用以模擬地形起伏面，該地形起伏面包含了近似的剝蝕和沉積作用(Li and Gerya, 2009).模型中除了底部邊界，其他邊界均采用自由滑動(dòng)速度邊界條件.

圖11 模型三示意圖(a) 初始溫度示意圖； (b) 初始黏度示意圖； (c) 初始密度示意圖.Fig.11 Schematic diagram of model 3(a) Initial temperature; (b) Initial viscosity; (c) Initial density.

圖11給出了模型三0.01, 0.05, 0.1 Ma三個(gè)不同時(shí)刻溫度、密度和黏度的分布.模型三的結(jié)果表明，模型運(yùn)行初期0.01 Ma內(nèi)，底辟流快速上涌冷卻，同時(shí)伴隨著周圍巖體快速熱傳導(dǎo)增溫，尤其以初始位置正上方的上地殼巖體升溫最為強(qiáng)烈.0.01～0.05 Ma內(nèi)，下地殼內(nèi)基性巖漿在巨大正浮力作用下緩慢而穩(wěn)定地上涌侵入上地殼底部并水平向外擴(kuò)展.0.05 Ma后，底侵巖漿初具規(guī)模，下地殼基性巖漿開始快速侵入，此后經(jīng)過約0.05 Ma，幾乎所有的巖漿侵入到上地殼底部，平鋪于下地殼界面之內(nèi).在快速上升熔體的沖擊及熔體和軟流圈物質(zhì)對巖石圈地幔底部熱侵蝕的共同作用下，對上覆巖石圈地幔底部造成一定程度的破壞.然后持續(xù)增多的板片熔體對上覆巖石圈地幔的破壞逐漸向巖石圈淺部蔓延，加深減薄程度；同時(shí)在熔體上升過程中誘發(fā)的地幔對流作用下，上覆巖石圈地幔的破壞也逐漸向陸內(nèi)方向擴(kuò)展，擴(kuò)展減薄范圍，從而導(dǎo)致對上覆大陸在垂向和水平向的雙重破壞，如圖10所示，半徑為10 km的底辟流上涌過程中，垂向地殼減薄約15 km，水平向陸內(nèi)侵蝕延伸約30 km，隨著底辟流持續(xù)增量上涌，在地殼內(nèi)部形成蘑菇狀侵入體.即具有穩(wěn)定巖漿通道的底辟流上涌過程中，從地幔深處90 km上涌至殼幔邊界時(shí)在水平方向的侵蝕范圍要比垂向的侵蝕范圍大一倍左右.

另外，模型黏度場和溫度場結(jié)果顯示，底辟流上涌誘發(fā)的向上流場強(qiáng)烈侵蝕作用導(dǎo)致上覆巖層溫度升高300～500 ℃，黏度降低為初始時(shí)刻的1/100～1/1000.隨著底辟流上涌在巖石圈地幔向淺部和向陸的雙重破壞，具有更低黏度，更高熱流的軟流圈物質(zhì)也會(huì)沿著已經(jīng)形成的巖漿通道上涌，同時(shí)雙向蔓延導(dǎo)致地幔相對較淺部位被弱化，產(chǎn)生劇烈熔融作用，使得上涌熔體在莫霍面附近的大陸巖石圈地幔底部逐漸冷卻并穩(wěn)定下來.實(shí)際情況不僅僅是單個(gè)底辟流上涌，而是伴隨含水量較高的海洋板塊俯沖，在地幔深處高溫高壓條件下板塊迅速發(fā)生脫水熔融，上涌熔體的體積不斷增大，底辟流體入侵地殼雙向侵蝕的范圍也隨之不斷擴(kuò)大，最終出露地表形成巖漿弧，該過程對古生代柴達(dá)木北緣地區(qū)熱地幔物質(zhì)的熔融溢出形成祁連巖漿弧提供理論參考意義.

3 討論

3.1 三個(gè)模型對比分析及意義

底辟流上涌實(shí)質(zhì)上對應(yīng)地幔對流過程中的上升流，不僅為地球板塊運(yùn)動(dòng)和洋脊擴(kuò)張?zhí)峁┝酥匾膭?dòng)力來源，同時(shí)為地球內(nèi)部各圈層物質(zhì)遷移提供了重要輸出.本文實(shí)驗(yàn)?zāi)Ｐ途哂幸欢ǖ钠者m性，討論了底辟流上涌的一般性問題，研究成果可以為古生代柴達(dá)木北緣地區(qū)熱地幔物質(zhì)的熔融溢出形成祁連巖漿弧的過程提供理論參考意義.三個(gè)模型的實(shí)驗(yàn)結(jié)論及優(yōu)缺點(diǎn)如下：

圖12 不同時(shí)刻半徑為10 km的底辟流上涌模型圖(a) 0.01 Ma； (b) 0.05 Ma； (c) 0.1 Ma.Fig.12 The upwelling model of the diapir with a radius of 10 km at different time

圖13 柴北緣都蘭南部區(qū)域單個(gè)底辟流上涌初始模型(a) 初始溫度示意圖； (b) 初始黏度示意圖； (c) 初始密度示意圖.Fig.13 Initial model of upwelling of a single diapir in South-Dulan region of North Qaidam(a) Initial temperature; (b) Initial viscosity; (c) Iinitial density.

圖14 不同時(shí)刻柴北緣都蘭南部區(qū)域單個(gè)底辟流上涌模型圖(a) 1.0 Ma； (b) 2.0 Ma； (c) 3.0 Ma.Fig.14 Diapir flow upwelling model diagram in South-Dulan region of North Qaidam at different times

模型一主要是從運(yùn)動(dòng)學(xué)的角度分析了底辟流上涌的速度變化，缺點(diǎn)是沒有考慮溫度場對上涌過程中底辟流的密度和黏度變化，結(jié)果表明在不考慮溫度場的情況下，底辟流上涌過程中形狀基本不變，單個(gè)底辟流上涌可以到達(dá)下地殼，但無法到達(dá)地表；為了更加真實(shí)地模擬底辟流上涌過程，模型二在模型一的基礎(chǔ)上增加了溫度場，模型結(jié)果表明單個(gè)底辟流無法上涌至地表，并且底辟流上涌過程中的形狀不能保持為原來的球形，而是逐漸擴(kuò)散，未到達(dá)下地殼便彌散和消亡在地幔中；而兩個(gè)底辟流上涌時(shí)，底辟流雙向侵蝕距離明顯加大，并且在相同的運(yùn)行時(shí)間內(nèi)(3 Ma)底辟流保持球形形狀，并未發(fā)生彌散現(xiàn)象，而且侵蝕作用也并未停止.由此我們推出，持續(xù)增量上涌的底辟流必然會(huì)形成直達(dá)地表的穩(wěn)定巖漿通道，巖漿通道的形成導(dǎo)致不斷補(bǔ)充的底辟流在水平和垂直兩個(gè)方向持續(xù)侵蝕周圍巖體，最終形成出露地表的巖漿弧.模型三在模型二的基礎(chǔ)上模擬了形成穩(wěn)定巖漿通道后底辟流上涌的情形，結(jié)果表明，底辟流上涌時(shí)間大大縮短，侵蝕效率也明顯加快，底辟流在0.1 Ma 左右即從地表下90 km 上涌至地表，最終形成出露地表的巖漿弧.正是由于巖漿通道形成后極大的縮短了上涌時(shí)間，導(dǎo)致北柴達(dá)木超高壓變質(zhì)帶和巖漿弧在形成時(shí)間和空間上一致.

綜上所述，模型一沒有考慮溫度的影響，只考慮底辟流半徑對底辟流上涌速度的影響；模型二考慮了溫度耦合的影響，但沒有考慮上涌通道的影響，實(shí)際上是海洋板塊俯沖時(shí)，底辟流上涌初期階段過程，其最終結(jié)果是形成上涌通道；模型三是考慮溫度耦合和上涌通道的底辟流上涌模型，伴隨著持續(xù)增量的底辟流上涌，在底辟流上方形成穩(wěn)定的巖漿通道時(shí)的過程，模擬的是底辟流最后階段快速上涌的過程.

3.2 柴北緣都蘭南部區(qū)域底辟流上涌模擬

根據(jù)柴北緣都蘭南部區(qū)域超高壓變質(zhì)巖石上涌現(xiàn)代觀測數(shù)據(jù)(Song et al., 2003)，采用和觀測值一致的底辟流初始上涌溫度，并且為了進(jìn)一步解釋說明多個(gè)底辟流上涌對實(shí)驗(yàn)結(jié)果的影響，設(shè)計(jì)雙底辟流上涌對照組模型(如圖15、16)，初始半徑均設(shè)置為10 km.底辟流上涌的初始溫度略低于周圍巖體的溫度，上涌時(shí)伴隨著劇烈熱交換，自身溫度會(huì)先升高，由于底辟流本身的黏度和密度值比周圍巖體小，溫度升高導(dǎo)致熔融比例加大，其黏度和密度會(huì)進(jìn)一步減小，由此導(dǎo)致侵蝕現(xiàn)象的發(fā)生.模型結(jié)果顯示，在相同的時(shí)間內(nèi)(3 Ma)，雙底辟流上涌時(shí)，底辟流向上侵蝕運(yùn)移的距離大于單個(gè)底辟流上涌距離，進(jìn)一步說明了當(dāng)持續(xù)增量的底辟流上涌時(shí)，底辟流向上侵蝕運(yùn)移的距離會(huì)不斷加大，最終形成穩(wěn)定的巖漿通道，從而形成出露地表的巖漿弧.

圖15 柴北緣都蘭南部區(qū)域雙底辟流上涌初始模型(a) 初始溫度示意圖； (b) 初始黏度示意圖； (c) 初始密度示意圖.Fig.15 Initial model of upwelling of double diapirs in South-Dulan region of North Qaidam(a) Initial temperature; (b) Initial viscosity; (c) Initial density.

圖16 不同時(shí)刻柴北緣都蘭南部區(qū)域雙底辟流上涌模型圖(a) 1.0 Ma； (b) 2.0 Ma； (c) 3.0 Ma.Fig.16 Diagram of the upwelling model of double diapir flow in South-Dulan region of North Qaidam at different times

3.3 與前人研究結(jié)果對比

模型一主要研究了單個(gè)底辟流上涌時(shí)單一變量半徑對上涌速度的影響，實(shí)驗(yàn)結(jié)果表明半徑越大，上涌速度越快.但實(shí)際情況是在海洋板塊俯沖時(shí)由于自身脫水熔融作用在地幔深度約90 km左右會(huì)導(dǎo)致大量的底辟上涌，這些底辟有可能是單個(gè)、多個(gè)，或者持續(xù)增量的形式同時(shí)上涌，而多個(gè)底辟上涌過程中由于大的底辟上涌速度快，在上涌過程中可能會(huì)追上半徑小的底辟，兩者產(chǎn)生黏結(jié)作用后，以兩者黏結(jié)前的速度中間值向地殼運(yùn)移.

針對單一底辟流的上涌速率，我們根據(jù)Hadamard Rybczynski方程計(jì)算牛頓流球體穿過冪次定律流進(jìn)行比較分析底辟流上涌時(shí)的速度(Weinberg and Podladchikov, 1994).

(18)

其中μeff是周圍環(huán)境流體的有效黏度，這里我們采用的常黏度進(jìn)行計(jì)算，取值1020 Pa·s，A是材料常數(shù)，μsph是球體的牛頓黏度，我們這里為了簡化模型采用的也是常黏度，取1018Pa·s，r是流體半徑，Δρ是流體與周圍地幔密度差，參數(shù)n是常數(shù)2(Miller et al.1988)，參數(shù)Xsol，M和G是m=1/n的函數(shù).

Xsol=1.3(1-m2)+m,

M=0.76+0.24m,

G=2.39-5.51m+3.77m2.

(19)

Xsol用于糾正流體的速度，Xsol,M,G都是和材料相關(guān)的參數(shù)，三個(gè)參數(shù)配合起來用于對流體下降速度的糾正.代入公式計(jì)算可得到Xsol=1.48,M=0.88,G=0.76.根據(jù)Hadamard Rybczynski公式計(jì)算出來的底辟流上升速度理論值為V=6.67 cm·a-1.通過實(shí)驗(yàn)?zāi)M三組不同半徑的底辟流上涌模型得到相關(guān)的實(shí)驗(yàn)結(jié)果為：半徑為5 km的底辟流平均上涌速度2.73 cm·a-1, 半徑為10 km的平均上涌速度3.33 cm·a-1, 半徑為15 km的平均上涌速度4 cm·a-1.實(shí)驗(yàn)結(jié)果和理論值在數(shù)量級(jí)上一致，由于考慮到地幔對流的影響，三組不同半徑的底辟流上涌模型模擬值略低于理論值.

另外，針對Hall和Kincaid(2001)提出的底辟流通過地幔楔的機(jī)理，我們的實(shí)驗(yàn)結(jié)果也進(jìn)行了相關(guān)對比分析.Hall和Kincaid(2001)提出的流體上涌方式和一個(gè)無量綱的常數(shù)Bn值有關(guān)，Bn值的計(jì)算方法如下：

(20)

Hall和Kincaid(2001)提供了一個(gè)針對海洋俯沖板塊流體上涌的概念框架，并且通過模型實(shí)驗(yàn)得出不同的Bn值范圍對應(yīng)不同上涌模式，但由于實(shí)驗(yàn)過程中采用的是等黏度，等溫度模型，忽略了模型上涌過程中溫度對流體黏度以及熔融的影響，與現(xiàn)實(shí)條件下底辟流上涌過程黏度隨溫度實(shí)時(shí)變化不符.另外，Hall和Kincaid(2001)模型中底辟流上涌速度假定是均勻的，存在以上兩點(diǎn)假設(shè)，而我們的模型設(shè)計(jì)做了相應(yīng)的改進(jìn)，不僅考慮溫度場對于底辟上涌過程密度和黏度的變化，上涌速度也是受到黏度和密度改變產(chǎn)生實(shí)時(shí)的變化.同時(shí)考慮了底辟流上涌過程中各個(gè)階段，各圈層之間巨大的黏滯度差異、復(fù)雜的巖體黏-彈-塑性本構(gòu)關(guān)系(Hall et al., 2003).

將本文模型二的實(shí)驗(yàn)參數(shù)帶入公式(20)計(jì)算所得的Bn值屬于高閾值范圍.根據(jù)模型二實(shí)驗(yàn)結(jié)果，單一底辟流無法以固定的直徑上涌，而是在上涌過程中半徑不斷增大，并伴隨物質(zhì)擴(kuò)散直至消亡在地幔中.所以Hall和Kincaid(2001)提出的針對上述不同的Bn值范圍，底辟流對應(yīng)不同上涌模式的結(jié)論還應(yīng)結(jié)合實(shí)際溫度場影響條件下，進(jìn)一步的完善與討論.

板塊俯沖導(dǎo)致的大量底辟流上涌，巖漿的持續(xù)上涌并由此形成的巖漿通道是導(dǎo)致地殼增厚熔融的重要因素(吳福元等，2014).在底辟流上涌的初始階段，單個(gè)底辟流上涌最終不能上涌至下地殼說明單個(gè)板片熔體不足以透過地幔巖石層向淺部侵蝕蔓延.因此，我們認(rèn)為實(shí)際情況可能是由于隨著祁連海洋板塊不斷俯沖至地幔深處，自身脫水熔融產(chǎn)生持續(xù)增量的底辟不斷上涌導(dǎo)致的脫水熔融作用會(huì)使地幔楔巖體的熔點(diǎn)降低，并且伴隨著巨量熔體向上遷移進(jìn)而產(chǎn)生大量幔源巖漿；另一方面，板片富水流體或含水熔體的加入導(dǎo)致熔融流體在向上遷移的過程中極大降低地幔楔內(nèi)巖石的流變強(qiáng)度和黏滯度，并強(qiáng)烈擾動(dòng)地幔對流場，進(jìn)而導(dǎo)致對巖石圈地幔的熱侵蝕作用以及多種構(gòu)造變形的發(fā)生.在地幔對流的共同作用下沿著巖石層的軟弱帶和斷層裂隙不斷向大陸方向融熔侵蝕，最終形成直達(dá)地表淺層的巖漿通道，即安第斯模型中的底辟流模型是最終導(dǎo)致祁連弧形成的重要機(jī)理(Yin et al., 2012).巖漿通道形成以后，隨著海洋板塊的俯沖，冷的變質(zhì)巖會(huì)沿著巖漿通道快速上涌出露地表，形成超高壓變質(zhì)帶.

3.4 現(xiàn)代觀測數(shù)據(jù)

根據(jù)綠梁山橄欖巖測定的SHRIMP鋯石U-Pb年齡數(shù)據(jù)，超高壓變質(zhì)巖產(chǎn)生在488～495 Ma(Zhang et al., 2005)，同樣的研究表明，錫鐵山超高壓變質(zhì)巖產(chǎn)生在478～488 Ma(Zhang et al., 2005)，都蘭區(qū)域產(chǎn)生在422～450 Ma(Song et al., 2005; Mattinson et al., 2006)，上述數(shù)據(jù)表明柴達(dá)木北部超高壓變質(zhì)作用可能持續(xù)了數(shù)十百萬年.稀土元素研究表明，綠梁山地區(qū)的石榴石橄欖巖來源于海洋地幔，這些巖石經(jīng)歷了早古生代的變質(zhì)作用(457±22 Ma)(Song et al., 2006).通過檢測榴輝巖的原巖和地質(zhì)年代，Zhang等(2005) 認(rèn)為綠梁山和錫鐵山具有截然相反的減壓路徑：前者減壓時(shí)間非常短，后者則在下地殼中停留了相對較長的時(shí)間(>80 Ma).減壓時(shí)間短可能是由于綠梁山地區(qū)相比錫鐵山地區(qū)地殼巖石層軟弱帶更集中，或者該區(qū)域巖漿上涌量巨大，侵蝕作用非常劇烈，從流體上涌至形成的巖漿通道經(jīng)歷的時(shí)間較短，那么對應(yīng)的巖漿通道形成后變質(zhì)巖石減壓出露地表的時(shí)間會(huì)相對短.

白云母的40Ar/39Ar熱力學(xué)分析表明綠梁山超高壓變質(zhì)巖從465～390 Ma運(yùn)移了<10 km的距離，即超高壓變質(zhì)巖在中地殼停留了超過70 Ma.我們的模擬結(jié)果表明，對于已經(jīng)形成巖漿通道的單個(gè)底辟或者超高壓變質(zhì)巖而言，向上減壓運(yùn)移的時(shí)間就是整個(gè)流體上涌過程經(jīng)歷的時(shí)間，該過程只需要數(shù)個(gè)百萬年就可以完成，而超高壓變質(zhì)巖在中地殼中停留時(shí)間漫長則是由于中地殼至上地殼的上涌通道并未形成.因此，巖漿上涌通道的形成無論是對于超高壓變質(zhì)巖出露地表所經(jīng)歷的減壓時(shí)間還是對于地表巖漿弧的形成時(shí)間均具有重要的影響.上涌通道形成時(shí)間越早，巖漿弧出露地表的年代越久遠(yuǎn)，同樣對應(yīng)的超高壓變質(zhì)巖石減壓時(shí)間越短.而從底辟流熔體開始上涌至上涌通道的形成時(shí)間則與俯沖區(qū)域上覆巖層的軟弱帶以及上涌流體的體積量有關(guān)，軟弱帶越集中，形成時(shí)間越短；上涌流體的體積越大，即半徑越大，上涌速度越快，流體地幔對流越劇烈，則垂向侵蝕的速度會(huì)越快對應(yīng)巖漿通道的形成時(shí)間越短.

柴達(dá)木北緣超高壓變質(zhì)巖地體榴輝巖的原巖可能來源于兩個(gè)途徑：一種是沿著大陸裂隙運(yùn)移至地表的鎂鐵質(zhì)高溫巖墻形成的火成巖，另一種是在島弧或洋中脊擴(kuò)張中心產(chǎn)生的大洋板塊碎片.綠梁山和錫鐵山榴輝巖Nd同位素?cái)?shù)據(jù)表明，其超高壓變質(zhì)巖石來源于高度虧損地幔，現(xiàn)代觀測表明榴輝巖來源于弧后和大洋中脊，而非來源于大陸(Zindler and Hart, 1986).柴北緣都蘭區(qū)域的稀土元素研究表明榴輝巖的原巖為洋中脊玄武巖或海山玄武巖(Song et al., 2003).以上元素測定結(jié)果支持柴達(dá)木北緣祁連弧底辟流上涌模式為具有穩(wěn)定巖漿通道的底辟流模型(Yin et al., 2012).因此，安第斯型底辟流模型可以較好地解釋早古生代青藏高原柴北緣祁連弧的形成過程.

4 結(jié)論

本文針對前人對超高壓變質(zhì)巖形成機(jī)制的不同地質(zhì)解釋，從簡單模擬到綜合模擬，建立三組數(shù)值模型，考慮了底辟流的半徑、溫度耦合、上涌通道等因素對超高壓變質(zhì)巖的形成過程進(jìn)行了系列數(shù)值模擬實(shí)驗(yàn)，研究結(jié)果對早古生代柴北緣祁連弧的形成過程有一定的啟發(fā)意義.取得了幾點(diǎn)認(rèn)識(shí)：

(1)單個(gè)底辟流很難直接上涌至下地殼，早古生代柴北緣超高壓變質(zhì)地體和巖漿弧可能是由于多個(gè)底辟流持續(xù)上涌，最終發(fā)育穩(wěn)定巖漿通道形成的.巖漿通道的形成對于超高壓變質(zhì)巖石的減壓時(shí)間以及巖漿弧的形成時(shí)間均具有重要影響.

(2)單個(gè)底辟很難直接上涌至下地殼，柴達(dá)木北緣超高壓變質(zhì)地體和巖漿弧可能是由于多個(gè)底辟流持續(xù)上涌，最終發(fā)育穩(wěn)定巖漿通道形成的，并且?guī)r漿通道的形成對于超高壓變質(zhì)巖石的減壓時(shí)間以及巖漿弧的形成時(shí)間均具有重要影響.

(3)具有穩(wěn)定巖漿通道的單個(gè)底辟流從地幔深處90 km上涌至殼幔邊界的過程中，在水平方向的侵蝕范圍要比垂向侵蝕范圍大一倍左右.

(4)安第斯型底辟流模型可以較好地解釋早古生代青藏高原柴北緣祁連弧的形成過程.

致謝感謝審稿專家對本文提出寶貴修改意見，同時(shí)感謝皇甫鵬鵬、王少坡、徐佳靜對本文的建議和幫助.