俯沖板塊的動力學(xué)性質(zhì)對平板俯沖的影響：數(shù)值模擬

朱致遠(yuǎn)，吳本君

南京大學(xué) 地球科學(xué)與工程學(xué)院，南京 210023

根據(jù)弧后應(yīng)力狀態(tài)，科學(xué)家們將俯沖帶分為馬里亞納型俯沖帶和秘魯—智利型俯沖帶兩類（Uyeda and Kanamori, 1979）。秘魯、智利附近的俯沖帶是秘魯、智利型俯沖帶的典型代表，納茲卡板塊向南美洲板塊下方快速俯沖。在墨西哥中部、秘魯和智利中部，俯沖板塊在南美大陸巖石圈下形成水平的俯沖區(qū)域（Manea et al., 2017）。把俯沖板塊在大陸巖石圈下出現(xiàn)水平片段的俯沖叫做平板俯沖（flat-slab subduction）。

平板俯沖與一系列特殊地質(zhì)現(xiàn)象有著緊密聯(lián)系，往往伴隨著強烈的上覆大陸巖石圈擠壓和地殼增厚，有規(guī)律的巖漿活動（Gutscher et al., 2000a;Kay and Mpodozis, 2001; Sigloch et al., 2008），以及俯沖帶上覆巖石圈的熱異常結(jié)構(gòu) （Dávila and Lithgow-Bertelloni, 2015）。例如，Gutscher等(2000a)認(rèn)為南美洲西部出現(xiàn)的與年輕大洋俯沖相關(guān)的第四紀(jì)埃達(dá)克火山巖可以通過平板俯沖模型進(jìn)行解釋。Li等（2007）使用平板俯沖模型解釋華南中生代的造山運動中淺海盆地的形成和相關(guān)的巖漿活動。Sigloch 等（2008）則通過長距離平板俯沖模型來解釋北美洲西部75 Ma左右大范圍巖漿活動和從海岸到大陸內(nèi)部的長距離造山運動。此外，相較于高角度的俯沖，平板俯沖還會阻礙超高壓巖石折返（Li et al., 2011）。

近年來，對平板俯沖形成機制與影響因素的研究認(rèn)為浮力是影響平板俯沖形成的第一要素（Pilger,1981; Skinner and Clayton, 2011），年輕的板塊與增厚的洋殼會使得俯沖板塊的平均密度較小，使其難以下沉從而使得平板俯沖更易于發(fā)生。此外，特殊的上覆巖石圈熱力學(xué)結(jié)構(gòu) （Rodríguez-González et al., 2012）、上覆巖石圈向海溝的運動（Lallemand et al., 2005）、俯沖板塊的流變性質(zhì)（Van Hunen et al.,2002; Billen and Hirth, 2007; Manea and Gurnis, 2007）、板塊與地幔的耦合作用（Boutelier and Cruden, 2008）等都影響平板俯沖形成。

俯沖板塊的厚度、粘度、密度等動力學(xué)性質(zhì)是決定其應(yīng)力應(yīng)變狀態(tài)的重要因素，很大程度上約束平板俯沖的形成（Huangfu et al., 2016）。然而，俯沖板塊的動力學(xué)性質(zhì)如何影響平板俯沖形成以及形態(tài)的變化仍需要進(jìn)一步的深入研究。下沉的俯沖板塊在上覆的大陸巖石圈下形成平板區(qū)域的長度和角度變化與俯沖板塊的動力學(xué)性質(zhì)的關(guān)系仍不明確。本文主要研究俯沖板塊的動力學(xué)性質(zhì)對于平板俯沖形態(tài)的制約。通過二維地球動力學(xué)數(shù)值模擬，著重探討俯沖板塊的厚度、粘度、密度等因素對俯沖過程中形成的上覆巖石圈下方的平板的長度和角度變化的影響，并結(jié)合智利中部的平板俯沖實例進(jìn)行了討論，從而加深了人們對于平板俯沖的動力學(xué)過程的認(rèn)識。

1 數(shù)值方法和模型設(shè)置

1.1 數(shù)值模擬方法

本文建立的俯沖動力學(xué)的二維數(shù)值模型是基于不可壓縮的斯托克斯流體，其控制方程為無量綱動量守恒方程和質(zhì)量守恒方程：

