皇甫鵬鵬, 王岳軍, 范蔚茗, 李忠海, 王喻鳴, 周永智

(1.中國科學(xué)院 廣州地球化學(xué)研究所 同位素地球化學(xué)國家重點實驗室, 廣東 廣州 510640; 2.中國科學(xué)院大學(xué) 計算地球動力學(xué)重點實驗室, 地球科學(xué)學(xué)院, 北京 100049; 3.中山大學(xué) 地球科學(xué)與地質(zhì)工程學(xué)院,廣東 廣州 510275; 4.中國科學(xué)院 青藏高原研究所 卓越研究中心, 北京 100101; 5.廣州市執(zhí)信中學(xué), 廣東廣州 510080)

大洋平板俯沖的數(shù)值模擬再現(xiàn): 洋–陸匯聚速率影響

皇甫鵬鵬1, 2, 3, 王岳軍3, 4*, 范蔚茗4, 李忠海2, 王喻鳴5, 周永智3

(1.中國科學(xué)院 廣州地球化學(xué)研究所 同位素地球化學(xué)國家重點實驗室, 廣東 廣州 510640; 2.中國科學(xué)院大學(xué) 計算地球動力學(xué)重點實驗室, 地球科學(xué)學(xué)院, 北京 100049; 3.中山大學(xué) 地球科學(xué)與地質(zhì)工程學(xué)院,廣東 廣州 510275; 4.中國科學(xué)院 青藏高原研究所 卓越研究中心, 北京 100101; 5.廣州市執(zhí)信中學(xué), 廣東廣州 510080)

利用地球動力學(xué)數(shù)值模擬方法探討了洋-陸匯聚時, 大洋巖石圈的絕對俯沖速率和上覆大陸巖石圈的向洋絕對逆沖速率對俯沖模式的影響, 尤其是上覆大陸的向洋絕對逆沖速率與平板俯沖之間的關(guān)系。模型結(jié)果顯示, 對于年齡為40 Ma的含正常洋殼厚度的大洋巖石圈, 在初始俯沖角度為現(xiàn)今洋–陸俯沖平均傾角的極小值(19°)條件下, 低速大洋俯沖(絕對俯沖速率≤3 cm/a)且上覆大陸巖石圈向洋絕對逆沖速率≥1 cm/a時, 具備形成平板俯沖的條件。當(dāng)中–高速大洋俯沖(絕對俯沖速度>3 cm/a)時, 在上覆大陸的絕對逆沖速率不小于俯沖速率時可以形成平板俯沖。當(dāng)增加初始俯沖角度到平均傾角的極大值(36°)時, 僅在低速大洋俯沖(絕對俯沖速率≤3 cm/a)且絕對逆沖速率達(dá)到 10 cm/a時(自然界中基本不存在), 才有可能出現(xiàn)平板俯沖, 其他情況均表現(xiàn)為陡俯沖。我們的模擬結(jié)果表明: (1)較高的大洋巖石圈絕對俯沖速率更容易克服板間耦合作用力而有利于陡俯沖形成; (2)較高的上覆大陸絕對逆沖速率更有利于俯沖板片彎曲而趨向于平板俯沖形成; (3)上覆大陸朝向海溝的逆沖速率會在俯沖板片下方產(chǎn)生水平向陸的地幔流, 絕對逆沖速率越大該地幔流越強烈,導(dǎo)致作用于板片下表面的水平剪切分量越大而有利于板片彎折和平板俯沖發(fā)生; (4)初始俯沖角度的增加對平板俯沖的形成起到強烈抑制作用。這些能被現(xiàn)今平板俯沖, 如具相似洋–陸匯聚速率條件的南美洲西海岸平板俯沖實例所驗證。

平板俯沖; 數(shù)值模擬; 俯沖速率; 大陸巖石圈向洋逆沖速率; 洋殼初始俯沖角度

0 引 言

根據(jù)洋–陸俯沖帶內(nèi)是否存在水平墊置于上覆大陸巖石圈底部的俯沖板片將大洋俯沖分為平俯沖(平板俯沖)和陡俯沖。但是現(xiàn)有資料顯示, 平板俯沖作為一種重要而又特殊的俯沖模式僅占現(xiàn)今全球洋–陸俯沖帶總長度的 10％左右(圖1, Gutscher et al., 2000a; Gutscher, 2002), 且集中發(fā)育于智利中部(Kay and Abbruzzi, 1996)、秘魯(Petford and Atherton, 1996)、厄瓜多爾(Beate et al., 2001)、哥斯達(dá)黎加(Grafe et al., 2002; Gardner et al., 2013)、墨西哥西南部(Suarez et al., 1990; Skinner and Clayton, 2011)、美國西海岸卡斯凱迪亞地區(qū)(Defant and Drummond, 1993)、阿拉斯加?xùn)|南部(Fuis et al., 2008)及日本南開海槽(Morris, 1995)。

與常見的大洋陡俯沖相比, 平板俯沖常伴隨一系列特殊地質(zhì)現(xiàn)象, 如上覆大陸巖石圈擠壓變形強烈、地殼增厚明顯, 從海溝向大陸內(nèi)部規(guī)律性出現(xiàn)大陸弧火山活動或間歇性巖漿活動(Gutscher et al., 2000a; Kay and Mpodozis, 2001; Sigloch et al., 2008),同時伴隨地表的動態(tài)抬升和俯沖帶上覆巖石圈的異常熱結(jié)構(gòu)(Dávila and Lithgow-Bertelloni, 2015)。正因如此, Gurscher et al. (2000a)用平板俯沖模型解釋了南美洲西海岸區(qū)域出現(xiàn)的與年輕大洋巖石圈俯沖相關(guān)的上新世–第四紀(jì)埃達(dá)克火山巖; Sigloch et al. (2008)提出包含長距離平板俯沖的雙階段俯沖模型解釋了北美西部 75～70 Ma的大范圍異常巖漿作用和長距離陸內(nèi)造山運動。

圖1　環(huán)太平洋俯沖帶大洋板片傾角示意圖(據(jù)Skinner and Clayton, 2011修改; 數(shù)據(jù)來自Lallemand et al., 2005; 圖中白色網(wǎng)格表示與平板俯沖區(qū)域相關(guān)聯(lián)洋脊)Fig.1 Map of circum-Pacific sunduction zone showing the dip angle of the subducting slabs

近 20年的研究對形成平板俯沖的影響因素取得了重要認(rèn)識, 認(rèn)為影響平板俯沖發(fā)生最重要的因素是俯沖板片的正浮力作用(Pilger, 1981; Brocher et al., 1994; Kay and Abbruzzi, 1996; Gutscher et al., 2000b; van Hunen et al., 2002b, 2004; Anderson et al., 2007; Espurt et al., 2008; Fuis et al., 2008; Sak et al., 2009; Skinner and Clayton, 2011), 同時上覆大陸巖石圈較冷的熱結(jié)構(gòu)和朝向海溝的絕對逆沖速率及板片吸力也有利于平板俯沖的發(fā)育(Stevenson and Turner, 1977; Tovish et al., 1978; Gutscher et al., 2000b; van Hunen et al., 2001, 2002a, 2002b, 2004; Perez-Gussinye et al., 2008; Li et al., 2011; Skinner and Clayton, 2011; Manea et al., 2012; Rodriguez-Gonzalez and Negredo, 2012)。盡管前人對平板俯沖形成機制研究較多, 但更多的研究仍集中于板片浮力因素。對于匯聚速率、尤其是上覆巖石圈向洋的絕對逆沖速率對俯沖模式的影響仍探究不多。

