李忠海，石耀霖

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

?

三維板塊幾何形態(tài)對(duì)大陸深俯沖動(dòng)力學(xué)的制約

李忠海，石耀霖

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

大陸深俯沖及超高壓變質(zhì)作用是大陸動(dòng)力學(xué)的重要研究?jī)?nèi)容，前人進(jìn)行了系統(tǒng)的地質(zhì)、地球物理觀測(cè)以及數(shù)值模擬研究.然而，自然界中大陸板塊的俯沖、碰撞及造山過程大部分具有明顯的沿走向的差異性，這種典型的三維特征可能很大程度上依賴于會(huì)聚大陸板塊的初始幾何學(xué)和運(yùn)動(dòng)學(xué)特征.本文采用三維高分辨率的動(dòng)力學(xué)數(shù)值模擬方法，建立了方形大陸板塊和楔形大陸板塊兩種不同的俯沖-碰撞模型，并且俯沖大陸板塊側(cè)面與大洋俯沖帶相鄰.數(shù)值模擬結(jié)果揭示大洋板塊可以持續(xù)地俯沖到地幔之中，而大陸板塊俯沖到一定深度處，其前端的俯沖板塊將發(fā)生斷離，并進(jìn)而造成殘余的大陸板塊俯沖角度的減小.方形大陸俯沖板塊的斷離深度約為150 km，而楔形大陸俯沖板塊的斷離深度較大，約250～300 km，這很大程度上取決于俯沖帶中大洋板塊的牽引力和大陸板塊的負(fù)浮力之間的競(jìng)爭(zhēng)關(guān)系.同時(shí)，無(wú)論方形還是楔形大陸板塊俯沖模型中，板塊斷離后，側(cè)向的大洋俯沖板塊仍可以拖曳約60～70 km寬的大陸邊緣巖石圈持續(xù)向下俯沖，揭示了新西蘭東部的洋-陸空間轉(zhuǎn)換俯沖帶的動(dòng)力學(xué)機(jī)制.并且，數(shù)值模型與喜馬拉雅造山帶和秦嶺—大別—蘇魯造山帶進(jìn)行了對(duì)比，進(jìn)而對(duì)其高壓-超高壓巖石空間展布沿走向的差異性特征和機(jī)制提供了一定的啟示.

大陸俯沖；超高壓變質(zhì)；板塊幾何形態(tài)；數(shù)值模擬；喜馬拉雅；大別—蘇魯

1　引言

在大陸板塊會(huì)聚過程中，大陸(邊緣)通常緊接著大洋俯沖的完結(jié)而開始俯沖和碰撞，大量密度較小的大陸地殼物質(zhì)進(jìn)入俯沖隧道，從而導(dǎo)致會(huì)聚速率的降低，并最終可能停止.在大陸碰撞的早期階段，最重要的事件之一就是高壓-超高壓變質(zhì)巖石的形成和折返，這是過去三十年里地球科學(xué)領(lǐng)域最重大的發(fā)現(xiàn)之一.這些發(fā)現(xiàn)揭示大陸地殼的物質(zhì)可以深俯沖至100 km以下，經(jīng)歷超高壓變質(zhì)作用，而后折返至地表(Chopin，2003；Zheng，2012).至今，全球的超高壓變質(zhì)帶已經(jīng)確認(rèn)的有20余處，主要分布于從歐洲的西阿爾卑斯山到中國(guó)大別—蘇魯造山帶的廣大地區(qū)(Liou et al.，2004；楊經(jīng)綏等，2009；鄭永飛，2009).但是，并不是所有的板塊碰撞造山帶都存在超高壓變質(zhì)巖；即使對(duì)于同一個(gè)造山帶，超高壓變質(zhì)巖也經(jīng)常只在特定的區(qū)域出現(xiàn)，例如喜馬拉雅造山帶的西端、阿爾卑斯造山帶的西南部以及秦嶺—大別造山帶東部的大別—蘇魯?shù)貐^(qū)等.

對(duì)于大陸深俯沖的研究，在系統(tǒng)的地質(zhì)學(xué)和地球物理學(xué)觀測(cè)基礎(chǔ)上，數(shù)值模擬越來越成為研究該地球動(dòng)力學(xué)過程的重要手段.前人的數(shù)值模型大多是基于垂直造山帶走向的二維模型(綜述如李忠海，2014).這些數(shù)值模型對(duì)大陸深俯沖和折返動(dòng)力學(xué)提供了一系列定量化的約束，從而極大地促進(jìn)了人們對(duì)于該過程基本動(dòng)力學(xué)機(jī)制的理解.然而，自然界中的大陸俯沖-碰撞帶很多都具有不規(guī)則的三維結(jié)構(gòu)，譬如喜馬拉雅造山帶和秦嶺—大別—蘇魯造山帶等.對(duì)于喜馬拉雅造山帶，大印度板塊最初的幾何學(xué)和運(yùn)動(dòng)學(xué)特征，以及與亞洲大陸的最終碰撞時(shí)間尚存比較大的爭(zhēng)議(Ali and Aitchinson，2005；Guillot et al.，2003).一個(gè)比較重要的推測(cè)在于印度板塊的西緣率先與亞洲碰撞，而后逐漸向東部的喜馬拉雅腹地拓展(Rowley，1996；Treloar and Coward，1991).對(duì)于秦嶺—大別—蘇魯造山帶，一般認(rèn)為揚(yáng)子大陸板塊東部首先向華北板塊之下俯沖，而后向西呈剪刀狀閉合(如，馬醒華和楊振宇，1993；Yin and Nie,1993).這些陸-陸碰撞的大致演化歷史表明，兩個(gè)大陸板塊的會(huì)聚并不一定沿著現(xiàn)今的縫合帶同時(shí)、同速俯沖和碰撞；相反，往往在某一位置首先接觸，而后逐步擴(kuò)展.這種三維特征的幾何學(xué)和運(yùn)動(dòng)學(xué)特征是二維模型所無(wú)法體現(xiàn)的，它也肯定對(duì)于大陸深俯沖及后期高壓-超高壓巖石的折返具有一定程度的制約.

