王瑞澤, 王建超, 金宗瑋, 劉駿標(biāo), 王雪純, 靳國棟, 邢會林, 2*

俯沖板塊脫水相變對板塊俯沖過程的影響

王瑞澤1, 王建超1, 金宗瑋1, 劉駿標(biāo)1, 王雪純1, 靳國棟1, 邢會林1, 2*

(1. 深海圈層與地球系統(tǒng)教育部前沿科學(xué)中心, 海底科學(xué)與探測技術(shù)教育部重點實驗室, 中國海洋大學(xué) 海洋地球科學(xué)學(xué)院, 山東 青島 266100; 2. 嶗山實驗室, 山東 青島 266237)

板塊俯沖是地球上最重要的動力學(xué)過程之一, 關(guān)于俯沖帶相變對地幔楔和上覆板塊影響前人已開展大量研究, 但關(guān)于相變反應(yīng)對俯沖板塊的動力學(xué)影響研究較少。本研究在拉格朗日積分點有限元代碼的基礎(chǔ)上, 實現(xiàn)了俯沖板塊的相變, 建立了二維自由俯沖模型, 研究了俯沖過程中俯沖板塊上礦物脫水相變對俯沖速度、俯沖板塊形態(tài)、板塊上的應(yīng)力變化的影響。同時我們研究了薄弱帶的黏度變化對俯沖起始與俯沖過程的影響。數(shù)值模擬實驗結(jié)果顯示: ①俯沖板塊脫水相變現(xiàn)象的發(fā)生降低了板塊俯沖的速度, 改變了俯沖板塊的形態(tài), 在俯沖板塊的相變界面處發(fā)生應(yīng)力累積; ②相較于礦物相變導(dǎo)致的黏度增大, 俯沖起始角度仍是影響俯沖板塊俯沖速度的主要因素; ③薄弱帶黏度越高越不利于俯沖板塊與上覆板塊的解耦, 抑制了俯沖過程的進行。該研究為探究板塊回卷的控制因素以及俯沖板塊應(yīng)力狀態(tài)變化的影響因素提供了新的視角。

俯沖板塊; 俯沖過程; 礦物相變; 有限元數(shù)值模擬

0 引 言

俯沖帶是地幔對流單元的下降分支, 是地球內(nèi)部占主導(dǎo)地位的物理和化學(xué)系統(tǒng), 也是地球上最大的“回收系統(tǒng)”(Stern, 2002)。俯沖板塊將遠(yuǎn)洋陸源沉積物、蝕變的和新鮮的玄武質(zhì)洋殼物質(zhì)以及地幔巖石圈等“原材料”運送到俯沖帶, 使海洋巖石圈、沉積物和海水與周圍的地幔物質(zhì)發(fā)生熔融, 重新平衡(Tatsumi, 2005)。

隨著俯沖深度的增加, 俯沖帶中溫度和壓力不斷增大。俯沖帶熱力學(xué)模擬是研究俯沖帶熱結(jié)構(gòu)最好的方法, 通過計算發(fā)現(xiàn)板塊在俯沖過程中被高溫的地幔加熱(Peacock and Wang, 1999; Gerya et al., 2002; Gorman et al., 2006; Syracuse et al., 2010)。同時, 研究證明, 在高溫高壓條件下, 俯沖板塊中礦物主要發(fā)生兩類相變反應(yīng): 一類為俯沖板塊上深度較淺的玄武巖?榴輝巖系列相變(Schmidt and Poli, 1998; Fryer et al., 1999; Forneris and Holloway, 2003; Forneris and Holloway, 2004); 另一類則是發(fā)生在俯沖地幔以及俯沖板塊上覆地幔楔中的橄欖巖?蛇紋巖?橄欖巖系列相變。蛇紋石化是一個重要的地球動力學(xué)過程, 它極大地改變了地幔和洋殼的物理和力學(xué)性質(zhì)。蛇紋巖是由超基性巖和基性巖在俯沖帶淺部、低溫(100 ℃)?中溫(700 ℃)條件下水合而成(Bose and Ganguly, 1995;Ulmer and Trommsdorff, 1995; Kawamoto and Holloway, 1997; Grove et al., 2006; Till et al., 2012; Green et al., 2014)。隨著俯沖深度的增加, 在高溫高壓條件下, 蛇紋巖失穩(wěn)發(fā)生脫水相變反應(yīng), 形成橄欖石與其他礦物組合(Ulmer and Trommsdorff, 1995; Wunder and Schreyer, 1997; Evans, 2004; Perrillat et al., 2005; Schwartz et al., 2013; Guillot et al., 2015)。