其中，p是動壓力，?ρ是俯沖板塊密度與地幔密度的差值，g是重力加速度，u是速度。應(yīng)力張量τij與粘度η和應(yīng)變率張量的關(guān)系如下：

其中，

簡單起見，由于模擬的俯沖過程的時間比較短（小于35 Ma），筆者的模型實驗設(shè)計中沒有考慮溫度和相變。使用聯(lián)合拉格朗日粒子法的有限元模擬程序Underworld2（Moresi et al., 2019;Stegman et al., 2006）來求解上述動力學(xué)方程。

俯沖板塊在俯沖過程中，受到的應(yīng)力超過屈服應(yīng)力時，會發(fā)生塑性形變。筆者的模型采用馮·米賽斯屈服準(zhǔn)則，將粘性流變和塑性流變結(jié)合，形成實際的有效粘—塑性流變關(guān)系：

其中，地幔的參考粘度η0= 1×1020Pa·s，屈服應(yīng)力τyield=48 Mpa，是應(yīng)變率張量的第二不變量。

1.2 模型設(shè)置

圖1 平板俯沖數(shù)值模型設(shè)置的示意圖（虛線方框放大顯示了初始俯沖帶的特征）Fig. 1 Setup of the numerical model of flat-slab subduction (Dash-line box shows a three-layer subducting slab and a viscous overriding plate)

如圖1所示，此次研究的模型長度為9000 km、深度為1000 km，被分為720×80個網(wǎng)格。每個網(wǎng)格內(nèi)隨機分布20個粒子，保存了物質(zhì)的密度、粘度等性質(zhì)，并可用于追蹤模型中物質(zhì)的運動歷史。上地幔的粘度為ηum=1×1020 Pa·s，與參考粘度η0相同。模型中，下地幔從660 km深度延伸至1000 km，其粘度ηlm為上地幔粘度的100倍。模型中的左右兩側(cè)邊界是周期性邊界，下邊界是無滑移邊界。上邊界的俯沖板塊和上覆板塊上施加了相對運動的速度。俯沖板塊以3 cm/a速率，上覆板塊以4.5 cm/a速率同時向海溝方向運動。主要的模型參數(shù)見表1。

參考模型中俯沖板塊的長度為4500 km，厚度為70 km，對應(yīng)大洋巖石圈年齡約為35 Ma。俯沖板塊密度比地幔密度多60 kg/cm3，大洋巖石圈在負(fù)浮力的作用下，向大陸巖石圈下方俯沖。開始時俯沖板塊和上覆板塊之間放置了50 km粘度與上地幔一致的薄弱區(qū)域，以便俯沖啟動。上覆板塊的長度為3500 km，厚度為120 km。假設(shè)上覆板塊的密度與地幔密度相同，粘度為上地幔粘度的4000倍。由于俯沖板塊在俯沖時受到強烈的拉張和擠壓作用，將俯沖板塊設(shè)置為三層：上下兩層為粘塑性層，中間為較硬的粘性層（圖1b），從而更好地模擬了板塊的應(yīng)力集中區(qū)域發(fā)生的塑性形變（Lyu et al., 2019; OzBench et al., 2008）。

表1 模型參數(shù)Table 1 Model parameters

2 結(jié)果

通過數(shù)值模擬實驗，此次研究了俯沖板塊的動力學(xué)性質(zhì)對平板俯沖形態(tài)的影響。首先，分析了參考模型中平板俯沖形態(tài)隨時間的變化。在參考模型的基礎(chǔ)上，考察了俯沖板塊厚度、粘度、密度等因素對于海溝運動速度，俯沖板片的平板長度和平板角度等參數(shù)的影響。數(shù)值模型的主要參數(shù)設(shè)置如表2所示。

2.1 參考模型

在參考模型REF_MODEL中，俯沖板塊的厚度為70 km，粘度為上地幔粘度2000倍。俯沖板塊密度比地幔密度大，Δρsp=60 kg/m3。為了更好的描述在上覆板塊下方形成的平板的形態(tài)，筆者定義了一些變量（表1）：vtr、tfs、Lfs、α、β和θ。vtr是海溝運動速度，向左運動為正，表示海溝后撤。tfs為平板俯沖發(fā)生時間，即在俯沖板片上首次出現(xiàn)較為明顯的平板區(qū)域的時間。平板區(qū)域擁有比它兩側(cè)更加平緩的傾角以及一定長度。Lfs是平板區(qū)域的水平投影距離，簡稱平板長度。α是平板與水平面的角度，順時針旋轉(zhuǎn)為正，簡稱平板角度。β和θ分別為俯沖角度和板片前段彎折角度，順時針旋轉(zhuǎn)為正。圖2是定義Lfs、α、β和θ的示意圖。