至今，關(guān)于大陸深俯沖的相關(guān)三維數(shù)值模型仍比較欠缺，因此本文旨在建立不同幾何學(xué)特征的三維大陸板塊俯沖-碰撞模型，探討方形大陸板塊俯沖和楔形大陸板塊俯沖過程的共性以及差異性，進(jìn)而與自然界中的特定俯沖-碰撞帶進(jìn)行對(duì)比研究.

2　數(shù)值模擬方法

2.1力學(xué)和熱學(xué)控制方程

對(duì)于板塊俯沖-碰撞相關(guān)的動(dòng)力學(xué)數(shù)值模型，一般對(duì)三組控制方程進(jìn)行求解，包括斯托克斯流體動(dòng)力學(xué)方程、物質(zhì)守恒方程以及熱量守恒方程.本文中，將對(duì)這些方程采用有限差分算法進(jìn)行求解(Gerya，2010).

(1)斯托克斯方程：

(1)

(2)不可壓縮流體的物質(zhì)守恒方程

(2)

(3)熱量守恒方程

(3)

2.2黏-塑性流變學(xué)性質(zhì)

數(shù)值模型中，偏應(yīng)力張量和應(yīng)變率張量的流變學(xué)關(guān)系采用整合的黏-塑性本構(gòu)關(guān)系.對(duì)于不可壓縮流體的黏滯性變形，流變學(xué)方程為：

(4)

這里，ηeff表示等效黏滯系數(shù)，它是一個(gè)關(guān)于溫度、壓力、物質(zhì)成分、應(yīng)變率以及部分熔融程度的函數(shù).

韌性流變的黏滯系數(shù)定義為：

(5)

上述的韌性流變需要與塑性流變相結(jié)合,從而形成實(shí)際的黏-塑性流變特征.為此，這里采用改進(jìn)的Drucker-Prager屈服準(zhǔn)則(Ranalli,1995)：

(6)

該方程組中，σyield表示屈服應(yīng)力.P是壓力.C0是P=0條件下的巖石剩余強(qiáng)度，φ是內(nèi)摩擦角.λ是孔隙流體系數(shù)，它控制著含流體孔隙介質(zhì)的脆性強(qiáng)度.因此，φeff可以解釋為等效內(nèi)摩擦角，它結(jié)合了內(nèi)摩擦角(φ)和孔隙流體系數(shù)(λ)的作用.

基于ηductile和ηplastic，黏-塑性流變關(guān)系的最終黏滯系數(shù)可以定義為韌性和塑性黏滯系數(shù)中的較小者(Ranalli，1995)：

(7)

2.3脫水-水化作用模型

俯沖帶中的水對(duì)于大洋俯沖和大陸碰撞動(dòng)力學(xué)具有非常重要的影響(李忠海等，2015；李忠海和許志琴，2015).數(shù)值模型中的水主要包括兩部分.一種是自由水，一種是礦物水.對(duì)于自由水，假設(shè)沉積物和大洋上地殼玄武巖在地表含2wt.%(質(zhì)量百分含量)的水，并且隨著深度的增加而減少，直至75 km的自由水含量降為零(Li et al.，2013).對(duì)此,假設(shè)線性的定量化關(guān)系式：

XH2O(P)(wt.%)=XH2O(P0)(1-Δz/75),

(8)

其中，XH2O(P0)代表地表的自由水含量2wt.%，Δz是深度(0～75 km)，對(duì)于其他的巖石屬性(大陸地殼，大洋下地殼，地幔)，假設(shè)其自由水含量為零.

對(duì)于礦物水，采用了PerpleX(Connolly，2005)所計(jì)算的四種典型巖性的含水量隨溫度和壓力的變化，這些數(shù)據(jù)是依據(jù)吉布斯自由能最小化原理，對(duì)特定巖石成分和熱力學(xué)數(shù)據(jù)計(jì)算而得(Gerya and Meilick，2011).隨著俯沖的進(jìn)行,礦物水析出后保存在新生成的示蹤點(diǎn)中，并獨(dú)立地運(yùn)移，當(dāng)其遇到由于水化或部分熔融作用而產(chǎn)生的水不飽和巖性的時(shí)候，該活動(dòng)的水就會(huì)被吸收.水的運(yùn)動(dòng)速度由該位置處的壓強(qiáng)梯度進(jìn)行計(jì)算，具體算法如下：

(9)

其中，vx(water)、vy(water)、vz(water)分別是水在x、y和z方向的速度.vx、vy、vz是該位置處的巖石速度.vy(percolation)代表流體-熔體向上的滲透速度，根據(jù)前人的研究，參考模型中采用vy(percolation)=10 cm·a-1(Gorczyk et al.，2007；Peacock，1990).在這里，模型假設(shè)水生成之后的水平速度與巖石速度相當(dāng)，由于水和巖石的密度差，導(dǎo)致其速度差異主要體現(xiàn)在垂向上.

2.4巖石部分熔融模型

數(shù)值模型中包含了多種巖石類型的部分熔融計(jì)算(Li et al.，2013).對(duì)于水化蝕變的地幔橄欖巖，采用了Katz等(2003)的算法和參數(shù).對(duì)于其他巖石類型，基于實(shí)驗(yàn)巖石學(xué)的約束條件，作為一種近似算法，假設(shè)部分熔融的體積比例與溫度之間存在一種線性關(guān)系：

(10)

其中，Tsolidus和Tliquidus分別代表特定巖性的固相線溫度和液相線溫度.

3　初始模型和邊界條件