由于俯沖帶巖石樣品獲取難度較大, 而實驗又難以解決俯沖過程中復(fù)雜的地質(zhì)問題。隨著計算機技術(shù)的進步, 數(shù)值模擬可以有效地解決地質(zhì)問題中存在的時空尺度局限性問題, 因此許多學(xué)者通過數(shù)值模擬手段對俯沖帶的相變反應(yīng)及其影響進行了研究。Gerya et al. (2002)結(jié)合前人巖石學(xué)成果, 基于有限差分方法, 針對俯沖通道的發(fā)展和地幔楔水化作用這一問題進行了模擬研究; Arcay et al. (2005)通過有限元數(shù)值模擬的方法, 研究俯沖板片脫水對地幔楔的影響; Wada et al. (2008)通過有限元數(shù)值模擬, 研究上覆板塊與俯沖板塊之間薄弱帶以及板塊脫水對地幔楔的影響; Gerya and Meilick (2011)通過有限差分程序, 研究俯沖帶熔體和流體對流變性質(zhì)的影響, 以及其對俯沖動力學(xué)的影響; Faccenda et al. (2012)基于有限差分方法, 研究板塊相變脫水對中?深地震活動、板塊弱化和深部水循環(huán)的影響; Liao et al. (2017)利用同樣的方法, 研究水化作用引起的上地幔弱化對克拉通巖石圈動力學(xué)的影響。

盡管前人對俯沖板塊的相變過程以及相變脫水對地幔物質(zhì)的影響進行了大量研究, 但主要側(cè)重于相變反應(yīng)對地幔楔以及上覆板塊的影響, 而很少考慮俯沖過程中礦物脫水相變對俯沖板塊動力學(xué)行為的影響(Schellart et al., 2010; Sandiford et al., 2019)。為了充分研究板塊俯沖至平衡深度時發(fā)生蛇紋巖巖石圈地幔相變對俯沖速度、俯沖板片形態(tài)的影響, 以及對俯沖板塊應(yīng)力應(yīng)變狀態(tài)的影響, 本次研究基于有限元數(shù)值模擬方法, 通過構(gòu)建俯沖板塊自由俯沖模型, 對俯沖板塊中各種材料參數(shù)進行設(shè)置, 對板塊俯沖過程進行模擬, 并結(jié)合前人研究成果, 進一步討論礦物相變對板塊俯沖的影響。

1 控制方程

1.1 有限元部分

通過求解穩(wěn)態(tài)不可壓縮蠕動流的斯托克斯方程進行數(shù)值模擬。

其中質(zhì)量守恒方程為:

動量守恒方程為:

對于不可壓縮流, 可通過公式(3)計算偏應(yīng)力張量:

本次研究采用拉格朗日積分點有限元法進行有限元數(shù)值模擬計算。拉格朗日積分點有限元法主要優(yōu)點在于: ①能在不顯著改變精度的情況下繼續(xù)進行極大變形; ②能夠跟蹤材料歷史和界面, 其精度高; ③在整個運行過程中保留規(guī)則網(wǎng)格, 可以使用快速數(shù)值求解器。求解域由歐拉網(wǎng)格和嵌入的拉格朗日追蹤點或粒子表示, 未知變量在網(wǎng)格節(jié)點上計算, 拉格朗日粒子在變形過程中攜帶歷史變量。該方法非常適合于模擬地質(zhì)環(huán)境中經(jīng)常遇到的連續(xù)介質(zhì)固體的類流體行為, 能夠?qū)Πǖ蒯α鳌鍓K俯沖和大陸碰撞在內(nèi)的構(gòu)造過程進行二維/三維熱?變形?化學(xué)反應(yīng)建模(Moresi et al., 2003, 2007)。

1.2 相變部分

葉蛇紋石是俯沖帶高溫高壓蛇紋石的變種, 是深度大于100 km俯沖帶中最主要的含水礦物, 俯沖帶中主要脫水相變過程(圖1)(Guillot et al., 2015)為

圖1 根據(jù)實驗數(shù)據(jù)和天然蛇紋石巖石學(xué)觀測得出的蛇紋石穩(wěn)定場(Perrillat et al., 2005; Guillot et al., 2015)

Fig.1 Serpentine stability fields based on experimental data and natural serpentine lithological observations