表2 數(shù)值模型Table 2 Numerical models

圖3顯示了參考模型REF_MODEL中平板俯沖的形成和演化的過程。在推動力與負(fù)浮力共同作用下，俯沖板塊首先向著上覆板塊下方俯沖（圖3b）。在約22 Ma時，如圖3c所示，俯沖板片前段在上地幔底部滑動，俯沖板塊中部出現(xiàn)一段與水平線的角度較小的區(qū)域，認(rèn)為平板俯沖在這時發(fā)生。從圖4a中可以看到，此時海溝后撤的速率vtr較大，約為4.0 cm/a。平板俯沖剛發(fā)生時，平板長度Lfs約為230 km（圖4b）。圖4還顯示了這時平板角度α約為20°，板塊俯沖角度β為35°，θ為44°。平板區(qū)域的角度與俯沖板片的前后段的彎折角度相差約15°～20°。

隨著時間的變化，平板俯沖進(jìn)一步發(fā)展。如圖4中實線所示，平板角度α變得更加平緩，平板區(qū)域的長度Lfs緩慢增加。海溝運動的速率vtr逐漸減小。板片俯沖角度β和板片前段彎折角度θ變化不大，約為5°～6°。圖3d顯示了32 Ma時平板俯沖的形態(tài)。這時，海溝運動速率減小至3.45 cm/a。平板區(qū)域長度Lfs約為245 km，變化不大。平板角度接近0°，與板片前后段有較明顯的角度差異（相差約為40°）。

圖2 模型參數(shù)示意圖Fig. 2 Definition sketch of the flat-slab

圖3 參考模型REF_MODEL的平板俯沖演化過程Fig. 3 The evolution of the reference model REF_MODEL

2.2 板塊厚度的影響

俯沖板塊的厚度直接影響了由于負(fù)浮力導(dǎo)致的板片拉力。筆者測試了不同的板片厚度對平板俯沖形態(tài)的影響（圖4，5）。改變厚度時，俯沖板塊內(nèi)部分層厚度保持比例。模型H_50中，俯沖板塊的厚度較薄，Hsp=50 km，由于板片負(fù)浮力產(chǎn)生的拉力較小，板片俯沖的速度較慢（圖5ac）。而且在上覆板塊的擠壓過程中，薄的俯沖板塊發(fā)生了較大變形，不利于平板俯沖的發(fā)育。板片在約20 Ma時開始發(fā)生平板俯沖，平板區(qū)域長度Lfs只有約180 km，且隨時間變化有較大起伏（圖4b）。由于薄板片易變形，平板的傾斜角度α在32 Ma時接近-10°（圖4c，負(fù)號表示此段板片向上彎折）。圖4e顯示，板片前段的彎折角度θ也較大，在55°～ 64°之間變化。

圖4 俯沖板塊厚度對板片平俯沖的影響Fig. 4 Influence of subducting plate thickness to flat slab subduction

圖5 模型H_50（a-c）和模型H_90（d-f）隨時間的演化過程Fig. 5 Time evolution for the numerical models H_50 (a-c) and H_90 (d-f)

當(dāng)俯沖板塊厚度為90 km時（模型H_90），較厚的板片產(chǎn)生了較大的負(fù)浮力，導(dǎo)致板塊的俯沖和后撤的速度也更快。因此，俯沖過程中沒有形成平板區(qū)域（圖5d-f）。同時，海溝運動的速率變化也不大，約為4 cm/a（圖4a）。如圖4c所示，17～ 33 Ma的時間里，板片傾斜角度α維持在30°左右。板片俯沖角度β隨時間基本沒有變化，與參考模型REF_MODEL相比，差別也不大。隨著板片前段在上地幔底部滑移，板片彎折角度θ從41°變化到30°。

通過以上分析，筆者發(fā)現(xiàn)太薄和太厚的板片都不利于平板區(qū)域的產(chǎn)生。參考模型REF_MODEL中的70 km厚度的俯沖板塊較有助于平板俯沖的發(fā)生。