自然界中，大陸板塊往往具有橫向不均一性，其中一個(gè)最顯著的特征在于幾何形態(tài)的復(fù)雜性，這對(duì)于大陸俯沖和碰撞的過程及地質(zhì)響應(yīng)具有重要的影響，因此這也是本文擬探討的核心問題.基于要研究的不同形態(tài)大陸板塊俯沖-碰撞動(dòng)力學(xué)背景，我們?cè)O(shè)計(jì)了兩組三維高分辨率的動(dòng)力學(xué)數(shù)值模型(圖1).初始模型中，上覆板塊是均一的大陸板塊，而俯沖板塊沿走向平分為兩個(gè)部分，一側(cè)為均一的大洋板塊，而另一側(cè)的大洋板塊之后與大陸板塊相連.兩組模型的區(qū)別在于大陸俯沖板塊的形態(tài)，一組模型中采用方形大陸俯沖板塊(圖1a)，而另一組模型中采用楔形大陸俯沖板塊(圖1b).盡管幾何學(xué)形態(tài)不同，但所有大陸板塊的上地殼厚度為20 km，下地殼厚度為15 km.大洋板塊的上地殼玄武巖厚度為3 km，下地殼輝長(zhǎng)巖厚度為5 km.在初始模型中，大洋板塊和上覆大陸板塊的結(jié)合部設(shè)計(jì)了一個(gè)傾角約為25°的初始薄弱帶(圖1)，該薄弱帶用于限定初始大洋俯沖的發(fā)生位置，此后的俯沖狀態(tài)和過程將根據(jù)模型的動(dòng)力學(xué)演化而進(jìn)行自適應(yīng)的調(diào)整.值得注意的是，在地殼表面之上，與自由滑動(dòng)的模型頂界面之間，設(shè)計(jì)有一層相對(duì)高黏滯度的偽空氣層，其與上地殼的接觸面用以模擬地形起伏面(Schmeling et al.,2008；Li et al.,2013).對(duì)于大陸巖石圈的溫度結(jié)構(gòu)，地表是0 ℃，巖石圈底部約為1345 ℃，其間采用線性插值的方法.大洋巖石圈的年齡約為30 Ma.巖石圈之下的軟流圈地幔的溫度梯度為0.5 ℃/km.數(shù)值模型中采用的各種不同巖石類型的屬性見表1和表2.

圖1　初始模型設(shè)計(jì)(a)和(b)分別為方形大陸板塊俯沖和楔形大陸板塊俯沖的三維初始模型，模型的空間尺度為1000 km×680 km×656 km，分辨率為2 km×2 km×2 km.為了更加清楚地顯示內(nèi)部結(jié)構(gòu)，最頂層的20 km被剪切掉.顏色代表不同的巖石類型：0—偽空氣層；1—水；2/3—沉積物；4—部分熔融沉積物；5/6—大陸上、下地殼；7/8—水化的大陸上、下地殼；9/10—部分熔融的大陸上、下地殼；11/12—大洋上、下地殼；13—部分熔融的大洋地殼；14—巖石圈地幔；15—軟流圈地幔；16/17—水化、蛇紋石化的地幔；18—部分熔融的地幔.水化的和部分熔融的巖石類型在圖1的初始模型中不存在，將隨著模型的演化而產(chǎn)生.Fig.1　Initial model configuration(a)and (b)are the initial models with rectangular and wedge-shaped continental subducting plates,respectively.The spatial scale of 3-D model domain is 1000 km×680 km×656 km with the resolution of 2 km×2 km×2 km.The top layer (y>-20 km)is cut off for clarity.Colors represent different rock types,with:0—air;1—water;2/3—sediment;4—partial molten sediment;5/6—upper/lower continental crust;7/8—hydrated upper/lower continental crust;9/10—partial molten upper/lower continental crust;11/12—upper/lower oceanic crust;13—partial molten oceanic crust;14—lithospheric mantle;15—asthenospheric mantle;16/17—hydrated/serpentinized mantle;18—partially molten mantle.The hydrated and partially molten rocks are not shown in Fig.1,but will appear during the evolution of the model.

表1　數(shù)值模型采用的黏滯性流變參數(shù)

注：η0是黏滯系數(shù)初始參考值，可通過下式計(jì)算:η0=(1/AD)×10(6n).部分熔融酸性巖(RF)的流變性質(zhì)用于模型中部分熔融的大陸地殼和沉積物.部分熔融基性巖(RG)的流變性質(zhì)用于模型中部分熔融的大洋地殼和地幔物質(zhì).參考文獻(xiàn)：Kirby，1983；Kirby 和 Kronenberg，1987；Ranalli 和 Murthy，1987；Ji 和 Zhao，1993；Ranalli，1995.

表2　數(shù)值模型采用的物質(zhì)參數(shù)

注：不同巖石類型的數(shù)字標(biāo)號(hào)如圖1所示.所有巖石具有相同的熱膨脹系數(shù)(α=3×10-5K-1)和壓縮系數(shù)(β=1×10-5MPa-1).熱傳導(dǎo)系數(shù)計(jì)算如下：K1=[0.64+807/(TK+77)]exp(0.00004PMPa)；K2=[1.18+474/(TK+77)]exp(0.00004PMPa)；K3=[0.73+1293/(TK+77)]exp(0.00004PMPa).巖石固相線溫度計(jì)算如下：TS1={889+17900/(P+54)+20200/(P+54)2,P<1200 MPa} 或 {831+0.06P,P>1200 MPa}；TS2={973-70400/(P+354)+778×105/(P+354)2,P<1600 MPa}或{935+0.0035P+0.0000062P2,P>1600 MPa}；TS3=KATZ2003.巖石液相線溫度計(jì)算如下：TL1=1262+0.09P；TL2=1423+0.105P；TL3=KATZ2003.巖石的黏滯性流變性質(zhì)如表1所示.模型中塑性形變主要受控于等效內(nèi)摩擦角.參考文獻(xiàn)分別為:1=Turcotte和Schubert，1982；2=Bittner和Schmeling，1995；3=Clauser和Huenges，1995；4=Ranalli，1995；5=Schmidt和Poli，1998；6=Katz等，2003.