現(xiàn)今主要平板俯沖區(qū)域的統(tǒng)計數(shù)據(jù)(表 1)顯示,除了所有平板俯沖區(qū)域都與年輕大洋巖石圈或浮力異常區(qū)相關(guān)外, 另一個重要特征是過半的區(qū)域, 其上覆大陸巖石圈向洋絕對逆沖速率都大于大洋巖石圈絕對俯沖速率。例如, 智利中部、秘魯、厄瓜多爾和卡斯凱迪亞地區(qū)的大洋巖石圈絕對俯沖速率分別為3.0 cm/a, 2.5 cm/a, 1.6 cm/a和0.8 cm/a, 而上覆大陸朝向海溝的絕對逆沖速率分別達(dá)到 4.5 cm/a, 4.5 cm/a, 3.8 cm/a和2.4 cm/a。結(jié)合前人研究成果可知, 俯沖區(qū)域的速度邊界條件是控制大洋俯沖模式的重要影響因素之一(Uyeda and Kanamori, 1979; Jarrard, 1986; Lallemand et al., 2005)。那么, 匯聚速率, 尤其是上覆大陸巖石圈絕對逆沖速率是如何影響洋–陸俯沖模式的？又是如何制約平板俯沖的形成？能否再現(xiàn)現(xiàn)今平板俯沖型式？這些問題對理解平板俯沖發(fā)生的動力條件尤為重要。為此, 本文利用二維地球動力學(xué)數(shù)值模擬方法研究了洋–陸匯聚時大洋巖石圈的絕對俯沖速率和上覆大陸巖石圈的向洋絕對逆沖速率對俯沖模式的影響效應(yīng), 限定了絕對逆沖速率的大小與平板俯沖形成之間的耦合關(guān)系, 并結(jié)合現(xiàn)今平板俯沖實例對模擬結(jié)果進(jìn)行了討論。

表1　現(xiàn)今平板俯沖區(qū)域的主要地球動力學(xué)參數(shù)Table 1 Main geodynamic conditions of the present-day flat subduction

1 數(shù)值模擬方法介紹

1.1 控制方程

本文模型計算采用二維地球動力學(xué)數(shù)值算法“I2VIS”(Gerya and Yuen, 2003a)。根據(jù)地球動力學(xué)特點, 一般主要針對三組控制方程進(jìn)行求解計算, 包括斯托克斯流體動力學(xué)方程、物質(zhì)守恒方程以及熱量守恒方程(Li, 2014)。

(1) 斯托克斯方程:

其中, g是重力加速度; 密度ρ依賴于溫度T、壓力P、巖石類型C和部分熔融比例M;xxσ＇,zzσ＇和xzσ＇是偏應(yīng)力張量。

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

其中,xν和zν表示水平速度和垂向速度。

(3) 熱量守恒方程:

其中, Cp是等壓熱容;表示溫度對時間的物質(zhì)導(dǎo)數(shù); qx和qz代表水平方向和垂直方向熱流值; k受熱傳導(dǎo)率, 受溫度T、壓力P和巖石類型C影響(Clauser and Huenges, 1995); Hr表示放射性生熱,是僅與巖石類型相關(guān)的常量; Ha表示絕熱變壓生熱, 與壓力變化相關(guān); Hs表示剪切生熱; ε.xx, ε.zz和ε.xz是應(yīng)變率張量。

1.2 粘–塑性流變學(xué)性質(zhì)

偏應(yīng)力張量和應(yīng)變率張量采用粘–塑性本構(gòu)方程, 對于不可壓縮的粘滯性變形可表示為:

式中,effη 為有效粘滯系數(shù), 是關(guān)于溫度、壓力、物質(zhì)成分、應(yīng)變率和熔融程度的函數(shù)。

韌性流變的粘滯系數(shù)與溫度、壓力、物質(zhì)成分、應(yīng)變率和熔融程度相關(guān), 可表示為:

實際巖體的粘–塑性流變是韌性流變和脆性/塑性流變兩者的結(jié)合。因此, 我們在數(shù)值算法中采用Mohr-Coulomb屈服準(zhǔn)則(Ranalli, 1995):

式中,yieldσ為屈服應(yīng)力;IIε.表示應(yīng)變率張量二階不變量;0C是 P=0條件下的巖石內(nèi)聚力; φ是內(nèi)摩擦角; λ是孔隙流體系數(shù)(Brace and Kohlstedt, 1980);表示有效內(nèi)摩擦角, 與內(nèi)摩擦角和孔隙流體系數(shù)相關(guān)。

1.3 部分熔融

基于實驗巖石學(xué)的約束條件, 該算法考慮了地殼巖石部分熔融行為(Gerya and Yuen, 2003b; Burg and Gerya, 2005)。并近似認(rèn)為部分熔融體積比例與溫度間存在如下線性關(guān)系:

部分熔融巖石的有效密度取決于熔融比例:

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

2 模型設(shè)計

根據(jù)現(xiàn)今全球洋–陸俯沖帶剖面的統(tǒng)計結(jié)果,本文設(shè)計初始俯沖角度為19°和36°的兩組模擬實驗,分別對應(yīng)于俯沖大洋板片平均傾角的極大值和極小值(Lallemand et al., 2005), 在上述兩組模型基礎(chǔ)上系統(tǒng)改變大洋巖石圈的絕對俯沖速率和上覆大陸巖石圈的向洋絕對逆沖速率, 以此探究上覆大陸巖石圈的向洋絕對逆沖速率對洋–陸俯沖模式的制約,以及初始俯沖角度對平板俯沖形成的影響。

設(shè)計初始模型長度為4000 km, 深度為670 km,使用不規(guī)則的有限差分單元對模型進(jìn)行離散化, 其中俯沖帶區(qū)域采用2 km×2 km單元, 而在模型邊界區(qū)域采用30 km×30 km 單元。模型設(shè)計大洋板塊和大陸巖石圈板塊的初始長度均為2000 km, 如圖2所示。其中, 大洋巖石圈年齡40 Ma, 厚度根據(jù)半空間冷卻模型(Turcotte and Schubert, 2002)計算設(shè)計為74 km,其中洋殼厚度8 km, 大陸巖石圈厚度120 km(包括35 km的大陸地殼和85 km的巖石圈地幔, 而大陸地殼部分則設(shè)計為6 km大陸沉積物地殼、14 km上地殼和15 km下地殼)。不同單元內(nèi)的巖石類型及其屬性詳見表 2和表 3。模型底部邊界設(shè)計為滲透性邊界, 采用近無限深度的外部自由滑動邊界條件(Burg and Gerya, 2005; Ueda et al., 2008; Li et al., 2010)以滿足在計算模型區(qū)域下方虛擬邊界上(例如 100 km 處)的自由滑動條件, 其他邊界則采用自由滑動的速度邊界條件。與普通的自由滑動邊界條件一樣, 外部自由滑動邊界條件同樣需要滿足計算區(qū)域內(nèi)的物質(zhì)守恒。模型計算過程中的洋–陸匯聚通過施加在大洋巖石圈朝向海溝的俯沖速率和大陸巖石圈朝向海溝的逆沖速率來共同實現(xiàn), 因為兩者均為模型內(nèi)施加的速度邊界條件, 并且在模型計算過程中保持恒定, 因此可以認(rèn)為是相當(dāng)于靜止參照物的絕對運動速率。根據(jù)統(tǒng)計資料(Lallemand et al., 2005), 大于10 cm/a的大陸巖石圈向洋絕對逆沖速率在現(xiàn)今全球俯沖帶處幾乎是不存在的, 因而在模型實驗設(shè)計中不予考慮, 從而限定其變化范圍為0～10 cm/a。模型中設(shè)計的大洋板片絕對俯沖速率和大陸巖石圈絕對逆沖速率見表4和表5。另外, 模型中的初始俯沖角度由薄弱帶角度形式來加以實現(xiàn)。