據(jù)Carminati et al. (1999)和Schwartz et al. (2001)的研究, 俯沖帶蛇紋巖脫水相變會導(dǎo)致巖石流變強度變大。本次研究基于有限元理論、固定的計算網(wǎng)格和拉格朗日積分點法的程序Underworld, 通過判斷板塊上的追蹤粒子是否發(fā)生相變進而改變材料黏度, 在程序中實現(xiàn)了俯沖板片上蛇紋巖脫水相變這一過程(圖2)。

圖2 相變部分算法計算流程圖

Fig.2 The algorithm flow chart of phase transition

圖3 不同俯沖起始角度模型幾何示意圖

2 模型設(shè)置

2.1 模型的幾何和流變學(xué)

在海洋板塊彎曲前端與大陸板塊之間設(shè)置寬度約10 km的薄弱帶, 初始模型中黏度設(shè)置為上地幔黏度的十分之一, 上覆板塊厚度設(shè)置為90 km, 與俯沖板塊俯沖起始深度相同。本次模擬實驗只考慮主要的相變過程(即蛇紋巖中葉蛇紋石失穩(wěn)脫水相變, 形成橄欖石為主以及少量蛇紋石的巖石組分)對俯沖板塊的流變學(xué)性質(zhì)的影響, 模型中各種材料參數(shù)如表1所示。

表1 模型使用材料參數(shù)表

2.2 邊界條件

本研究假設(shè)成熟洋板塊和上地幔之間的重力不穩(wěn)定性導(dǎo)致了俯沖板塊的下沉, 俯沖是由板塊前端的彎曲尖端觸發(fā)。在模型的邊界施加自由滑動速度邊界。因此, 只有負(fù)浮力(板塊拉張)才能驅(qū)使海洋板塊俯沖, 誘發(fā)地幔流動, 沒有速度或應(yīng)力施加到俯沖板塊或上覆板塊上。

2.3 模型設(shè)計

本模擬設(shè)置了俯沖起始角度為30°的對照組模型1～2, 不同俯沖起始角度的實驗組模型3～6(表2)。對照組為在重力作用下板塊發(fā)生俯沖, 俯沖過程中俯沖板塊、上地幔等不同材料的黏度保持不變; 實驗組為在重力作用下板塊發(fā)生俯沖, 在俯沖至固定深度時, 俯沖板塊會發(fā)生相應(yīng)的相變, 從而導(dǎo)致黏度變化。本次研究主要探究板塊俯沖過程中相變對俯沖板塊內(nèi)應(yīng)力場的變化、俯沖板塊的動力學(xué)行為的影響。同時根據(jù)模擬結(jié)果, 增設(shè)兩組模型以探討不同的薄弱帶黏度對俯沖板塊俯沖的影響(模型7、8), 以及兩組模型以探討不同薄弱帶寬度對俯沖板塊俯沖過程的影響(模型9、10)。

表2 模型參數(shù)設(shè)置表

3 結(jié)果分析

3.1 俯沖過程分析(模型1、2)

板塊俯沖運動是由板塊自身重力和克服板塊底部黏滯阻力的能力決定的, 其中下降板的負(fù)浮力是板塊運動的主要動力。板塊的自由俯沖模型結(jié)果顯示俯沖具有兩個階段。

階段一: 由于地幔和大洋板塊之間的重力不穩(wěn)定, 大洋板塊從預(yù)先設(shè)置的俯沖彎曲前端擾動處開始俯沖(圖4a、e)。隨著板塊俯沖, 上覆板塊持續(xù)向海溝方向后退。上地幔在俯沖作用下發(fā)生擾動, 以平行于俯沖板塊下方海溝的方向運動(圖4b、d、f、g)。俯沖板塊后撤引起的水平向地幔角流既驅(qū)動上覆板塊向海溝側(cè)運動, 同時也拖拽俯沖板塊持續(xù)向下俯沖(Zhang et al., 2017)。階段二: 在俯沖進行一段時間后, 自由俯沖模型中的俯沖板塊在重力作用下以接近90°向地幔深處俯沖(圖4d、h)。

(a)～(d) 為板塊未發(fā)生相變俯沖至50 Ma結(jié)果圖; (e)～(h) 為板塊發(fā)生相變俯沖至50 Ma結(jié)果圖。白色箭頭標(biāo)記為均方根速度(Vrms)。

當(dāng)俯沖起始角度為30°時, 巖石圈地幔未發(fā)生相變時, 整個俯沖板塊黏度均勻, 因此整個板塊俯沖趨勢保持一致。俯沖彎曲前端在自身重力的作用下, 俯沖角度由起始時的30°逐漸變?yōu)?0°垂直向下俯沖。若在模型中添加俯沖過程中蛇紋石化巖石圈地幔的相變反應(yīng)后, 俯沖過程則變得有所不同。由于大洋巖石圈在俯沖過程中發(fā)生了相變, 導(dǎo)致黏度變大, 進而使得俯沖更難以進行, 因此俯沖角度從30°至90°所需時間明顯少于未發(fā)生相變的模型。