兩組數(shù)值模型的邊界條件完全相同.對(duì)于速度邊界條件，頂部(y=0)和前后兩個(gè)側(cè)面(z=0和z=656 km)是自由滑動(dòng)邊界.底部邊界是滲透性邊界，采用近無(wú)限深度的外部自由滑動(dòng)邊界條件(Li et al.，2013).模型的左右兩個(gè)邊界是會(huì)聚邊界，自初始條件開始，兩個(gè)邊界之間的會(huì)聚速率保持4 cm·a-1不變.對(duì)于溫度邊界條件，模型頂部為固定溫度(0 ℃)；四周的側(cè)面邊界條件為水平方向溫度梯度為零(即零熱流).底部邊界采用的是外部邊界固定溫度條件，即在模型底邊界(680 km)之下的200 km處，假設(shè)一個(gè)固定的地幔溫度，這樣就可以使得在底部滲透邊界處，溫度和熱流都可以據(jù)模型演化而動(dòng)態(tài)調(diào)整(Li et al.，2013).

4　模擬結(jié)果

4.1方形大陸板塊俯沖-碰撞模型

方形大陸板塊俯沖-碰撞的初始模型如圖1a所示，兩個(gè)大陸板塊之間的大洋板塊的長(zhǎng)度約為450 km，板塊的會(huì)聚速率保持4 cm·a-1不變.

該模型的演化結(jié)果揭示大洋俯沖帶和大陸碰撞帶產(chǎn)生了差異巨大的演化特征.在大陸碰撞帶，大陸巖石圈緊隨著大洋板塊的俯沖而進(jìn)入俯沖隧道(圖2a)，但是由于大陸地殼的密度較小，因此伴隨著大陸俯沖的進(jìn)行，一些地殼物質(zhì)會(huì)從俯沖板塊拆離并發(fā)生折返(圖2c).同時(shí)，先期俯沖的大洋板塊在碰撞后發(fā)生斷離，斷離深度約為150 km，斷離時(shí)間為20 Ma左右(圖2c).板塊斷離后，殘留的俯沖板塊向上彎曲，造成大陸板塊俯沖角度變小，并進(jìn)一步推動(dòng)了俯沖大陸地殼的折返(圖2e).在大洋俯沖一側(cè)，海溝持續(xù)地后撤，伴隨著俯沖大洋板塊的脫水，以及地幔楔中的水化和部分熔融作用(圖2b，2d，2f).同時(shí)，大洋海溝的后撤產(chǎn)生了弧后的擴(kuò)張，進(jìn)而促進(jìn)了上覆地殼物質(zhì)從陸-陸碰撞一側(cè)向洋-陸俯沖一側(cè)逃逸.

4.2楔形大陸板塊俯沖-碰撞模型

楔形大陸板塊俯沖-碰撞的初始模型如圖1b所示，兩個(gè)大陸板塊之間的大洋板塊的長(zhǎng)度在模型側(cè)面的最窄處約為250 km，在模型中部的最寬處約為600 km.與方形大陸板塊會(huì)聚模型相同，板塊的會(huì)聚速率仍然保持4 cm·a-1不變.

該楔形大陸板塊俯沖-碰撞的數(shù)值模擬結(jié)果也揭示了大洋俯沖一側(cè)和大陸碰撞一側(cè)截然不同的演化特征.大洋板塊持續(xù)地俯沖到地幔之中，伴隨著俯沖板塊的脫水以及地幔楔的部分熔融.同時(shí)，由于海溝的后撤導(dǎo)致上覆巖石圈的伸展以及弧后的擴(kuò)張(圖3b，3d，3f).在大陸碰撞一側(cè)，俯沖的大陸地殼物質(zhì)從不同的深度處折返(圖3a).楔形大陸板塊的犄角持續(xù)地深俯沖，直至約250～300 km深度處發(fā)生斷離(圖3c)，而后俯沖大陸板塊的角度變緩，并逐漸轉(zhuǎn)變?yōu)榈徒嵌认虏迥Ｊ?圖3e).楔形大陸板塊與方形大陸板塊俯沖的主要差別在于，前者楔形板塊犄角處的大陸板塊較窄，相對(duì)容易地俯沖到較大的深度，而后發(fā)生物質(zhì)的折返以及板塊的斷離(圖2，3).

4.3不同形態(tài)大陸板塊俯沖-碰撞帶的內(nèi)部結(jié)構(gòu)

圖2和圖3中兩種不同形態(tài)大陸板塊的俯沖-碰撞模型揭示二者在大尺度動(dòng)力學(xué)方面具有一定的共性，都展示了大洋板塊俯沖的連續(xù)性以及大陸碰撞之后俯沖板塊的斷離.但是，二者在大陸板塊深俯沖、巖石折返及板塊斷離的深度方面又具有很大的差異性.為了進(jìn)一步揭示這種深部結(jié)構(gòu)的差異性，我們將模型內(nèi)部結(jié)構(gòu)的演化進(jìn)行對(duì)比(圖4)，該圖中只顯示大洋板塊和大陸板塊的地殼物質(zhì)，而過濾掉地幔物質(zhì)，從而能夠更加清晰地展示模型深部結(jié)構(gòu)的演化特征.

從兩個(gè)模型內(nèi)部結(jié)構(gòu)演化的對(duì)比發(fā)現(xiàn)，方形大陸板塊的俯沖和斷離深度較淺，結(jié)合圖2可得，板塊的初始斷離深度約為150 km，并從模型側(cè)面邊緣向洋-陸轉(zhuǎn)化帶快速擴(kuò)展.而后大陸板塊俯沖角度變緩，從而形成一個(gè)斷離窗口(圖4a—4c).在楔形大陸板塊的俯沖-碰撞模型中，犄角處的大陸板塊可以俯沖至更大的深度處，而后發(fā)生斷離(250～300 km).斷離窗口沿著楔形的俯沖大陸板塊向內(nèi)并向上擴(kuò)展，從而中部洋-陸轉(zhuǎn)換區(qū)域的大陸俯沖和板塊斷離的深度較淺(圖4d—4f).但是兩個(gè)模型都顯示，板塊斷離并非沿著洋殼和陸殼的轉(zhuǎn)換帶而發(fā)生，而是在大陸板塊內(nèi)部的一定位置處，從而俯沖的大洋板塊可以拖曳一定寬度(約60～70 km)的大陸邊緣板塊持續(xù)向下俯沖(圖4).