(a) 模型主要區(qū)域(1500 km×670 km)的初始物質(zhì)與溫度場。圖中白色線條是等溫線, 黃色字體和圖示代表模型邊界條件; (b) 初始俯沖帶特征。模型中的顏色代表不同巖石組成, 如圖例所示。圖例中部分熔融巖石在初始模型中并未顯示, 但在模擬運算過程中由于溫度增加至巖石熔點時而逐漸出現(xiàn)。圖2　初始模型和邊界條件Fig.2 Initial model configuration and boundary conditions

表2　數(shù)值模擬采用的粘滯性流變參數(shù)Table 2 Parameters of viscous flow laws in the numerical experiments

關(guān)于模型的初始熱邊界條件, 設(shè)計模型頂部為固定溫度(0 ℃), 大洋和大陸巖石圈底邊界溫度為1300 ℃ (Turcotte and Schubert, 2002), 軟流圈地幔的溫度梯度為 0.5 ℃/km, 兩側(cè)邊界的水平方向溫度梯度為零(即零熱流)。底部邊界采用的是外部邊界固定溫度條件, 即在模型底邊界之下1000 km處假設(shè)一個固定的地幔溫度(根據(jù)地幔固定溫度梯度計算得出), 這樣就可以使得在 670 km處的底部滲透邊界上的溫度和熱流可以隨模型計算而動態(tài)調(diào)整。

同時, 為了更真實反映俯沖碰撞的地表變化,在地殼表面和模型頂界面之間設(shè)計了一層相對高粘滯度的偽空氣層, 其與上地殼的接觸面被用來模擬模型計算過程中的地貌起伏變化, 且該地貌起伏面模擬包含近似的地表剝蝕和沉積作用(Gerya and Yuen, 2003a; Burg and Gerya, 2005)。偽空氣層的初始厚度在大陸巖石圈之上為 10 km, 在大洋地殼之上為12 km。

表3　數(shù)值模型中的主要材料參數(shù)Table 3 Parameters of the materials in the numerical models

表4　初始俯沖角度為19°的模型實驗結(jié)果Table 4 Numerical results at the initial subduction angle of 19°

3 模擬結(jié)果

首先根據(jù)初始俯沖角度的不同, 將模型實驗分為低初始俯沖角度(19°)和高初始俯沖角度(36°)兩組,模擬結(jié)果見表 4和表 5。依大洋巖石圈的絕對俯沖速率(subυ), 將歸類為低速大洋俯沖,歸類為中速大洋俯沖,歸類為高速大洋俯沖。

3.1 低初始俯沖角度(19°)

現(xiàn)今全球洋–陸俯沖帶的統(tǒng)計資料顯示, 在深度<125 km的大洋巖石圈淺部, 其俯沖角度平均值為 27.5°±8.5°(Lallemand et al., 2005), 本組模型實驗的初始俯沖角度設(shè)計為這一平均值的極小值 19°來代表低的初始俯沖角度情況, 以此探究其在低速、中速和高速洋–陸俯沖條件下的俯沖型式。

表5　初始俯沖角度為36°的模型實驗結(jié)果Table 5 Numerical results at the initial subduction angle of 36°

模型結(jié)果顯示(表 4), 低速大洋俯沖條件下, 只要上覆大陸巖石圈存在朝向海溝的絕對逆沖速率(≥1 cm/a), 就能導(dǎo)致平板俯沖的發(fā)生。反之, 若無絕對逆沖速率, 則為陡俯沖。圖3和圖4分別展示俯沖速率為2 cm/a, 絕對逆沖速率為0和1 cm/a兩個模型的計算結(jié)果, 前者為陡俯沖, 后者為平板俯沖。

以絕對俯沖速率為2 cm/a和絕對逆沖速度為1 cm/a的模型為例, 論述低速大洋俯沖條件下平板俯沖的演化過程(如圖 4所示)。洋-陸俯沖開始后,大洋巖石圈沿著初始俯沖界面進(jìn)入軟流圈, 由于較小初始俯沖角度導(dǎo)致的較大板間耦合力、較小的大洋俯沖速率不利于上下板塊的解耦、以及上覆大陸巖石圈朝向海溝的逆沖速率更有利于俯沖板片的幾何彎折的三種條件下, 大洋板片粘附在上覆大陸巖石圈地幔底部以近水平角度向前運動(圖4a)。另外, 俯沖帶處的擠壓應(yīng)力使上覆大陸巖石圈前緣發(fā)生彎折, 隨后遭強烈剝蝕的大陸沉積物和地殼物質(zhì)沉積于海溝處并進(jìn)入俯沖隧道內(nèi), 形成“俯沖侵蝕”現(xiàn)象。模型運行到40.4 Ma時, 大洋俯沖平板的長度能達(dá)到600 km (圖4b), 俯沖板片前端近乎垂直進(jìn)入軟流圈地幔。此時大洋沉積物隨著俯沖板片近水平進(jìn)入到俯沖板片與上覆巖石圈地幔之間的俯沖通道內(nèi), 并伴隨的俯沖沉積物部分熔融。同時在俯沖通道內(nèi)的沉積物和周圍物質(zhì)間的巨大密度差所產(chǎn)生的正浮力作用促使部分大洋沉積物沿俯沖通道折返至地表, 并伴隨俯沖洋殼物質(zhì)的減壓熔融。隨著俯沖作用的持續(xù)進(jìn)行(圖 4c), 在俯沖板片的拉力作用下發(fā)生回撤而導(dǎo)致板片上方地幔楔內(nèi)速度場產(chǎn)生強烈擾動。

當(dāng)絕對俯沖速率和絕對逆沖速度均為5 cm/a時,中速大洋俯沖條件下平板俯沖的演化過程如圖6a所示。洋–陸俯沖開始后, 在較大的俯沖速率下, 大洋板片并不像低速大洋俯沖條件下所表現(xiàn)出的粘附于上覆大陸巖石圈地幔底部的水平向前運動, 而是以大約30°的俯沖角度斜向進(jìn)入軟流圈(圖6a)。大洋和大陸巖石圈的相向匯聚導(dǎo)致俯沖帶處的擠壓應(yīng)力促使上覆大陸巖石圈前緣向上彎折并遭受強烈剝蝕,此時部分大陸沉積物和地殼物質(zhì)沉積于海溝處并進(jìn)入俯沖隧道。隨著大洋俯沖的持續(xù)進(jìn)行(圖6b), 上覆大陸巖石圈在大洋板片之上產(chǎn)生巨大壓力, 促使板片在上覆巖石圈底部向上彎折而引起板片傾角陡然變小, 誘發(fā)平板俯沖發(fā)生。此時部分大洋沉積物從俯沖板片上剝離而堆積于海溝處形成增生楔, 且伴隨著部分沉積物的俯沖折返。大洋俯沖19 Ma后(圖6c), 500 km 長的大洋平板俯沖正式形成, 但在俯沖平板和上覆大陸巖石圈地幔之間仍殘余有部分軟流圈物質(zhì), 該軟流圈物質(zhì)將在隨后大洋俯沖平板的持續(xù)向上運動直至墊付于上覆巖石圈底部而被逐漸擠出。之前形成的增生楔進(jìn)入俯沖隧道而形成“俯沖侵蝕”現(xiàn)象, 并隨著溫度的上升而發(fā)生部分熔融。