3.2 蛇紋石化巖石圈地幔相變對俯沖速度的影響

為了研究不同俯沖起始角度對俯沖速度的影響, 以及蛇紋石化巖石圈地幔發(fā)生相變對俯沖速度帶來的影響, 本次研究設(shè)置了模型3～6。首先, 設(shè)定俯沖起始角度為30°, 俯沖至同一時刻發(fā)生相變與未發(fā)生相變的模型(模型1、2), 并進行了均方根速度成圖(圖5)。結(jié)果顯示, 板塊俯沖至相同時刻下, 未發(fā)生相變的模型中俯沖板塊處的均方根速度明顯大于發(fā)生相變模型。由于板塊俯沖是在板塊的負(fù)浮力以及板塊與地幔之間的黏滯阻力共同作用下完成的, 而發(fā)生相變俯沖板塊底部需要克服的黏滯阻力變大, 導(dǎo)致俯沖板塊俯沖速度變小。

圖5 俯沖起始角度為30°模型俯沖至20.3 Ma時均方根速度圖

模擬中通過追蹤不同俯沖起始角度以及是否發(fā)生相變影響下, 各組模型俯沖板塊最底部粒子所在深度, 得到模型1～6俯沖深度與時間的關(guān)系。通過對未發(fā)生相變、俯沖起始角度不同的模型1、3、5進行分析(圖6), 結(jié)果顯示,在相同俯沖時間內(nèi), 俯沖起始角度為30°的模型俯沖速度最快, 隨著俯沖起始角度的增加, 俯沖速度相應(yīng)變慢。在相同的俯沖時間內(nèi), 俯沖起始角度越小, 俯沖板塊俯沖所能達到的深度越深。

圖6 不同俯沖起始角度未發(fā)生相變模型俯沖速度圖

發(fā)生相變、俯沖起始角度不同的模型2、4、6俯沖深度隨俯沖時間變化的結(jié)果(圖7)顯示, 板塊相變在一定程度上影響了板塊俯沖的速度, 導(dǎo)致板塊俯沖速度變慢。但是隨著俯沖起始角度的增大, 板塊相變對俯沖速度的影響逐漸減小, 因此, 俯沖起始角度是控制俯沖板塊俯沖速度的主要因素。

圖7 不同俯沖起始角度發(fā)生相變模型俯沖速度圖

3.3 俯沖過程應(yīng)力變化分析

圖8 不同俯沖起始角度發(fā)生相變與未發(fā)生相變模型板塊上最大主應(yīng)力的相對變化

3.4 薄弱帶作用結(jié)果

為了研究薄弱帶的黏度和寬度在俯沖過程中的作用, 本次研究設(shè)置了模型7～10。當(dāng)模型中薄弱帶固定寬度為10 km, 薄弱帶黏度為低于俯沖板塊與上覆板塊的黏度1×1016Pa·s時, 俯沖板塊彎曲前端可以較快速的與上覆板塊發(fā)生解耦, 在自身重力的作用下, 耗時約23 Ma進入俯沖角度近乎90°的垂直俯沖階段, 并在俯沖后期發(fā)生板塊回卷現(xiàn)象(圖4h)。而增設(shè)兩組更大黏度的模型7和模型8則表明, 隨著薄弱帶流變強度增大, 俯沖角度從起始的30°變化至90°所需的時間變長。其中薄弱帶黏度為1×1019Pa·s的模型7(圖9a、b); 即薄弱帶不存在, 俯沖洋殼與流變強度相同的上覆地殼直接接觸時, 耗時約28 Ma變?yōu)楦_角度接近90°的垂直向下俯沖; 而薄弱帶黏度為1×1022Pa·s的模型8(圖9c、d), 則在俯沖角度俯沖至約45°時, 俯沖趨勢就被上部高黏度的薄弱帶所抑制, 黏度比俯沖板塊和上覆板塊高的薄弱帶不利于俯沖板塊與上覆板塊之間的解耦, 抑制了俯沖過程的進行。

圖9 不同薄弱帶黏度設(shè)置下板塊俯沖模型黏度場與速度場