5　討論

5.1大洋板塊拖曳窄條大陸板塊深俯沖

大陸深俯沖的驅(qū)動(dòng)力既包括大陸板塊前端已經(jīng)俯沖的大洋板塊的“拉力”和后端洋中脊的“推力”(Turcotte and Schubert，2002)，也包括三維空間中俯沖大陸板塊側(cè)面的大洋板塊的拉力.基于對(duì)后者進(jìn)行的數(shù)值模擬研究，前人提出了大洋板塊拖曳其側(cè)面的窄條陸殼俯沖的模型，并且得到當(dāng)陸殼厚度為30 km時(shí)，大洋俯沖可拖曳的大陸板塊極限寬度約為70～90 km；當(dāng)陸殼厚度為20 km時(shí)，可拖曳大陸板塊極限寬度為100～150 km(石耀霖和范桃園，2001，2002).我們新的數(shù)值模型中大陸地殼的厚度為35 km，包含20 km的上地殼和15 km的下地殼，同時(shí)模型包含了更加實(shí)際的地質(zhì)過程，例如脫水-水化作用、部分熔融作用等.模擬結(jié)果顯示，俯沖的大洋板塊可以拖曳約60～70 km寬的大陸板塊(無(wú)論是方形還是楔形)持續(xù)向下俯沖(圖4).這與前人的模型結(jié)果具有比較好的吻合(石耀霖和范桃園，2001).

圖2　方形大陸板塊俯沖-碰撞數(shù)值模型顏色代表不同的巖石類型，色標(biāo)如圖1.(a，c，e)大陸碰撞一側(cè)的模型演化；(b，d，f)大洋俯沖一側(cè)的模型演化.Fig.2　Numerical model with rectangular subducting continental plateColors represent different rock types as in Fig.1.(a，c，e)Evolution of the continental collision side；(b，d，f)Evolution of the oceanic subduction side.

圖3　楔形大陸板塊俯沖-碰撞數(shù)值模型顏色代表不同的巖石類型，色標(biāo)如圖1.(a，c，e)大陸碰撞一側(cè)的模型演化；(b，d，f)大洋俯沖一側(cè)的模型演化.Fig.3　Numerical model with wedge-shaped subducting continental plateColors represent different rock types as in Fig.1.(a，c，e)Evolution of the continental collision side；(b，d，f)Evolution of the oceanic subduction side.

圖4　方形和楔形大陸板塊俯沖-碰撞模型的內(nèi)部結(jié)構(gòu)演化對(duì)比顏色代表不同的巖石類型，色標(biāo)如圖1.(a—c)方形大陸板塊俯沖碰撞模型；(d—f)楔形大陸板塊俯沖碰撞模型.Fig.4　Comparison of internal structures between subduction models with rectangular and wedge-shaped continental platesColors represent different rock types as in Fig.1.(a—c)Model evolution of rectangular continental plate subduction；(d—f)Model evolution of wedge-shaped continental plate subduction.

圖5　太平洋板塊向新西蘭島之下俯沖的構(gòu)造簡(jiǎn)圖太平洋板塊主要在新西蘭北島和南島北端向澳大利亞板塊下俯沖.俯沖帶南端的Chatham隆起雖然在海平面下，但為大陸地殼，而其北面為典型的大洋地殼.地震P波速度成像的觀測(cè)結(jié)果也揭示南島北端的俯沖帶為低速的陸殼俯沖(Eberhart-Phillip and Reyners，1997).從而，該地區(qū)是一個(gè)從大洋俯沖帶沿走向過渡為大陸俯沖-碰撞帶的典型實(shí)例.Fig.5　Tectonic sketch of Pacific plate subduction under New Zealand islandThe Pacific plate subducts mainly under the North Island and northern part of the South Island of the New Zealand.The Chatham rise,located at the southern end of the subduction zone,is mainly composed of continental crust,although it is below the sea level.In contrast,the northern section of the subduction zone is definitely oceanic subduction.The seismic P-wave velocity structures also reveal the continental plate subduction under the northern part of the South Island (Eberhart-Phillip and Reyners,1997).Thereby,this is a typical region for along-strike transition from oceanic subduction to continental subduction and collision.

同時(shí)，這種大洋板塊拖曳其側(cè)面的窄條陸殼俯沖的模型在自然界中也有很好的觀測(cè)，例如新西蘭南島北端的陸殼俯沖帶(Allis，1981；Reyners and Cowan，1993；Eberhart-Phillip and Reyners,1997；Walcott，1998；石耀霖和范桃園，2001).圖5顯示在新西蘭島東海岸，太平洋的大洋巖石圈板塊向西俯沖到澳大利亞板塊之下，而大洋俯沖帶南端與大陸組分的Chatham隆起相鄰.冷而密度高的大洋俯沖板塊可以攜帶其南部東西向延伸的Chatham隆起及附近南北寬度約100多公里的窄條低密度陸殼巖石俯沖到新西蘭南島北端之下.相反，Chatham隆起更南的大陸地殼的主體飄浮在地幔之上而不能俯沖下去.因此，這里產(chǎn)生了一種特殊的現(xiàn)象，在南島北端的大陸俯沖帶與一般由于溫度較冷而地震波速較高的大洋板塊俯沖帶不同，該處大陸俯沖帶因陸殼波速低于地幔波速，因此俯沖帶反而呈低速(Eberhart-Phillip and Reyners,1997)(圖5b—5d).無(wú)論是多種獨(dú)立的數(shù)值模型，還是自然界中的觀測(cè)都證實(shí)大洋俯沖板塊可以拖曳一定寬度(幾十至一百多千米左右)的大陸(邊緣)板塊深俯沖至地幔之中.