圖3　初始俯沖角度為19°, 板片絕對俯沖速率為2 cm/a, 上覆大陸巖石圈絕對逆沖速率為0時運行到44 Ma的模擬結(jié)果Fig.3 Numerical result at 44 Ma when assuming that the initial subduction angle is 19°, slab absolute subduction velocity is 2 cm/a and continental absolute trenchward velocity is 0 cm/a

圖4　初始俯沖角度為19°, 板片絕對俯沖速率為2 cm/a, 上覆大陸巖石圈絕對逆沖速率為1 cm/a的模型分別在16.7 Ma (a), 40.4 Ma (b)和57.6 Ma (c)時的演化結(jié)果(結(jié)果顯示模型從開始階段即為平板俯沖)Fig.4 Numerical results at 16.7 Ma (a), 40.4 Ma (b), and 57.6 Ma (c) when assuming that the initial subduction angle is 19°, slab absolute subduction velocity is 2 cm/a and continental absolute trenchward velocity is 1 cm/a

圖5　初始俯沖角度為19°, 板片絕對俯沖速率為5 cm/a, 上覆大陸巖石圈絕對逆沖速率為4 cm/a的模型在運行到19.1 Ma的結(jié)果Fig.5 Numerical result at 19.1 Ma when assuming that initial subduction angle is 19°, slab absolute subduction velocity is 5 cm/a and continental absolute trenchward velocity is 4 cm/a

圖6　初始俯沖角度為19°, 板片絕對俯沖速率為5 cm/a, 大陸巖石圈絕對逆沖速率為5 cm/a的模型分別在 5.9 Ma (a), 13.9 Ma (b)和19.1 Ma (c)時的演化結(jié)果, 模型逐漸從陡俯沖轉(zhuǎn)換為平板俯沖Fig.6 Numerical results at 5.9 Ma (a), 13.9 Ma (b) and 19.1 Ma (c) when assuming that initial subduction angle is 19°, slab absolute subduction velocity is 5 cm/a and continental absolute trenchward velocity is 5 cm/a

圖7　初始俯沖角度為19°, 板片絕對俯沖速率為9 cm/a, 上覆大陸巖石圈絕對逆沖速率為7 cm/a的模型在運行到11.4 Ma的結(jié)果Fig.7 Numerical result at 11.4 Ma when assuming that initial subduction angle is 19°, slab absolute subduction velocity is 9 cm/a and continental absolute trenchward velocity is 7 cm/a

我們以絕對俯沖速率為9 cm/a和絕對逆沖速度為8 cm/a的模型為例, 闡述平板俯沖的形成演化過程。隨著大洋俯沖開始(圖 8a), 在應(yīng)力作用下, 上覆大陸巖石圈在俯沖帶處向上彎折變形并遭受強烈剝蝕, 部分侵蝕的大陸物質(zhì)沉積在海溝處并隨著大洋巖石圈俯沖進(jìn)入俯沖隧道, 并以約 30°傾角進(jìn)入軟流圈。隨后大洋板片逐漸在大陸巖石圈底部深度處向斜下方彎折(圖 8b)而形成平板俯沖型式, 隨后大洋板片近水平墊付于上覆大陸巖石圈地幔底部(圖8c)。

3.2 高初始俯沖角度(36°)

根據(jù)現(xiàn)今全球洋–陸俯沖帶的統(tǒng)計資料, 36°為大洋巖石圈淺部俯沖角度平均值的極大值(Lallemand et al., 2005), 因而可代表初始俯沖角度較大的一類洋–陸俯沖。該組模型實驗相對于前述模型僅增大初始俯沖角度, 其余參數(shù)均保持不變。計算結(jié)果顯示(表5), 增大初始俯沖角度會極大地限制平板俯沖的形成。此時只有在絕對俯沖速率≤3 cm/a、且絕對逆沖速率達(dá)到10 cm/a時才能產(chǎn)生平板俯沖。

圖9展示了絕對俯沖速率為3 cm/a且絕對逆沖速率依次為6 cm/a、8 cm/a和10 cm/a經(jīng)過15 Ma時的模擬結(jié)果。模擬結(jié)果顯示, 前兩個模型中俯沖板片分別以40°和25°的角度傾斜進(jìn)入軟流圈(圖9a、b), 且后者相對于前者具有更小的俯沖角度, 但由于上覆大陸巖石圈的絕對逆沖速率所產(chǎn)生的壓力還不足以使俯沖板片發(fā)生彎折仍表現(xiàn)出陡俯沖型式。當(dāng)大陸巖石圈絕對逆沖速率增大到 10 cm/a時, 在上覆大陸巖石圈的擠壓應(yīng)力作用下, 俯沖大洋板片于 120 km深度處發(fā)生幾何彎折而導(dǎo)致平板俯沖產(chǎn)生(圖 9c), 此時俯沖平板段的長度 300 km, 但其俯沖板片前鋒仍以約20°角度傾斜俯沖。上述三個模型的模擬結(jié)果都顯示在俯沖期間存在“俯沖侵蝕”現(xiàn)象導(dǎo)致上覆大陸巖石圈物質(zhì)進(jìn)入俯沖帶內(nèi)。同時隨大洋巖石圈俯沖進(jìn)入俯沖隧道的大洋沉積物由于熱傳導(dǎo)作用增溫弱化而部分熔融。

4 討 論

4.1 大洋巖石圈絕對俯沖速率對俯沖模式的影響

在大洋俯沖動力學(xué)研究中, 能否形成平板俯沖主要受板片正浮力作用、板間耦合力作用以及外部運動學(xué)條件等因素的影響和制約(Gutscher et al., 2000b; van Hunen et al., 2002b, 2004; Anderson et al., 2007; Perez-Gussinye et al., 2008; Rodriguez-Gonzalez and Negredo, 2012)。在其他動力學(xué)參數(shù)不變的情況下, 大洋板片絕對俯沖速率的大小在很大程度上決定了大洋板片能否克服初始板間耦合力作用。如果其絕對俯沖速率較小, 板間初始耦合力作用占據(jù)主導(dǎo)地位, 則大洋板片會持續(xù)粘附于上覆大陸巖石圈底部, 從而為平板俯沖產(chǎn)生創(chuàng)造有利條件。當(dāng)大洋俯沖速率增大到足以克服板塊間耦合力作用時, 則大洋板片會與上覆大陸巖石圈發(fā)生解耦而進(jìn)入軟流圈形成陡俯沖。因此, 大洋巖石圈絕對俯沖速率在一定程度上制約著洋–陸俯沖模式的轉(zhuǎn)換和演化。