為探究不同的薄弱帶寬度對俯沖過程的影響, 本次研究增設(shè)了薄弱帶寬度為5 km的模型9和薄弱帶寬度為20 km的模型10。結(jié)果顯示, 在薄弱帶黏度固定為1×1016Pa·s時, 薄弱帶寬度增大時, 板塊俯沖速度明顯增大(圖10)。通過對比不同薄弱帶寬度模型中俯沖起始角度由30°變?yōu)榇怪备_階段所用的時間, 顯示薄弱帶寬度為20 km時, 所需的時間約為15 Ma; 薄弱帶寬度為10 km時, 所需時間約為23 Ma; 薄弱帶寬度為5 km時, 所需時間約為27 Ma。因此, 薄弱帶黏度固定且為較低黏度時, 薄弱帶寬度越大, 俯沖板塊與上覆板塊之間解耦越快。

圖10 不同薄弱帶寬度下發(fā)生相變模型俯沖速度圖

4 討 論

本次研究在拉格朗日積分點有限元方法的基礎(chǔ)上, 添加了與俯沖板塊蛇紋巖脫水相變相關(guān)的代碼, 研究了礦物相變對俯沖過程產(chǎn)生的一系列影響。但是由于本次研究的模型中只考慮材料黏度變化對俯沖板塊的影響, 因此只考慮了俯沖板塊中蛇紋巖脫水相變導(dǎo)致的黏度差異對俯沖過程的影響, 而認(rèn)為俯沖板塊中大洋中脊玄武巖系列的相變反應(yīng)對俯沖板塊的影響較小。

4.1 板塊俯沖形態(tài)

在前人研究的基礎(chǔ)上(Zhang et al., 2017), 選擇的基礎(chǔ)模型為自由俯沖模型, 沒有額外恒定的邊界條件, 所以板塊不能保持傾角恒定的持續(xù)俯沖, 表現(xiàn)出俯沖角度從30°逐漸向90°增加。而在模型2中, 由于俯沖板塊120 km以下板塊發(fā)生相變, 俯沖角度變化更大, 出現(xiàn)類似于Perchuk et al. (2020)模擬結(jié)果中的板塊回卷現(xiàn)象(圖4h)。前人研究認(rèn)為榴輝巖化導(dǎo)致的重力異常是引起板塊回卷的主要原因之一(孫圣思和嵇少丞, 2011; 皇甫鵬鵬等, 2016; Dai et al., 2020)。而在本次模擬中由于設(shè)定密度固定不變, 只有板塊俯沖至一定深度時黏度發(fā)生變化, 俯沖板塊才會發(fā)生了類似板塊回卷現(xiàn)象, 所以俯沖板塊相變導(dǎo)致板塊黏度變大可能也是導(dǎo)致板塊回卷現(xiàn)象的原因之一。

4.2 板塊應(yīng)力變化

模擬結(jié)果顯示, 板塊俯沖至不同深度時, 在板塊發(fā)生相變的界面處出現(xiàn)應(yīng)力集中導(dǎo)致應(yīng)力陡變。應(yīng)力集中出現(xiàn)的原因是在俯沖過程中, 俯沖板塊相變導(dǎo)致黏滯力增大, 因此板塊自相變界面以下部分俯沖速度降低, 發(fā)生了短暫的停滯。Kita et al. (2006)發(fā)現(xiàn)日本東北部下方70～100 km深度處雙地震帶上平面存在帶狀地震活動集中區(qū), 并將其解釋為與Hacker et al. (2003)模擬的相邊界附近的脫水反應(yīng)有關(guān); Tsuji et al. (2008)利用雙差層析成像技術(shù), 估算出日本東部太平洋板塊下一個厚度約為10 km的含水低速帶在70～90 km處逐漸發(fā)生脫水消失。本次模擬結(jié)果也證實了相變界面處發(fā)生了應(yīng)力集中。因此, 俯沖板塊黏度變化導(dǎo)致的應(yīng)力集中可能是導(dǎo)致雙地震帶上形成帶狀地震活動集中區(qū)的原因之一。

4.3 薄弱帶的影響

俯沖帶的數(shù)值模擬實驗中, 俯沖洋殼和上覆板塊之間的薄弱帶黏度的大小會影響板塊速度(Androvi?ová et al., 2013; Holt et al., 2015), 較高黏度的薄弱帶會抑制板塊回卷現(xiàn)象的發(fā)生(Arredondo and Billen, 2017)。因此, 俯沖洋殼和上覆板塊的解耦對俯沖板塊深部變形具有很強的敏感性。如圖4h、9b、9d所示, 較低的薄弱帶黏度導(dǎo)致板塊回卷較容易發(fā)生, 而較強的薄弱帶黏度導(dǎo)致板塊回卷現(xiàn)象難以發(fā)生。同時, 薄弱帶的寬度也會對俯沖板塊與上覆板塊的解耦產(chǎn)生影響, 在薄弱帶黏度較低時, 越寬的薄弱帶越有利于俯沖板塊與上覆板塊的解耦。