5.2三維大陸板塊的幾何形態(tài)對(duì)深俯沖的制約

方形和楔形大陸板塊俯沖-碰撞模型的對(duì)比發(fā)現(xiàn)，在會(huì)聚速率等其他條件一定的情況下，楔形俯沖板塊的前端(即比較窄的大陸板塊)更容易俯沖到比較大的深度(圖4)，而后才發(fā)生斷離，這很大程度上歸功于與之相鄰的側(cè)向大洋板塊的牽引.在模型總寬度一定的情況下，如果大陸板塊的寬度較小，那么大洋板塊相對(duì)較寬；進(jìn)而俯沖洋-陸板塊所受到的牽引力較大，而負(fù)浮力的積累較小，從而可以俯沖到相對(duì)較大深度處.該實(shí)驗(yàn)結(jié)果表明，大陸板塊的俯沖、斷離和折返不僅受控于自身所在的二維垂向剖面的應(yīng)力演化，而且受控于三維空間中相鄰板塊的作用力.

自然界中的超高壓變質(zhì)巖石往往出露于造山帶的一端，而非沿整個(gè)造山帶展布，例如：喜馬拉雅造山帶的西構(gòu)造結(jié)附近、阿爾卑斯造山帶的西南端，挪威的西部片麻巖省，以及秦嶺—大別造山帶東部的大別和蘇魯?shù)貐^(qū)(Liou et al.,2004).這種分布特征與大陸俯沖板塊的幾何形態(tài)可能具有一定的關(guān)聯(lián)，下面我們以喜馬拉雅造山帶和秦嶺—大別造山帶的超高壓變質(zhì)巖石的分布為例進(jìn)行探討.

對(duì)于喜馬拉雅造山帶而言，前人推測(cè)印度板塊的西緣可能率先與亞洲碰撞，而后逐漸向喜馬拉雅腹地拓展(Rowley，1996；Treloar and Coward，1991).同時(shí)，整個(gè)喜馬拉雅造山帶已經(jīng)發(fā)現(xiàn)的超高壓變質(zhì)巖石都集中于造山帶的西緣，而整個(gè)喜馬拉雅腹地并無(wú)超高壓巖石發(fā)現(xiàn)(Yin，2006；Li et al.,2011).在這種構(gòu)造背景下，即使印度—亞洲會(huì)聚的速率保持恒定(如數(shù)值模型)，由于其側(cè)面大洋俯沖板塊的拖曳，西喜馬拉雅地區(qū)印度板塊也能夠俯沖到比較大的深度，從而相對(duì)容易形成超高壓變質(zhì)巖.同時(shí)，前人的數(shù)值模擬結(jié)果也表明，在較高的板塊會(huì)聚速率下，大陸巖石圈更容易俯沖至比較大的深度(Li et al.,2011)，因此如果伴隨印度—亞洲大陸碰撞的進(jìn)行，會(huì)聚速率降低，那么將進(jìn)一步導(dǎo)致率先高速俯沖碰撞的西喜馬拉雅地區(qū)相對(duì)于后期慢速俯沖的喜馬拉雅腹地更容易產(chǎn)生大陸深俯沖和超高壓折返.

對(duì)于秦嶺—大別—蘇魯造山帶，一般認(rèn)為揚(yáng)子大陸板塊東部首先向華北板塊之下俯沖，而后向西呈剪刀狀閉合(如，馬醒華和楊振宇，1993；Yin and Nie,1993).在此構(gòu)造背景之下，基于數(shù)值模擬的結(jié)果可得在東部地區(qū)大陸板塊可以俯沖到很大的深度處，而西部地區(qū)大陸俯沖深度較淺.該模型也可以比較好地解釋，東部蘇魯—大別地區(qū)產(chǎn)生大陸板塊的深俯沖以及超高壓變質(zhì)巖石的折返(Cong and Wang，1999；鄭永飛等,2009)，而在西部的秦嶺地區(qū)主要發(fā)生陸-陸碰撞造山，不形成超高壓變質(zhì)帶.

6　結(jié)論

本文中，我們建立了不同三維幾何形態(tài)(方形和楔形)的大陸板塊俯沖-碰撞數(shù)值模型.主要結(jié)論如下：

(1)洋-陸空間轉(zhuǎn)換數(shù)值模型揭示深部結(jié)構(gòu)沿走向的差異性，大陸碰撞帶發(fā)生俯沖板塊斷離，而大洋俯沖板塊持續(xù)下插.同時(shí)，地殼物質(zhì)發(fā)生從陸-陸碰撞帶向洋-陸俯沖帶的側(cè)向逃逸.

(2)相對(duì)于方形大陸，楔形大陸板塊的前端更容易俯沖到比較大的深度(250～300 km)，而后發(fā)生斷離；這很大程度上歸功于俯沖帶中大洋板塊較大的牽引力和大陸板塊相對(duì)較小的負(fù)浮力之間的競(jìng)爭(zhēng)效應(yīng).

(3)大陸碰撞帶板塊斷離后，側(cè)向的大洋俯沖板塊仍可以拖曳約60～70 km寬的大陸板塊(無(wú)論是方形還是楔形)持續(xù)向下俯沖.

(4)這些模擬結(jié)果有助于理解新西蘭南島北端大陸板塊俯沖帶與大洋板塊俯沖帶側(cè)面相接時(shí)，陸殼俯沖帶地震波速為低速的現(xiàn)象；顯示該處可能正是現(xiàn)今活動(dòng)的陸殼俯沖場(chǎng)所.同時(shí)，數(shù)值模型還對(duì)喜馬拉雅和秦嶺—大別造山帶中高壓-超高壓巖石展布沿走向的差異性特征和機(jī)制提供啟示.

致謝感謝張立飛老師提供的有益討論和幫助.感謝兩位評(píng)審人的建設(shè)性意見和建議.