表 4的模擬結(jié)果表明, 在初始俯沖角度相對較低情況下(如19°), 上覆大陸巖石圈向洋絕對逆沖速率保持不變的條件下, 大洋巖石圈的絕對俯沖速率越大則越不利于平板俯沖的形成, 在當(dāng)大洋巖石圈的絕對俯沖速率大于上覆大陸巖石圈向洋絕對逆沖速率時則只能導(dǎo)致陡俯沖型式。原因如下: (1)絕對俯沖速率越大越利于俯沖板片與上覆大陸巖石圈間的解耦, 因而越不利于平板俯沖產(chǎn)生; (2)相對較小的大洋俯沖速率能使俯沖板片滯留于軟流圈地幔的時間較長而有利于俯沖板片增溫弱化, 致使其在同樣的上覆大陸巖石圈壓力作用下更容易發(fā)生幾何彎折而利于平板俯沖產(chǎn)生; (3)更重要的是, 上覆巖石圈向洋逆沖速率所引發(fā)的板片下方的向陸方向的局部地幔流的強度與平板俯沖形成關(guān)系密切, 如果該地幔流占據(jù)主導(dǎo)地位則有利于平板俯沖的產(chǎn)生, 反之則不利于。如圖10所示, 在初始俯沖角度為19°、上覆大陸巖石圈逆沖速率恒定為5 cm/a時, 大洋板片俯沖速率為4 cm/a或5 cm/a時, 由上覆板塊引發(fā)的板片下方的向陸地幔流占據(jù)主導(dǎo)而有利于平板俯沖的發(fā)育; 反之, 如果板片俯沖速率增大至 6 cm/a 或7 cm/a時, 則由板片俯沖引發(fā)的順板片向下的地幔流占據(jù)主導(dǎo)地位而有利于陡俯沖的發(fā)生。

4.2 大陸巖石圈絕對逆沖速率對平板俯沖的制約

有利于大洋板片發(fā)生彎折的外部運動學(xué)條件是導(dǎo)致平板俯沖發(fā)生的重要因素之一, 而上覆大陸巖石圈的向洋絕對逆沖速率則是控制大洋俯沖動力學(xué)的重要外部運動學(xué)參數(shù)之一。因此, 上覆大陸絕對逆沖速率的大小對平板俯沖形成與否至關(guān)重要。

圖8　初始俯沖角度為19°, 板片絕對俯沖速率為9 cm/a, 上覆大陸巖石圈絕對逆沖速率為8 cm/a的模型分別在5.0 Ma (a), 9.0 Ma (b)和11.4 Ma (c)時的演化結(jié)果(結(jié)果顯示模型逐漸從陡俯沖轉(zhuǎn)換為平板俯沖)Fig.8 Numerical results at 5.0 Ma (a), 9.0 Ma (b) and 11.4 Ma (c) when assuming that initial subduction angle is 19°, slab absolute subduction velocity is 9 cm/a and continental absolute trenchward velocity is 8 cm/a

圖9　初始俯沖角度為36°, 板片絕對俯沖速率為3 cm/a的模型在大陸巖石圈絕對逆沖速率分別為6 cm/a (a), 8 cm/a (b)和10 cm/a (c)運算到15 Ma時的結(jié)果(顯示隨上覆巖石圈絕對逆沖速率增加, 大洋俯沖逐漸從陡俯沖轉(zhuǎn)變?yōu)槠桨甯_)Fig.9 Numerical results for continental absolute trenchward velocity ranging from 6 cm/a (a) to 8 cm/a (b) and then to 10 cm/a (c), when assuming that initial subduction angle is 36° and slab absolute subduction velocity is 3 cm/a

表4的模擬結(jié)果顯示, 在恒定大洋俯沖速率時,上覆大陸巖石圈的向洋絕對逆沖速率不小于、且與大洋板片俯沖速率差值越大則越有利于平板俯沖的形成。主要原因有: (1)上覆大陸巖石圈的絕對逆沖速率越大, 則直接作用于俯沖板片朝向斜下方的應(yīng)力越大, 使之越有利于板片在上覆大陸巖石圈底部強制性彎折而導(dǎo)致俯沖板片傾角減小; (2)上覆巖石圈向洋逆沖速率會導(dǎo)致俯沖板片下方產(chǎn)生水平向的向陸的擾動地幔流。逆沖速率越大則局部地幔流擾動越強烈, 作用于板片的水平剪切分量越大而越有利于板片彎折和平板俯沖發(fā)生。圖11所示, 在初始俯沖角度為36°、大洋板片絕對俯沖速率恒定為3 cm/a時,上覆巖石圈的向洋絕對逆沖速率為6 cm/a和8 cm/a時,俯沖板片的擾動地幔流指向深部而利于陡俯沖; 而當(dāng)逆沖速率增至 10 cm/a時, 則由上覆板塊向洋運動所引起的指向內(nèi)陸一側(cè)的地幔流占據(jù)主導(dǎo)而誘使平板俯沖發(fā)生。

在如上所述之平板俯沖能夠發(fā)生的特定情況下,上覆大陸巖石圈的絕對逆沖速率大小一定程度上也影響著平板俯沖平板段的發(fā)育長度。在恒定大洋板片絕對俯沖速率條件下, 較大的大陸巖石圈絕對逆沖速率可有利于平板俯沖到更長距離。圖12展示了初始俯沖角度為19°、低速大洋板片俯沖(絕對俯沖速率為3 cm/a)情況下, 大陸巖石圈絕對逆沖速率從1 cm/a增至6 cm/a時平板俯沖長度的模擬結(jié)果。該模擬結(jié)果顯示上覆大陸巖石圈絕對逆沖速率大小與俯沖平板長度間具有正相關(guān)性。模型結(jié)果也表明, 由于大洋板片負(fù)浮力作用, 即使大陸巖石圈絕對逆沖速率不斷增加, 但俯沖平板長度僅能保持在 1100 km以內(nèi)。

通過進(jìn)一步分析, 本文模型實驗結(jié)果也表明上覆板塊的絕對逆沖速率和俯沖板塊的絕對俯沖速率之間的絕對速率比與俯沖模式之間也存在一定的聯(lián)系。當(dāng)初始俯沖角度較小時(～19°), 無論大洋巖石圈的俯沖速率為多少, 只要此絕對速率比不小于 1.0均能導(dǎo)致平板俯沖的發(fā)生; 對于低速大洋俯沖情況,絕對速率比甚至小于 1.0也能發(fā)生, 因此我們可以把 1.0作為絕對速率比能否導(dǎo)致平板俯沖發(fā)生的閾值的下限。但需要注意的是, 該結(jié)論僅適用于大洋板塊初始俯沖角度較低的情況, 而對于初始俯沖角度較高的情況, 平板俯沖本身就極難發(fā)育, 因而也無此規(guī)律可循。

4.3 初始俯沖角度對洋–陸俯沖模式的影響

初始俯沖角度決定了上覆大陸巖石圈和俯沖板塊之間的初始接觸面積, 進(jìn)而直接約束板間耦合力大小。而板間耦合力, 在微觀上指的是兩板塊間可以積累和承受的物質(zhì)顆粒間的粘結(jié)力, 宏觀表現(xiàn)為兩板塊發(fā)生相對運動而不使板塊分離的摩擦力。當(dāng)板間耦合力相對較大時, 俯沖板片更容易粘附在上覆巖石圈底部向前水平運動一段距離后再傾斜進(jìn)入軟流圈地幔。如果被動大陸邊緣是緩長型的, 那么初始俯沖時的俯沖角度可能較小; 反之, 如果被動大陸邊緣是陡窄型的, 那么初始俯沖角度也就較大。