5 結(jié) 論

本次研究通過在拉格朗日積分點有限元代碼中添加俯沖板塊礦物相變部分, 建立俯沖板塊自由俯沖相變模型, 模擬了相變對俯沖板塊的動力學(xué)行為與應(yīng)力狀態(tài)的影響, 以及不同的俯沖角度與薄弱帶對俯沖過程的影響, 得到以下幾點認(rèn)識:

(1) 礦物相變導(dǎo)致了俯沖板塊上發(fā)生相變的部分黏度變大, 進而使板塊俯沖需要克服的板塊底部黏滯阻力變大, 從而降低了板塊俯沖的速度, 改變了俯沖板塊的俯沖形態(tài), 同時在俯沖板塊的相變界面處發(fā)生應(yīng)力累積。

(2) 雖然俯沖板塊上相變降低了其俯沖速度, 但俯沖起始角度才是影響俯沖板塊俯沖速度的主要因素。

(3) 當(dāng)薄弱帶寬度一定時, 黏度較低的薄弱帶有利于上覆板塊與俯沖板塊之間的解耦, 薄弱帶黏度越高越不利于俯沖板塊與上覆板塊的解耦, 在一定程度上抑制俯沖過程的進行; 當(dāng)薄弱帶黏度(較低黏度)一定時, 薄弱帶的寬度越大越有利于俯沖板塊與上覆板塊的解耦, 反之則會對俯沖過程起到抑制作用。

致謝:中國科學(xué)院廣州地球化學(xué)研究所郭鋒研究員及另一名匿名審稿人對本文提出了建設(shè)性修改意見, 在此致以特別感謝。

皇甫鵬鵬, 王岳軍, 范蔚茗, 李忠海, 王喻鳴, 周永智. 2016. 大洋平板俯沖的數(shù)值模擬再現(xiàn): 洋?陸匯聚速率影響. 大地構(gòu)造與成礦學(xué), 40(3): 429–445.

李三忠, 王光增, 索艷慧, 李璽瑤, 戴黎明, 劉一鳴, 周潔, 郭玲莉, 劉永江, 張國偉. 2019. 板塊驅(qū)動力: 問題本源與本質(zhì). 大地構(gòu)造與成礦學(xué), 43(4): 605–643.

孫圣思, 嵇少丞. 2011. 大洋板塊俯沖帶地震波各向異性及剪切波分裂的成因機制. 大地構(gòu)造與成礦學(xué), 35(4): 628–647.

Androvi?ová A, ?í?ková H, van den Berg A. 2013. The effects of rheological decoupling on slab deformation in the Earth’s upper mantle., 57(3): 460–481.

Arcay D, Tric E, Doin M. 2005. Numerical simulations of subduction zones: Effect of slab dehydration on the mantle wedge dynamics., 149(1–2): 133–153.

Arredondo K M, Billen M I. 2017. Coupled effects of phase transitions and rheology in 2-D dynamical models of subduction.:, 122(7): 5813–5830.

Bose K, Ganguly J. 1995. Experimental and theoretical studies of the stabilities of talc, antigorite and phase A at high pressures with applications to subduction processes., 136(3–4): 109–121.

Carminati E, Giunchi C, Argnani A, Sabadini R, Fernandez M. 1999. Plio-Quaternary vertical motion of the Northern Apennines: Insights from dynamic modeling., 18(4): 703–718.

Cloos M. 1993. Lithospheric buoyancy and collisional orogenesis: Subduction of oceanic plateaus, continental margins, island arcs, spreading ridges, and seamounts., 105(6): 715–737.

Dai L M, Wang L L, Lou D, Li Z H, Dong H, Ma F F, Li F K, Li S Z, Yu S Y. 2020. Slab rollback versus delamination: Contrasting fates of flat-slab subduction and implicationsfor South China evolution in the Mesozoic.:, 125(4), e2019JB019164.

Duretz T, Gerya T V, May D A. 2011. Numerical modelling of spontaneous slab breakoff and subsequent topographic response., 502(1): 244–256.

Evans B W. 2004. The serpentinite multisystem revisited: Chrysotile is metastable., 46(6): 479–506.