References

Ali J R,Aitchison J C.2005.Greater India.Earth-Science Reviews,72(3-4):169-188.Allis R G.1981.Continental underthrusting beneath the Southern Alps of New Zealand.Geology,9(7):303-307.

Bittner D,Schmeling H.1995.Numerical modelling of melting processes and induced diapirism in the lower crust.Geophysical Journal International,123(1):59-70.

Chopin C.2003.Ultrahigh-pressure metamorphism:tracing continental crust into the mantle.Earth and Planetary Science Letters,212(1-2):1-14.

Clauser C,Huenges E.1995.Thermal conductivity of rocks and minerals.//Ahrens T J ed.Rock Physics and Phase Relations.Washington DC:American Geophysical Union,105-126.

Cong B L,Wang Q C.1999.The Dabie-Sulu UHP rocks belt:Review and prospect.Chinese Science Bulletin,44(12):1074-1085.

Connolly J A D.2005.Computation of phase equilibria by linear programming:A tool for geodynamic modeling and its application to subduction zone decarbonation.Earth and Planetary Science Letters,236(1-2):524-541.Eberhart-Phillips D,Reyners M.1997.Continental subduction and three-dimensional crustal structure:The northern South Island,New Zealand.Journal of Geophysical Research,102(B6):11843-11861.

Gerya T.2010.Introduction to Numerical Geodynamic Modelling.Cambridge:Cambridge University Press.

Gerya T V,Meilick F I.2011.Geodynamic regimes of subduction under an active margin:effects of rheological weakening by fluids and melts.Journal of Metamorphic Geology,29(1):7-31.

Gorczyk W,Willner A P,Gerya T V,et al.2007.Physical controls of magmatic productivity at Pacific-type convergent margins:Numerical modelling.Physics of the Earth and Planetary Interiors,163(1-4):209-232.

Guillot S,Garzanti E,Baratoux D,et al.2003.Reconstructing the total shortening history of the NW Himalaya.Geochemistry,Geophysics,Geosystems,4(7):doi:10.1029/2002GC000484.

Ji S C,Zhao P L.1993.Flow laws of multiphase rocks calculated from experimental data on the constituent phases.Earth and Planetary Science Letters,117(1-2):181-187.Katz R F,Spiegelman M,Langmuir C H.2003.A new parameterization of hydrous mantle melting.Geochemistry,Geophysics,Geosystems,4(9):1073.Kirby S H.1983.Rheology of the lithosphere.Review of Geophysics,21(6):1458-1487.Kirby S H,Kronenberg A K.1987.Rheology of the lithosphere:selected topics.Reviews of Geophysics,25(6):1219-1244.

Li Z H,Xu Z Q,Gerya T V.2011.Flat versus steep subduction:contrasting modes for the formation and exhumation of high- to ultrahigh-pressure rocks in continental collision zones.Earth and Planetary Science Letters,301(1-2):65-77.

Li Z H,Xu Z Q,Gerya T,et al.2013.Collision of continental corner from 3-D numerical modeling.Earth and Planetary Science Letters,380:98-111.

Li Z H.2014.A review on the numerical geodynamic modeling of continental subduction,collision and exhumation.Science China:Earth Sciences,57(4):47-69.Li Z H,Xu Z Q.2015.Dynamics of along-strike transition between oceanic subduction and continental collision:Effects of fluid-melt activity.Acta Petrologica Sinica (in Chinese),31(12):3524-3530.

Li Z H,Liu M Q,Gerya T.2015.Material transportation and fluid-melt activity in the subduction channel:numerical modeling.Science China:Earth Sciences,58(8):1251-1268.

Liou J G,Tsujimori T,Zhang R Y,et al.2004.Global UHP metamorphism and continental subduction/collision:The Himalayan model.International Geology Review,46(1):1-27.Ma X H,Yang Z Y.1993.The collision and suturing of the three major blocks in China and the reconstruction of the paleo-Eurasia continent.Chinese Journal of Geophysics (Acta Geophysica Sinica)(in Chinese),36(4):476-488.

Peacock S A.1990.Fluid processes in subduction zones.Science,248(4953):329-337.

Ranalli G,Murphy D C.1987.Rheological stratification of the lithosphere.Tectonophysics,132(4):281-295.

Ranalli G.1995.Rheology of the Earth.2nd ed.Netherlands:Springer.

Reyners M,Cowan H.1993.The transition from subduction to continental collision:crustal structure in the North Canterbury region,New Zealand.Geophysical Journal International,115(3):1124-1136.

Rowley D B.1996.Age of initiation of collision between India and Asia:A review of stratigraphic data.Earth and Planetary Science Letters,145(1-4):1-13.

Schmeling H,Babeyko A Y,Enns A,et al.2008.A benchmark comparison of spontaneous subduction models-Towards a free surface.Physics of the Earth and Planetary Interiors,171(1-4):198-223.

Schmidt M W,Poli S.1998.Experimentally based water budgets for dehydrating slabs and consequences for arc magma generation.Earth and Planetary Science Letters,163(1-4):361-379.Shi Y L,Fan T Y.2001.Analysis on maximum scale of continental crust sliver driven by subducting oceanic slab:A case study of south island of New Zealand and UHPM at Dabie,China.Chinese Journal of Geophysics (in Chinese),44(6):754-760.

Shi Y L,Fan T Y.2003.Continental sliver subduction driven by oceanic slab:tectonic settings for large scale ultra-high pressure metamorphism.Acta Geologica Sinica (in Chinese),76(1):77-82.Treloar P J,Corward M P.1991.Indian Plate motion and shape:constraints on the geometry of the Himalayan orogen.Tectonophysics,191(3-4):189-198.Turcotte D L,Schubert G.2002.Geodynamics.2nd ed.Cambridge:Cambridge University Press.Walcott R I.1998.Modes of oblique compression:Late Cenozoic tectonics of the south island of New Zealand.Reviews of Geophysics,36(1):1-26.