表4和5的模擬結(jié)果清晰展示了初始俯沖角度對俯沖型式的影響, 較大的初始俯沖角度會極大地抑制平板俯沖的發(fā)生。當(dāng)初始俯沖角度較小時: 低速大洋俯沖時, 只要存在大陸巖石圈絕對逆沖速率即可以產(chǎn)生平板俯沖;而在中高速大洋俯沖時, 只有在大陸巖石圈絕對逆沖速率超過大洋板片絕對俯沖速率才可導(dǎo)致平板俯沖發(fā)生, 且平板俯沖長度最長不會超過1100 km。在初始俯沖角度較高時, 平板俯沖多難孕育, 如初始俯沖角度增至36°, 則平板俯沖僅在大陸巖石圈絕對逆沖速率增至10 cm/a才有可能發(fā)生。

4.4 現(xiàn)今平板俯沖模式的實例驗證

現(xiàn)有資料顯示, 平板俯沖區(qū)域只占現(xiàn)今全球俯沖帶區(qū)域總長度的 10 ％左右(如 Gutscher et al., 2000a; Gutscher et al., 2000b; van Hunen et al., 2002b; Skinner and Clayton., 2011)。表1總結(jié)了主要平板俯沖區(qū)域的重要地球動力學(xué)參數(shù), 從中可以看出, 現(xiàn)今大部分平板俯沖區(qū)域的一個重要特征是: 與大洋板片絕對俯沖速率相比, 上覆大陸巖石圈一般都具有相對較大的絕對逆沖速率。

南美洲西海岸是現(xiàn)今全球平板俯沖最集中分布的區(qū)域(占全球平板俯沖區(qū)域總長度 50％以上), 也是平板俯沖研究最深入的地區(qū)。地震學(xué)證據(jù)(Gutscher, 2002; Hayes et al., 2012)顯示, 納茲板塊向南美洲大陸俯沖時, 在智力中部、秘魯以及厄瓜多爾等區(qū)域下方俯沖的大洋板片幾乎水平底墊于上覆南美大陸巖石圈之下, 形成向陸內(nèi)延伸幾百公里的低角度近水平俯沖。同樣, GPS觀測數(shù)據(jù)(DeMets et al., 1994)顯示, 南美大陸巖石圈在大西洋中脊向西的推力和納茲板塊偏東北向的俯沖與回撤作用下高速運動(～5 cm/a), 導(dǎo)致南美西海岸上覆大陸巖石圈向洋絕對逆沖速率普遍大于納茲板塊的絕對俯沖速率。例如, 在平板俯沖發(fā)育的智利中部、秘魯和厄瓜多爾地區(qū)的大洋板片絕對俯沖速率分別為3.0 cm/a, 2.5 cm/a和1.6 cm/a, 而上覆大陸巖石圈向洋絕對逆沖速率分別達(dá)到4.5 cm/a, 4.5 cm/a和3.8 cm/a。北美洲西海岸的卡斯凱迪亞地區(qū)發(fā)育的平板俯沖則是由Farallon板塊的殘余板片俯沖(Juan de Fuca)所致, 其觀測數(shù)據(jù)(DeMets et al., 1994; Gripp and Gordon, 2002)表明, Juan de Fuca板片以北北東方向斜向俯沖入北美大陸, 其垂直于俯沖帶的絕對俯沖速率較小,而其上覆北美洲大陸巖石圈則以>2 cm/a的速度向西逆沖, 其上覆大陸巖石圈較大洋板片具有更高的絕對速率。如此特征與本文數(shù)值模擬結(jié)論一致。

圖10　初始俯沖角度為19°, 上覆大陸巖石圈絕對逆沖速率為5 cm/a時在板片絕對俯沖速率分別為4 cm/a (a), 5 cm/a (b), 6 cm/a (c)和7 cm/a (d)俯沖1000 km時的結(jié)果及相應(yīng)的粘滯度場和速度場(結(jié)果顯示隨著大洋板片俯沖速率的增加, 大洋俯沖逐漸從平板俯沖轉(zhuǎn)變?yōu)槎父_)Fig.10 Numerical results for slab absolute subduction velocity ranging from 4 cm/a (a), 5 cm/a (b), 6 cm/a (c), to7 cm/a (d), and the corresponding viscosity and calculated velocity fields (a′), (b′), (c′), (d′)respectively, when assuming that initial subduction angle is 19°, continental absolute trenchward velocity is 5 cm/a, the simulations stop when convergence reaches 1000 km

圖11　初始俯沖角度為36°, 板片絕對俯沖速率為3 cm/a的模型在上覆大陸巖石圈絕對逆沖速率分別為6 cm/a (a), 8 cm/a (b) 和10 cm/a (c)的15 Ma時的粘滯度和速度場結(jié)果(結(jié)果顯示上覆大陸巖石圈絕對逆沖速率的增加可促進(jìn)大洋俯沖由陡俯沖轉(zhuǎn)變?yōu)槠桨甯_)Fig.11 Numerical results of viscosity and calculated velocity fields at the conditions of an initial subduction angle of 36°, slab absolute subduction velocity of 3 cm/a and continental absolute overthrusting velocity increasing from 6 cm/a (a) to 8 cm/a (b) and then to 10 cm/a (c), respectively, and the simulations stop at ～15 Ma

圖 12 平板俯沖長度與上覆大陸巖石圈向洋絕對逆沖速率關(guān)系圖(假設(shè)初始俯沖角度為 19°, 板片絕對俯沖速率 3 cm/a, 大陸巖石圈的絕對逆沖速率從 1 cm/a增至6 cm/a的模擬結(jié)果)Fig.12 Correlation between the absolute overthrusting velocity of continental lithosphere and the slab length of the flat subduction

同時, 如表 1所示, 現(xiàn)今平板俯沖幾乎均發(fā)育于大洋板片較為年輕的太平洋東海岸匯聚邊界, 所有平板俯沖區(qū)域的大洋板片年齡均小于 45 Ma, 且絕大多數(shù)小于 30 Ma, 這表明平板俯沖發(fā)生的重要條件是俯沖大洋巖石圈年齡遠(yuǎn)小于正常大洋巖石圈平均年齡 (60 Ma, Condie, 1997)。大洋俯沖板片正浮力增加(如年輕的大洋板片或增厚洋殼等)可能是導(dǎo)致平板俯沖發(fā)生的主要原因, 而上覆板塊朝向海溝的絕對逆沖速率同樣也制約著平板俯沖形成。

5 結(jié) 論

本文通過數(shù)值模擬研究再現(xiàn)了洋–陸匯聚速率對平板俯沖形成演化的影響效應(yīng), 獲得了如下認(rèn)識:

(1) 對年齡為 40 Ma的大洋板片, 當(dāng)初始俯沖角度為洋–陸匯聚時平均俯沖角度的極小值(19°)時,在低速大洋俯沖(≤3 cm/a)條件下, 較小的絕對逆沖速率(≥1 cm/a)即可導(dǎo)致平板俯沖; 而在中高速大洋俯沖(>3 cm/a)條件下, 平板俯沖的發(fā)生要求大陸巖石圈絕對逆沖速率不小于大洋板片絕對俯沖速率。當(dāng)初始俯沖角度增至洋-陸匯聚時平均俯沖角度的極大值(36°)時, 平板俯沖只有在大洋低速俯沖, 且大陸巖石圈絕對逆沖速率達(dá)到 10 cm/a以上時才可發(fā)生。初始俯沖角度的增加對平板俯沖的形成起到重要抑制作用。