2.3 板塊粘度的影響

俯沖板塊的粘度直接決定了板片的強度，控制了板片變形的難易程度。筆者考察了俯沖板塊粘度分別為600ηum，1000ηum和4000ηum時，平板俯沖形態(tài)的變化。雖然隨著俯沖過程的演化，所有模型中海溝運動速率都在減小，如圖6a所示，但是俯沖板塊的粘度直接影響了海溝后撤的速率。板片粘度越小，海溝后撤速率越大。當(dāng)ηsp/ηum=600時（模型V_600），在32 Ma時海溝運動速率為3.7 cm/a，比模型V_4000（ηsp/ηum= 4000）中的速率大0.4 cm/a。

此外，從20 Ma到30 Ma的平板俯沖演化過程中，四個模型的平板長度都隨著時間而增加（圖6b）。模型V_4000中，平板長度Lfs從235 km增加到了275 km。而模型V_600中的平板長度隨時間僅有少量變化，從190 km 增加到了200 km。不同粘度的模型中，平板角度α都隨平板俯沖的演化而趨于水平；板片俯沖角度β和前段彎折角度θ隨時間變化不大，約在5°以內(nèi)。

2.4 板塊密度差的影響

根據(jù)半空間冷卻模型，板塊厚度和密度都會隨著年齡增加而增大?？紤]到板片的密度同時受到應(yīng)力、組分等多方面因素的影響，在未考慮溫度場和相變的情況下，將密度和厚度分別進(jìn)行討論。前人的研究表明，俯沖板塊的密度比地幔密度重，差值約為50～100 kg/m3（Skinner and Clayton,2011; OzBench et al., 2008）。 圖7顯示了俯沖板塊與地幔的密度差對平板俯沖演化過程的影響。從圖7a可以看出，如果發(fā)生了平板俯沖，如模型D_50和D_70，俯沖板片的密度差的變化對海溝后撤的速率影響不大。板塊密度差越小，越容易發(fā)生平板俯沖。如圖8a-c所示，當(dāng)Δρsp=50 kg/m3（模型D_50）時，平板俯沖在約15 Ma就發(fā)生了，隨著時間的演化，在32 Ma時，平板長度Lfs達(dá)到313 km。但是，當(dāng)板片密度差很大時，以Δρsp=80 kg/m3（模型D_80）為例，由于負(fù)浮力增大，板片難以被地幔流托起。如圖8d-f所示，俯沖過程中沒有形成平板區(qū)域，板片中段的傾斜角度α在25°～32°之間變化。與模型H_90類似，海溝運動速率的變化也不大，在3.9 cm/a上下浮動。圖7d和7e顯示不同的俯沖板片密度差對板片的俯沖角度β和前段彎折角度θ基本沒有影響。

圖6 俯沖板塊粘度的影響Fig. 6 Influence of subducting plate viscosity

圖7 俯沖板塊密度的影響Fig. 7 Influence of subducting plate density

圖8 模型D_50（a-c）和模型D_80（d-f）隨時間的演化過程Fig. 8 Time evolution for the numerical models D_50 (a-c) and D_80 (d-f)

這些結(jié)果與前人的研究一致（Espurt et al.,2008; Huangfu et al., 2016）。

3 討論

3.1 智利中部的平板俯沖

南美俯沖帶是納茲卡板塊向南美洲板塊西緣俯沖而形成的。俯沖過程的一個構(gòu)造結(jié)果是板間和板內(nèi)地震的發(fā)生。這些地震形成了貝尼奧夫帶，描述了俯沖板片的形狀。通過對南美洲板塊下方的貝尼奧夫帶的研究，科學(xué)家們發(fā)現(xiàn)俯沖板片在秘魯和智利下方形成了數(shù)百米的近水平區(qū)域（Gutscher et al., 2000b）。Manea等（2017）總結(jié)和比較了秘魯和智利平板俯沖帶，認(rèn)為海溝后撤、俯沖板塊的構(gòu)造演化和上覆板塊結(jié)構(gòu)的不連續(xù)性明顯有助于平板俯沖的形成。Chen等（2019）的納茲卡板塊重建模型發(fā)現(xiàn)，納茲卡俯沖在80 Ma前先從南美洲北部（5°S）開始。然后向南傳播，大約55 Ma前，納茲卡俯沖到達(dá)南美洲南部（40°S）。該研究還發(fā)現(xiàn)，前淵沉積和安第斯山脈的形成與納茲卡板片和下地幔之間的相互作用有關(guān)。