Yang J S,Xu Z Q,Zhang J X,et al.2009.Tectonic setting of main high-and ultrahigh-pressure metamorphic belts in China and adjacent region and discussion on their subduction and exhumation mechanism.Acta Petrologica Sinica (in Chinese),25(7):1529-1560.Yin A,Nie S Y.1993.An indentation model for the north and south China collision and the development of the Tan-Lu and Honam fault systems,eastern Asia.Tectonics,12(4):801-813.

Yin A.2006.Cenozoic tectonic evolution of the Himalayan orogen as constrained by along-strike variation of structural geometry,exhumation history,and foreland sedimentation.Earth-Science Reviews,76(1-2):1-131.

Zheng Y F,Ye K,Zhang L F.2009.Developing the plate tectonics from oceanic subduction to continental collision.Chinese Science Bulletin,54(15):2549-2555.

Zheng Y F.2012.Metamorphic chemical geodynamics in continental subduction zones.Chemical Geology,328:5-48.

附中文參考文獻(xiàn)

李忠海.2014.大陸俯沖-碰撞-折返的動(dòng)力學(xué)數(shù)值模擬研究綜述.中國(guó)科學(xué):地球科學(xué),44(5):817-841.

李忠海,許志琴.2015.大洋俯沖和大陸碰撞沿走向的轉(zhuǎn)換動(dòng)力學(xué)及流體-熔體活動(dòng)的作用.巖石學(xué)報(bào),31(12):3524-3530.

李忠海,劉明啟,Gerya T.2015.俯沖隧道中物質(zhì)運(yùn)移和流體-熔體活動(dòng)的動(dòng)力學(xué)數(shù)值模擬.中國(guó)科學(xué):地球科學(xué),45(7):881-899.

馬醒華,楊振宇.1993.中國(guó)三大地塊的碰撞拼合與古歐亞大陸的重建.地球物理學(xué)報(bào),36(4):476-488.

石耀霖,范桃園.2001.大洋巖石層拖曳窄條陸殼俯沖的極限尺度分析——以新西蘭南島和大別山超高壓變質(zhì)帶為例.地球物理學(xué)報(bào),44(6):754-760.

石耀霖,范桃園.2003.大洋板塊拖曳窄條陸殼俯沖——大規(guī)模超高壓變質(zhì)形成的構(gòu)造條件.地質(zhì)學(xué)報(bào),76(1):77-82.

楊經(jīng)綏,許志琴,張建新等.2009.中國(guó)主要高壓-超高壓變質(zhì)帶的大地構(gòu)造背景及俯沖/折返機(jī)制的探討.巖石學(xué)報(bào),25(7):1529-1560.

鄭永飛,葉凱,張立飛.2009.發(fā)展板塊構(gòu)造:從洋殼俯沖到大陸碰撞.科學(xué)通報(bào),54(13):1799-1803.

(本文編輯何燕)

Constraints of 3-D plate geometry on the dynamics of continental deep subduction

LI Zhong-Hai,SHI Yao-Lin

Key Laboratory of Computational Geodynamics,College of Earth Sciences,University of Chinese Academy of Sciences,Beijing 100049,China

Investigations of continental deep subduction and Ultra-High Pressure (UHP)metamorphism play significant roles in better understanding the continental geodynamics.Thereby,a number of geological and geophysical observations,as well as numerical simulations,have already been conducted systematically.However,along-strike variations are commonly observed in continental subduction/collision zones and mountain belts in nature.These three-dimensional (3-D)characteristics may be strongly dependent on the initial geometry and kinematics of converging continental plates.In order to study this problem,we have conducted 3-D high-resolution numerical models of subduction with rectangular or wedge-shaped continental plates,in which the subducting continental plate is laterally neighbored by oceanic subduction zone.The results demonstrate that oceanic subduction is continuous;however,slab break-off occurs in the continental subduction zone,after which the dip angle of the residue continental slab is decreased.In the model with rectangular continent,slab break-off occurs at around 150 km depth;whereas in the wedge-shaped continental subduction model,the depth of slab break-off is larger at about 250～300 km.The different phenomena are mainly attributed to the competing effects between the upward negative buoyancy of subducting continental slab and the downward drag from along-strike neighboring oceanic slab.In addition,after slab break-off,no matter in the model with a rectangular or wedge-shaped continent,a narrow continental marginal lithosphere (60～70 km in width)can keep going with the subducting oceanic slab to the mantle.This provides a good explanation for the continental sliver subduction dragged by neighboring oceanic subduction in New Zealand.The numerical models are also compared with the Himalaya and Qinling-Dabie-Sulu mountain belts,which provide significant implications for the along-strike variation in distributions of HP-UHP metamorphic rocks.

Continental subduction;UHP metamorphism;Plate geometry;Numerical modeling;Himalaya;Dabie-Sulu

李忠海，石耀霖.2016.三維板塊幾何形態(tài)對(duì)大陸深俯沖動(dòng)力學(xué)的制約.地球物理學(xué)報(bào)，59(8):2806-2817,

10.6038/cjg20160808.

Li Z H,Shi Y L.2016.Constraints of 3-D plate geometry on the dynamics of continental deep subduction.Chinese J.Geophys.(in Chinese),59(8):2806-2817,doi:10.6038/cjg20160808.

國(guó)家重點(diǎn)基礎(chǔ)研究發(fā)展計(jì)劃項(xiàng)目(2015CB856106),國(guó)家自然科學(xué)基金項(xiàng)目(41590860，41304071),中國(guó)科學(xué)院戰(zhàn)略性科技先導(dǎo)專項(xiàng)B(XDB18020104)和中組部青年千人計(jì)劃項(xiàng)目聯(lián)合資助.

李忠海，男，1982年生，博士，研究員，從事計(jì)算地球動(dòng)力學(xué)研究.E-mail：li.zhonghai@ucas.ac.cn

10.6038/cjg20160808

P541，P313

2015-08-02，2016-04-26收修定稿