(2) 較大的大洋板片絕對俯沖速率更易于克服板間的耦合作用力而不利于平板俯沖形成。相反,較大的大陸巖石圈絕對逆沖速率更易于促進(jìn)俯沖板片彎折而孕育平板俯沖型式。

致謝： 衷心感謝審稿專家中國科學(xué)技術(shù)大學(xué)冷偉老師和中國科學(xué)院地質(zhì)與地球物理研究所陳林老師對文章修改所提出的寶貴意見和建議！

Anderson M, Alvarado P, Zandt G and Beck S. 2007. Geometry and brittle deformation of the subducting Nazca Plate, Central Chile and Argentina. Geophysical Journal International, 171: 419–434.

Beate B, Monzier M, Spikings R, Cotten J, Silva J, Bourdon E and Eissen J P. 2001. Mio-Pliocene adakite generation related to flat subduction in southern Ecuador: The Quimsacocha volcanic center. Earth and Planetary Science Letters, 192: 561–570.

Bijwaard H. 1999. Seismic travel-time tomography for detailed global mantle structure. Utrecht University PhD thesis.

Bittner D and Schmeling H. 1995. Numerical modeling of melting processes and induced diapirism in the lower crust. Geophysical Journal International, 123: 59–70.

Bostock M G and Vandecar J C. 1995. Upper mantle structure of the northern Cascadia subduction zone. Canadian Journal of Earth Sciences, 32: 1–12.

Brace W F and Kohlstedt D L. 1980. Limits on lithospheric stress imposed by laboratory experiments. Journal of Geophysical Research: Solid Earth, 85: 6248–6252.

Brocher T M, Fuis G S, Fisher M A, Plafker G, Moses M J, Taber J J and Christensen N I. 1994. Mapping the Megathrust beneath the Northern Gulf of Alaska using wide-angle seismic data. Journal of Geophysical Research: Solid Earth, 99: 11663–11685.

Burg J P and Gerya T V. 2005. The role of viscous heating in Barrovian metamorphism of collisional orogens: Thermomechanical models and application to the Lepontine Dome in the Central Alps. Journal of Metamorphic Geology, 23: 75–95.

Clauser C and Huenges E. 1995. Thermal conductivity of rocks and minerals. AGU Reference Shelf, 3: 105-126.

Condie K C. 1997. Plate Tectonics (Forth edition). Butterworth-Heinemann.

Dávila F M and Lithgow-Bertelloni C. 2015. Dynamic uplift during slab flattening. Earth and Planetary Science Letters, 425: 34–43.

Defant M J and Drummond M S. 1993. Mount St. Helens: Potential example of the partial melting of the subducted lithosphere in a volcanic arc. Geology, 21: 547–550.

DeMets C, Gordon R G, Argus D F and Stein S. 1994. Effect of recent revisions to the geomagnetic reversal time scale on estimates of current plate motions. Geophysical Research Letters, 21: 2191–2194.

Espurt N, Funiciello F, Martinod J, Guillaume B, Regard V, Faccenna C and Brusset S. 2008. Flat subduction dynamics and deformation of the South American plate: Insights from analog modeling. Tectonics, 27.

Fuis G S, Moore T E, Plafker G, Brocher T M, Fisher M A, Mooney W D, Nokleberg WJ, Page R A, Beaudoin B C, Christensen N I, Levander A R, Lutter W J, Saltus R W and Ruppert N A. 2008. Trans-Alaska Crustal Transect and continental evolution involving subduction underplating and synchronous foreland thrusting. Geology, 36: 267–270.

Fukao Y, Widiyantoro S and Obayashi M. 2001. Stagnant slabs in the upper and lower mantle transition region. Reviews of Geophysics, 39: 291–323.

Gardner T W, Fisher D M, Morell K D and Cupper M L. 2013. Upper-plate deformation in response to flat slab subduction inboard of the aseismic Cocos Ridge, Osa Peninsula, Costa Rica. Lithosphere, 5: 247–264.

Gerya T V and Yuen D A. 2003a. Characteristics-based marker-in-cell method with conservative finite-differences schemes for modeling geological flows with strongly variable transport properties. Physics of the Earth and Planetary Interiors, 140: 293–318.

Gerya T V and Yuen D A. 2003b. Rayleigh-Taylor instabilities from hydration and melting propel ‘cold plumes’ at subduction zones. Earth and Planetary Science Letters, 212: 47–62.

Gorbatov A, Widiyantoro S, Fukao Y and Gordeev E. 2000. Signature of remnant slabs in the North Pacific from P-wave tomography. Geophysical Journal International, 142: 27–36.

Grafe K, Frisch W, Villa I M and Meschede M. 2002. Geodynamic evolution of southern Costa Rica related to low-angle subduction of the Cocos Ridge: Constraints from thermochronology. Tectonophysics, 348: 187–204.

Gripp A E and Gordon R G. 2002. Young tracks of hotspots and current plate velocities. Geophysical Journal International, 150: 321–361.

Gutscher M A and Lallemand S. 1999. Birth of a major strike-slip fault in SW Japan. Terra Nova, 11: 203–209.

Gutscher M A. 2002. Andean subduction styles and their effect on thermal structure and interplate coupling. Journal of South American Earth Sciences, 15: 3–10.

Gutscher M A, Maury R, Eissen J P and Bourdon E. 2000a. Can slab melting be caused by flat subduction? Geology, 28: 535–538.

Gutscher M A, Spakman W, Bijwaard H and Engdahl E R. 2000b. Geodynamics of flat subduction: Seismicity and tomographic constraints from the Andean margin. Tectonics, 19: 814–833.

Hayes G P, Wald D J and Johnson R L. 2012. Slab1.0: A three-dimensional model of global subduction zone geometries. Journal of Geophysical Research: Solid Earth, 117.

Jarrard R D. 1986. Relations among subduction parameters. Reviews of Geophysics, 24: 217–284.

Ji S C and Zhao P L. 1993. Flow laws of multiphase rocks calculated from experimental-data on the constituent phases. Earth and Planetary Science Letters, 117: 181–187.

Kay S M and Abbruzzi J M. 1996. Magmatic evidence for Neogene lithospheric evolution of the central Andean “flat-slab” between 30 degrees S and 32 degrees S. Tectonophysics, 259: 15–28.

Kay S M and Mpodozis C. 2001. Central Andean ore deposits linked to evolving shallow subduction systems and thickening crust. GSA today, 11: 4–9.

Kirby S H. 1983. Rheology of the Lithosphere. Reviews of Geophysics, 21: 1458–1487.

Kirby S H and Kronenberg A K. 1987. Rheology of the lithosphere-selected topics. Reviews of Geophysics, 25: 1219–1244.

Lallemand S, Heuret A and Boutelier D. 2005. On the relationships between slab dip, back-arc stress, upper plate absolute motion, and crustal nature in subduction zones. Geochemistry Geophysics Geosystems, 6(9): Q09006, doi: 10.1029/2005GC000917.

Li Z H, Gerya T V and Burg J P. 2010. Influence of tectonic overpressure on P-T paths of HP-UHP rocks incontinental collision zones: Thermomechanical modelling. Journal of Metamorphic Geology, 28: 227–247.

Li Z H, Xu Z Q and 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: 65–77.

Li Z H. 2014. A review on the numerical geodynamic modeling of continental subduction, collision and exhumation. Science China: Earth Sciences, 57, 47–69.