Faccenda M, Gerya T V, Mancktelow N S, Moresi L. 2012. Fluid flow during slab unbending and dehydration: Implications for intermediate-depth seismicity, slab weakening and deep water recycling.,,, 13(1), Q01010.

Forneris J F, Holloway J R. 2003. Phase equilibria in subductingbasaltic crust: Implications for H2O release from the slab., 214(1–2): 187–201.

Forneris J F, Holloway J R. 2004. Evolution of mineral compositions during eclogitization of subducting basaltic crust., 89(10): 1516–1524.

Fryer P, Wheat C G, Mottl M J. 1999. Mariana blueschist mud volcanism: Implications for conditions within the subduction zone., 27(2): 103–106.

Gerya T V, Meilick F I. 2011. Geodynamic regimes of subduction under an active margin: Effects of rheologicalweakening by fluids and melts., 29(1): 7–31.

Gerya T V, St?ckhert B, Perchuk A L. 2002. Exhumation of high-pressure metamorphic rocks in a subduction channel: A numerical simulation., 21(6), 1056.

Gerya T V, Yuen D A. 2003. Rayleigh-Taylor instabilities from hydration and melting propel ‘cold plumes’ at subduction zones., 212(1): 47–62.

Gorman P J, Kerrick D M, Connolly J. 2006. Modeling open system metamorphic decarbonation of subducting slabs.,,, 7(4), Q04007.

Green D H, Hibberson W O, Rosenthal A, Kovács I, Yaxley G M, Falloon T J, Brink F. 2014. Experimental study of the influence of water on melting and phase assemblagesin the upper mantle., 55(10): 2067– 2096.

Grove T L, Chatterjee N, Parman S W, Médard E. 2006. The influence of H2O on mantle wedge melting., 249(1–2): 74–89.

Guillot S, Schwartz S, Reynard B, Agard P, Prigent C. 2015. Tectonic significance of serpentinites., 646: 1–19.

Hacker B R, Peacock S M, Abers G A, Holloway S D. 2003. Subduction factory 2. Are intermediate-depth earthquakes in subducting slabs linked to metamorphic dehydration reactions?:, 108(B1), ESE11.

Holt A F, Becker T W, Buffett B A. 2015. Trench migration and overriding plate stress in dynamic subduction models., 201(1): 172–192.

Kawamoto T, Holloway J R. 1997. Melting temperature and partial melt chemistry of H2O-saturated mantle peridotite to 11 gigapascals., 276(5310): 240–243.

Kita S, Okada T, Nakajima J, Matsuzawa T, Hasegawa A. 2006. Existence of a seismic belt in the upper plane of the double seismic zone extending in the along-arc direction at depths of 70–100 km beneath NE Japan., 33(24), L24310.1- L24310.5.

Korenaga T, Korenaga J. 2016. Evolution of young oceanic lithosphere and the meaning of seafloor subsidence rate.:, 121(9): 6315–6332.

Liao J, Wang Q, Gerya T, Ballmer M D. 2017. Modeling craton destruction by hydration-induced weakening of the upper mantle.:, 122(9): 7449–7466.

Moresi L, Dufour F, Mühlhaus H. 2003. A Lagrangian integrationpoint finite element method for large deformation modeling of viscoelastic geomaterials., 184(2): 476–497.

Moresi L, Quenette S, Lemiale V, Meriaux C, Appelbe B, Mühlhaus H B. 2007. Computational approaches to studying non-linear dynamics of the crust and mantle., 163(1–4): 69–82.

Peacock S M, Wang K, 1999. Seismic consequences of warm versus cool subduction metamorphism: Examples from southwest and northeast Japan., 286(5441): 937–939.

Perchuk A L, Gerya T V, Zakharov V S, Griffin W L. 2020. Building cratonic keels in Precambrian plate tectonics., 586(7829): 395–401.

Perrillat J P, Daniel I, Koga K T, Reynard B, Cardon H, Crichton W A. 2005. Kinetics of antigorite dehydration: A real-time X-ray diffraction study., 236(3–4): 899–913.

Sandiford D, Moresi L, Sandiford M, Yang T. 2019. Geometric controls on flat slab seismicity., 527: 115787.

Schellart W P, Stegman D R, Farrington R J, Freeman J, Moresi L. 2010. Cenozoic tectonics of western North America controlled by evolving width of Farallon slab., 329(5989): 316–319.

Schmidt M W, Poli S. 1998. Experimentally based water budgets for dehydrating slabs and consequences for arc magma generation., 163(1–4): 361–379.

