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

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

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

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

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

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

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

0 引 言

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

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

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

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

1 控制方程

1.1 有限元部分

通過求解穩態不可壓縮蠕動流的斯托克斯方程進行數值模擬。

其中質量守恒方程為:

動量守恒方程為:

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

本次研究采用拉格朗日積分點有限元法進行有限元數值模擬計算。拉格朗日積分點有限元法主要優點在于: ①能在不顯著改變精度的情況下繼續進行極大變形; ②能夠跟蹤材料歷史和界面, 其精度高; ③在整個運行過程中保留規則網格, 可以使用快速數值求解器。求解域由歐拉網格和嵌入的拉格朗日追蹤點或粒子表示, 未知變量在網格節點上計算, 拉格朗日粒子在變形過程中攜帶歷史變量。該方法非常適合于模擬地質環境中經常遇到的連續介質固體的類流體行為, 能夠對包括地幔對流、板塊俯沖和大陸碰撞在內的構造過程進行二維/三維熱?變形?化學反應建模(Moresi et al., 2003, 2007)。

1.2 相變部分

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

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

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

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

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

Fig.2 The algorithm flow chart of phase transition

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

2 模型設置

2.1 模型的幾何和流變學

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

表1 模型使用材料參數表

2.2 邊界條件

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

2.3 模型設計

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

表2 模型參數設置表

3 結果分析

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

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

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

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

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

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

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

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

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

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

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

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

3.3 俯沖過程應力變化分析

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

3.4 薄弱帶作用結果

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

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

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

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

4 討 論

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

4.1 板塊俯沖形態

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

4.2 板塊應力變化

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

4.3 薄弱帶的影響

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

5 結 論

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

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

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

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

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

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

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

孫圣思, 嵇少丞. 2011. 大洋板塊俯沖帶地震波各向異性及剪切波分裂的成因機制. 大地構造與成礦學, 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

國家自然科學基金重大計劃重點支持項目(92058211)、國家自然科學基金面上項目(52074251)、中央高?；究蒲袠I務經費(842012003)和高等學校學科創新引智計劃項目(B20048)聯合資助。

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

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

P541

A

1001-1552(2023)06-1232-010

10.16539/j.ddgzyckx.2023.06.003