Manea V C, Pérez-Gussinyé M and Manea M. 2012. Chilean flat slab subduction controlled by overriding plate thickness and trench rollback. Geology, 40: 35–38.

Morris P A. 1995. Slab melting as an explanation of Quaternary volcanism and aseismicity in Southwest Japan. Geology, 23: 395–398.

Perez-Gussinye M, Lowry A R, Morgan J P and Tassara A. 2008. Effective elastic thickness variations along the Andean margin and their relationship to subduction geometry. Geochemistry Geophysics Geosystems, 9(2): Q02003, doi: 1020/2007GC001786.

Petford N and Atherton M. 1996. Na-rich partial melts from newly underplated basaltic crust: The Cordillera Blanca Batholith, Peru. Journal of Petrology, 37: 1491–1521.

Pilger R H. 1981. Plate reconstructions, aseismic ridges, and low-angle subduction beneath the Andes. Geological Society of America Bulletin, 92: 448–456.

Preece S J. 1991. Tephrostratigraphy of the Late Cenozoic Gold Hill Loess, Fairbanks Area, Alaska. University of Toronto PhD thesis.

Ranalli G. 1995. Rheology of the Earth. Springer.

Ranalli G and Murphy D C. 1987. Rheological stratification of the lithosphere. Tectonophysics, 132: 281–295.

Rodriguez-Gonzalez J and Negredo A M. 2012. The role of the overriding plate thermal state on slab dip variability and on the occurrence of flat subduction. Geochemistry Geophysics Geosystems, 13(1): Q01002, doi: 10.1029/ 2011GC003859.

Sak P B, Fisher D M, Gardner T W, Marshall J S and LaFemina P C. 2009. Rough crust subduction, forearc kinematics, and Quaternary uplift rates, Costa Rican segment of the Middle American Trench. Geological Society of America Bulletin, 121: 992–1012.

Schmidt M W and Poli S. 1998. Experimentally based water budgets for dehydrating slabs and consequences for arc magma generation. Earth and Planetary Science Letters, 163: 361–379.

Sigloch K, McQuarrie N and Nolet G. 2008. Two-stage subduction history under North America inferred from multiple-frequency tomography. Nature Geoscience, 1: 458–462.

Skinner S M. 2013. Plate Tectonic Constraints on Flat Subduction and Paleomagnetic Constraints on Rifting. California Institute of Technology PhD thesis.

Skinner S M and Clayton R W. 2011. An evaluation of proposed mechanisms of slab flattening in Central Mexico. Pure and Applied geophysics, 168: 1461–1474.

Stevenson D J and Turner J S. 1977. Angle of subduction. Nature, 270: 334–336.

Suárez G, Monfret T, Wittlinger G and David C. 1990. Geometry of subduction and depth of the seismogenic zone in the Guerrero Gap, Mexico. Nature, 345: 336–338.

Tovish A, Schubert G and Luyendyk B P. 1978. Mantle flow pressure and the angle of subduction: Non-Newtonian corner flows. Journal of Geophysical Research, 83: 5892–5898.

Turcotte D L and Schubert G. 2002. Geodynamics. Cambridge University Press.

Ueda K, Gerya T and Sobolev S V. 2008. Subduction initiation by thermal-chemical plumes: Numerical studies. Physics of the Earth and Planetary Interiors, 171: 296–312.

Uyeda S and Kanamori H. 1979. Back-arc opening and the mode of subduction. Journal of Geophysical Research, 84: 1049–1061.

van Hunen J, van den Berg A P and Vlaar N J. 2001. Latent heat effects of the major mantle phase transitions on low-angle subduction. Earth and Planetary Science Letters, 190: 125–135.

van Hunen J, van den Berg A P and Vlaar N J. 2002a. The impact of the South-American plate motion and the Nazca Ridge subduction on the flat subduction below South Peru. Geophysical Research Letters, 29: 35– 1–35–4.

van Hunen J, van den Berg A P and Vlaar N J. 2002b. On the role of subducting oceanic plateaus in the development of shallow flat subduction. Tectonophysics, 352: 317–333.

van Hunen J, van den Berg A P and Vlaar N J. 2004. Various mechanisms to induce present-day shallow flat subduction and implications for the younger Earth: A numerical parameter study. Physics of the Earth and Planetary Interiors, 146: 179–194.

Numerical Modeling of Flat Subduction: Constraints from the Ocean-continent Convergence Velocity

HUANGFU Pengpeng1, 2, 3, WANG Yuejun3, 4*, FAN Weiming4, LI Zhonghai2, WANG Yuming5and ZHOU Yongzhi3
(1. State Key Laboratory of Isotope Geochemistry, Guangzhou Institute of Geochemistry, Chinese Academy of Sciences, Guangzhou 510640, Guangdong, China; 2. Key Laboratory of Computational Geodynamics, College of Earth Sciences, University of Chinese Academy of Sciences, Beijing 100049, China; 3. School of Earth Science and Geological Engineering, Sun Yat-Sen University, Guangzhou 510275, Guangdong, China; 4. CAS Center for Excellence in Tibetan Plateau Earth Sciences, Institute of Tibetan Plateau Research, Chinese Academy of Sciences, Beijing 100101, China; 5. Zhixin High School, Guangzhou 510080, Guangdong, China)

2-D numerical modeling is used to investigate the effects of ocean-continent convergence velocity on the flat subduction dynamics. It aims to better understand how the absolute subduction velocity of the oceanic lithosphere and absolute overthrusting velocity of the overriding continental lithosphere influence the subduction mode, in particular how the absolute overthrusting velocity regulates the development of flat subduction. Numerical models are divided into two groups with the initial subduction angle of 19° and 36°, respectively. Numerical models with 40 Ma-aged oceanic lithosphere and initial subduction angle of 19° show that (1) under the condition of low oceanic subduction velocity (≤3 cm/a), the flat subduction can be easily formed when the overthrusting velocity is no less than 1 cm/a; (2) under the condition of moderate to high oceanic subduction velocity (>3 cm/a), the flat subduction can occur only when the absolute overriding velocity is no less than the absolute subduction velocity. However, for models with initial subduction angle of 36°, the flat subduction occurs only when the absolute subduction velocity is no larger than 3 cm/a and the absolute overthrusting velocity reaches 10 cm/a. Large absolute subduction velocity can promote decoupling between the downgoing and overriding plates and thus favors the development of steep subduction. In contrast, larger absolute overthrusting velocity facilitates the slab bending and thereby contributes to the development of flat subduction. Additionally, flat subduction at a high velocity of the overriding plate is significantly correlated with the landward mantle flow beneath the slab, which can facilitate the upward bending of the slab. Our results are evidenced by the kinematic models of the flat subduction beneath the South America.

flat subduction; numerical simulation; absolute subduction velocity; absolute overthrusting velocity; initial subduction angle

P542; TE352

A

1001-1552(2016)03-0429-017

2015-11-16; 改回日期: 2016-02-05
項目資助: 973項目(2014CB440901和 2015CB856106)、國家自然科學(xué)基金項目(41190073, 41402165和 41304071)和中山大學(xué)中央高?；究蒲袠I(yè)務(wù)費專項資金聯(lián)合資助。

皇甫鵬鵬(1985–), 男, 博士研究生, 從事數(shù)值模擬研究。Email: huangfu@gig.ac.cn

王岳軍(1969–), 男, 教授, 從事大地構(gòu)造研究。Email: wangyuejun@mail.sysu.edu.cn