Schwartz S, Allemand P, Guillot S. 2001. Numerical model of the effect of serpentinites on the exhumation of eclogitic rocks: Insights from the Monviso ophiolitic massif (Western Alps)., 342(1–2): 193–206.

Schwartz S, Guillot S, Reynard B, Lafay R, Debret B, Nicollet C, Lanari P, Auzende A L. 2013. Pressure- temperature estimates of the lizardite/antigorite transition in high pressure serpentinites., 178: 197–210.

Stern R J. 2002. Subduction zones., 40(4): 3–1.

Syracuse E M, van Keken P E, Abers G A. 2010. The global range of subduction zone thermal models., 183(1–2): 73–90.

Tatsumi Y. 2005. The subduction factory: How it operates in the evolving Earth., 15(7): 4.

Till C B, Grove T L, Withers A C. 2012. The beginnings of hydrous mantle wedge melting., 163(4): 669–688.

Tsuji Y, Nakajima J, Hasegawa A. 2008. Tomographic evidence for hydrated oceanic crust of the Pacific slab beneath northeastern Japan: Implications for water transportation in subduction zones., 35(14), L14308.

Ulmer P, Trommsdorff V. 1995. Serpentine stability to mantle depths and subduction-related magmatism., 268(5212): 858–861.

Wada I, Wang K, He J, Hyndman R D. 2008. Weakening of the subduction interface and its effects on surface heat flow, slab dehydration, and mantle wedge serpentinization.:, 113(B4), B04402.

Wunder B, Schreyer W. 1997. Antigorite: High-pressure stability in the system MgO-SiO2-H2O (MSH)., 41(1–3): 213–227.

Zhang Q W, Guo F, Zhao L, Wu Y M. 2017. Geodynamics of divergent double subduction: 3-D numerical modeling of a Cenozoic example in the Molucca Sea region, Indonesia.:, 122(5): 3977–3998.

Influence of Dehydration Phase Transition of Subducting Plate on the Process of Plate Subduction

WANG Ruize1, WANG Jianchao1, JIN Zongwei1, LIU Junbiao1, WANG Xuechun1, JIN Guodong1, XING Huilin1, 2*

(1. Frontiers Science Center for Deep Ocean Multispheres and Earth System, MOE Key Lab of Submarine Geosciences and Prospecting Techniques, College of Marine Geosciences, Ocean University of China, Qingdao 266100, Shandong, China; 2. Laoshan Laboratory, Qingdao 266237, Shandong, China)

Plate subduction is one of the most important dynamic processes on the earth. The influences of phase change in subduction zones on the mantle wedge and overriding plates have been studied extensively, but researches on the influences of phase change on the dynamics of subducting plates are rare. Based on the Lagrangian integration point finite element code, we implemented the phase change of the plate into the code and performed 2D free subduction modelling to simulate the process of subduction. We investigated the influence of dehydration-related phase change on the subduction speed, shape of the subducting plate, and change of the stress field of the plate during the subducting process. Meanwhile we studied the influence of the viscosity change of the weak zone on the initiation and process of subduction. Based on these experiments we concluded that: (1) The phase change caused by mineral dehydration of the subducting plate reduces the speed of plate subducting and changes the shape of the subducting plate, and leads to stress accumulation at the phase change interface of the subducting plate. (2) Compared with the increase in viscosity caused by mineral phase transitions, initial subducting angle is still the main factor affecting the subduction velocity of subducting plates. (3) The higher the viscosity of the weak zone, the more difficult the decoupling of the subducting plate and the overlying plate, and thus inhibiting the progress of the subduction process. The results of simulation provide a new perspective for us to explore the controlling factors of rollback and the influencing factors of changes in the stress state of the subducting plate.

subducting plate; subducting process; phase change; finite element simulation

2021-12-05;

2022-03-28

國家自然科學(xué)基金重大計劃重點支持項目(92058211)、國家自然科學(xué)基金面上項目(52074251)、中央高?；究蒲袠I(yè)務(wù)經(jīng)費(842012003)和高等學(xué)校學(xué)科創(chuàng)新引智計劃項目(B20048)聯(lián)合資助。

王瑞澤(1997–), 男, 碩士研究生, 海洋地質(zhì)專業(yè)。E-mail: ruize1997@163.com

邢會林(1965–), 男, 教授, 博士生導(dǎo)師, 主要從事超級計算地球科學(xué)理論、軟件研發(fā)及其應(yīng)用等研究。E-mail: h.xing@ouc.edu.cn

P541

A

1001-1552(2023)06-1232-010

10.16539/j.ddgzyckx.2023.06.003