在秘魯附近，由于俯沖開始時間較早，俯沖板片已經(jīng)進(jìn)入下地幔并堆疊（Chen et al., 2019）。在智利中部，地震層析成像顯示納茲卡板片還停留在下地幔頂部（圖9），與筆者的參考模型REF_MODEL得到的俯沖板片的形態(tài)相似。參考模型REF_MODEL預(yù)測的平板片長度約為240 km，與前人研究結(jié)果相近（Manea et al., 2017）。智利中部的納茲卡板塊俯沖時的年齡約為38 Ma （Manea et al., 2017），對應(yīng)的板塊厚度約為75 km。由前面的分析可知，這個板片厚度有利于形成平板俯沖。納茲卡板片前段的彎折角度約為48°（Li et al.,2019），參考模型REF_MODEL得到的彎折角度θ與之相一致。本文的數(shù)值模擬實驗可以為研究南美洲智利地區(qū)的平板俯沖過程和機制提供參考。

3.2 板塊運動速度和海溝后撤

皇甫鵬鵬等（2016）的研究發(fā)現(xiàn)俯沖大洋板塊和上覆大陸板塊相對運動的速率對平板俯沖的形成有著決定性的影響。上覆板塊朝向海溝的運動速率越大越有利于形成平板俯沖。Coltice等（2017）根據(jù)HS3-NUVEL1A模型（Gripp and Gordon, 2002）的結(jié)果，統(tǒng)計了全球主要俯沖帶兩側(cè)板塊的運動速度和海溝的運動速度。其中，在智利附近的俯沖帶，俯沖的納茲卡板塊的運動速率約為3.3 cm/a，而上覆板塊的運動速率約為4.6 cm/a。本文數(shù)值模型所用的板塊運動速度邊界條件也是根據(jù)以上結(jié)論設(shè)置的。

根據(jù)Li等（2019）的總結(jié)，海溝后撤的模型中俯沖板塊與地幔轉(zhuǎn)換帶呈較小角度，最終導(dǎo)致板片滯留在地幔轉(zhuǎn)換帶處。本文的模型設(shè)置中，上覆板塊的向海溝逆沖速率大于俯沖板塊的運動速率，海溝后撤的速率約為3～4 cm/a，與前人的研究結(jié)果一致（Coltice et al., 2017）。值得注意的是，模型H_90和模型D_80沒有出現(xiàn)水平板片，它們的海溝運動速率在3.9～4.0 cm/a 之間浮動。其他發(fā)生平板俯沖的模型中，海溝運動的速率都隨著平板片的形成而減小，有的模型（V_4000）的海溝后撤速度減小到3.2 cm/a。

本文主要關(guān)注了俯沖板塊的厚度、粘度和密度差等性質(zhì)分別對平板片形態(tài)的影響。簡單起見，溫度相變等因素并沒有被考慮。上覆板塊的厚度粘度等性質(zhì)也有可能對平板俯沖的動力學(xué)過程有顯著影響（Manea et al., 2012; Huangfu et al., 2016），筆者將在后續(xù)的工作中進(jìn)行更加深入的研究。

4 結(jié)論

本文研究了平板俯沖的動力學(xué)特征，以及俯沖板塊的動力學(xué)性質(zhì)對俯沖平板片形成的影響。通過數(shù)值模擬實驗，得出了以下初步結(jié)論：

（1） 俯沖板塊的厚度和密度差對平板片的形成有著決定性的影響。厚度為70 km左右的板片最容易發(fā)生平板俯沖。太厚的俯沖板塊（Hsp=90 km）比較難發(fā)生變形，阻礙平板俯沖的發(fā)生。俯沖板塊與上地幔的密度差越小，越容易發(fā)生平板俯沖。但是，當(dāng)板片密度差很大時（Δρsp=80 kg/m3），俯沖過程中不能形成平板區(qū)域。

圖9 智利中部的板片平俯沖Fig. 9 Central Chile flat slab

（2） 俯沖板塊的粘度也是平板片形態(tài)的主要影響因素。較大的板片強度有助于形成平板俯沖。板片粘度越大，形成的平板片長度越長。

（3）平板俯沖的過程伴隨著海溝后撤速率的減小。