俞霄，王世民

(中國科學院大學地球與行星科學學院 中國科學院計算地球動力學重點實驗室， 北京 100049) (2020年4月28日收稿； 2020年6月2日收修改稿)

在經典的板塊構造理論中，巖石圈板塊是剛性的，板塊在地球表面的運動可近似看作是剛體的定點轉動，由角速度向量(Euler vector)描述[1]。描述板塊運動的模型分為板塊相對運動模型和板塊絕對運動模型。板塊相對運動模型基于地球物理和空間大地測量數據(包括地震滑移矢量、洋中脊擴張速率、轉換斷層方位角及GPS等空間觀測的結果)，確定板塊的劃分并得到板塊之間相對運動的角速度[2-6]。板塊絕對運動模型研究板塊相對于深部地幔的運動。前人研究中常用的代表深部地幔的參考系包括熱點參考系[7-9]和無合轉動(no-net-rotation)參考系[5,10]。由于板塊邊界受力的不對稱及軟流圈黏度的不均勻性，巖石圈可能存在整體的合轉動[11-13]，基于熱點參考系的板塊絕對運動模型更符合板塊動力學研究的需要。本文基于相對起源于深部地幔柱的熱點定義的板塊絕對運動模型T25M[13]研究板塊運動的驅動力。

前人關于剛性板塊運動的研究通常以下列兩個基本假設為前提：一是板塊運動的角速度在至少幾百萬年內是保持恒定的，可以給出一套板塊的角速度向量描述相當長時間內各板塊的運動[4]；二是各板塊處于力矩平衡狀態(tài)，可以用板塊力矩平衡擬合得出板塊受到的各種力矩的相對大小[14]。事實上，以上兩個假設并不等價。當剛體旋轉軸不沿其慣量主軸時，其角速度與角動量方向存在夾角，由剛性板塊的力矩平衡可以導出其角動量恒定，但不能同時保證角速度恒定。本文將對這一問題作深入分析。

對于板塊運動的驅動力，早期的研究認為地幔對流推動著板塊運動，就像在傳送帶上運送物品一樣[15]。然而，考慮到巖石圈作為邊界層成為地幔對流系統的一部分時，板塊底面的剪切力總體上起著黏滯板塊運動的阻力作用[16-17]。后續(xù)大量研究工作將巖石圈處理為獨立的子系統，研究該子系統中的各種作用力[18-19]，主要包括：與俯沖板片有關的力、洋中脊推力、板底拖曳力和碰撞帶產生的力等。

有關板塊運動的驅動機制，現有的研究主要從以下兩方面開展：第1類研究從板塊力矩平衡出發(fā)，設定一些力矩的參數，擬合求解各種力矩在力矩平衡條件下的取值[14,18-21]。在計算板塊受到的力矩時，通常以地心作為矩心(支點)，其優(yōu)勢是可以近似認為板塊受到的重力指向地心而軟流圈對板塊的支持力背離地心，這兩種力的力臂為零，從而簡化力矩平衡的計算過程。這類研究的主要結論是：俯沖板片負浮力與板片所受阻力比其他作用力約大一個數量級，但俯沖板片上的合力與其他驅動力量級相當；地幔對巖石圈底部的剪切力阻礙板塊運動，且大陸板塊板底剪切力比大洋板塊更大。由于問題的超定性，雖然不同作者假設的驅動力模式設定有所不同，但均能為各板塊建立較好的力矩平衡[22]。第2類研究直接根據板塊邊界形狀計算板塊所受驅動力矩，比較驅動力矩與板塊絕對運動角速度的方向[23-24]。由角動量定理可知，力矩與角動量而不是角速度之間存在相關性，將驅動力矩與板塊角速度直接對比缺乏足夠的物理基礎。

在計算板塊邊界力矩方面，源于大洋板塊地形的洋脊推力的計算方法較為成熟，而由于俯沖板片所受地幔阻力不易計算，俯沖板片凈拉力的計算存在困難[25]。在計算板底剪切力矩方面，前人普遍采用均勻拖曳系數的方法[14,18-19]，僅能考慮大陸和大洋下軟流圈的黏度差異。

本文采用板塊驅動上地幔流動的模型，并利用由上地幔物理條件約束的地幔黏度分布[25]，計算上地幔物質的剪切流動狀態(tài)，進而計算作用于各板塊的板底剪切力及其力矩。結合洋脊推力矩的定量計算結果，研究現今巖石圈板塊的力矩平衡關系，進而對板底剪切和洋脊推力以外的板塊作用力進行討論。

1 板塊運動穩(wěn)定性分析

考慮板塊運動為剛體的定點轉動，角動量定理要求板塊所受合外力矩T等于角動量L對時間的導數

(1)

其中角動量L定義為轉動慣量I與角速度ω的乘積

L=Iω.

(2)

剛體的轉動慣量I是一個二階張量，在直角坐標系中可以按以下形式積分計算

I=

(3)

如果建立在固結于剛體慣量主軸的坐標系，其中的剛體轉動動力學歐拉方程可以寫為以下形式

(4)

其中:Ixx、Iyy、Izz為3個方向的主慣量，ωx、ωy、ωz為3方向角速度分量，Tx、Ty、Tz為3方向合力矩分量[26]。

如果板塊受到方向與運動速度相反的板底阻力約束，考慮簡單的板塊驅動上地幔層流，板底阻力矩Tpb為

(5)

其中:r為地球半徑矢量，μ為上地幔黏度，A為板塊底面積，hum為上地幔厚度。

對于均勻密度ρ、厚度hlp的巖石圈板塊，其主慣量為

I=ρr2Ahlp,

(6)

結合式(2)、式(5)和式(6)可以得到

(7)

式(7)表明，板底阻力矩與板塊角動量方向相反，大小成正比，比例系數為

(8)

取上地幔黏度μ量級1021Pa·s、巖石圈密度ρ量級3×103kg/m3、巖石圈厚度hlp量級100 km和上地幔厚度hum量級600 km的典型值計算[25]，比例系數c量級為106/s。

帶有板底力矩的歐拉方程如下

(9)

這里Tx、Ty、Tz代表除板底力矩以外的合力矩分量。由數量級估計，板塊轉動慣量數量級為1035kg·m2，角速度數量級為10-16rad/s，所受力矩數量級為1025N·m。由此，歐拉方程中第2項遠小于其他項，可以忽略。方程(9)近似為

(10)

解得

(11)

其中：ωx0、ωy0、ωz0為初始角速度的3個分量。式(11)解中右側第2項在10-6秒量級時間內迅速衰減，即由于板底剪切力矩對于板塊運動的約束非常強大，對于板塊邊界上施加力矩的變化，板塊角速度的響應十分迅速，能在很短的時間內達到新的平衡。

如果認為一般情況下，板塊處于力矩平衡狀態(tài)，其角動量保持恒定，板底阻力矩近似與板塊角動量方向相反，則其他力矩的合力矩與角動量方向近似相同。這種板塊角動量方向與驅動力矩之間的內在聯系并不存在于角速度方向與驅動力矩之間。由于板塊在運動過程中轉動慣量可以不斷變化，引起其角速度與角動量間夾角的變化，角速度方向與驅動力矩間沒有簡單的依賴關系。因此，不同于Richardson[24]比較力矩與角速度方向的方法，本文提出更符合物理規(guī)律的比較板塊角動量與板塊受到的各類力矩以識別驅動力的方法：驅動板塊運動的力矩方向應接近板塊角動量的方向，而阻礙板塊運動的力矩則大致與角動量方向相反。這樣可以有效識別板塊運動的驅動力與阻力。

2 板塊角動量計算

計算板塊的角動量需已知巖石圈板塊的幾何劃分、密度分布和運動狀態(tài)。NNR-MORVEL56模型[10]給出了全球56個板塊的邊界劃分，這一模型在PB2002模型[27]的基礎上新增4條板塊邊界，PB2002模型給出了全球各條板塊邊界的類型。最新的基于深部地幔柱參考系的板塊絕對運動模型T25 M[13]給出了56個板塊絕對運動的角速度向量。Litho1.0巖石圈結構模型[28]給出了全球巖石圈的三維結構信息，該模型將地球表面劃分為40 962個點組成的81 920個三角形的網格，給出了每個坐標點下巖石圈的深度分層與密度信息。為簡化運算，本文將模型中巖石圈的9層劃分簡化為4層并對每一層的密度取一個確定值(圖1)。

選取全球面積最大的19個板塊作為研究對象(占地球表面積的97%)，通過尋找每個板塊內的Litho1.0模型網格，將其質量分布做式(3)中的積分，得到每個板塊的轉動慣量張量，用慣量張量乘以T25 M模型給出的板塊角速度可以得到板塊角動量向量。計算得到的各板塊的角動量向量見表1。

圖1 本文簡化后的Litho1.0模型[28]分層示意Fig.1 Schematics of the Litho 1.0 model[28] simplified in this paper

表1 板塊角動量Table 1 Angular momentum of the plates

3 板塊作用力矩計算

3.1 板塊邊界力矩計算

洋脊推力由洋底地形與大洋巖石圈的密度結構決定。Turcotte和Schubert[25]得出，洋脊推力隨大洋巖石圈年齡成正比地增大，直到年齡100 Ma的大洋板塊洋脊推力達到約為4×1012N/m的典型值后，更老的洋底基本不再繼續(xù)下沉。本文結合大洋巖石圈年齡分布[29]對老于100 Ma的大洋板塊均設定洋脊推力為4×1012N/m，而對較為年輕的大洋板塊根據其年齡做比例調整，計算得到洋脊推力(圖2)。采用MORVEL56模型板塊邊界類型的劃分，以地球半徑作為力臂，對每一段洋脊的推力求得力矩，對這些力矩求和可以得到一個板塊受到的洋脊合力矩。

俯沖板片因其熱結構受到負浮力，因俯沖運動受到周圍地幔物質的流體壓力，表現為板片下方的托力和上方的吸力[30-31]。負浮力垂直于板片的分量與板片下方地幔的托力和上方地幔的吸力平衡，維持了俯沖角度[32]。負浮力沿板片分量與周圍地幔對板片阻力的合力作為板片凈拉力作用在板塊上[33-34]。板片負浮力比洋脊推力大一個量級，但板片受到周圍地幔阻力與地幔流動和黏度狀態(tài)相關，使得俯沖板片的凈拉力可能與洋脊推力量級相同[25]。Schellart[34]通過沙箱實驗方法得到俯沖板片凈拉力大約是負浮力的8%～12%。俯沖板片負浮力沿板片方向分量的大小可以由以下公式[14]計算：

Fsp=ΔρgLHsinθ,

(12)

其中:Δρ為俯沖板片與周圍地幔物質的密度差，McNutt[35]給出的參考值為80 kg/m3，g為重力加速度,L為俯沖板片長度，H為俯沖板片厚度，θ為俯沖傾角?？紤]巖石圈冷卻模型，H可由以下公式[25]計算：

(13)

其中:熱擴散系數κ取值1 mm2/s,t為俯沖處的板塊年齡。Lallemand等[36]整理了全球俯沖帶上159個地點的俯沖板片的各項參數，包含板片俯沖處的大洋巖石圈年齡、板片長度和俯沖傾角。由其文中的參數和上述公式可以計算俯沖板片負浮力沿板片方向的分量，本文假設該力在俯沖板片彎折處直接傳遞作為對板塊的水平拉力。將159個點的取值在同一條板塊俯沖邊界上做平均處理，得到全球各條俯沖邊界上單位長度負浮力對板塊水平拉力的大小(圖2)。以地球半徑為力臂，估算出各板塊的負浮力力矩。

紅色邊界為洋中脊，藍色邊界為俯沖帶，板塊簡稱參照MORVEL56模型[10]圖2 洋脊推力與俯沖板片負浮力(單位1012 N/m)的全球分布Fig.2 The global distribution of the ridge push and slab pull forces measured in 1012 N/m

為簡化計算，假設全球所有轉換斷層邊界上不受力，大陸裂谷受到1012N/m的擴張力，碰撞帶邊界受到1012N/m的擠壓力。以地球半徑為力臂計算裂谷和碰撞帶上的力矩，這兩種力矩的大小只作為數量級上的參考。

3.2 板底剪切力矩計算

定義板底剪切力為板塊運動的阻力，其方向與當地的板底速度方向相反。由牛頓黏滯定律，對于板底面積元Ai，其板底力如下

(14)

應用干上地幔(dry upper mantle)位錯蠕變下的黏度隨溫度、壓力變化的計算公式[25，37]：

(15)

(16)

上地幔黏度可以表達為

(17)

其中:a為上地幔物質的活化能Ea、地表的熔點Tm0和普適氣體常數R決定的常數，取值為29.4，Tm是上地幔物質的熔點。

對于上地幔溫度分布，設定660 km深度溫度為1 875 K，上地幔溫度梯度為絕熱溫度梯度[25]：

(18)

其中:αv=3×10-5K-1，cp=1 000 J·kg-1·K-1，g=9.8 m·s-2，可以得到溫度T(K)與深度h(km)的關系

T=1 544.3e2.94×10-4 h.

(19)

對于上地幔物質的熔點Tm，Turcotte和Schubert[25]給出的參考值為地表Tm0=2 140 K，Tm隨深度增加的合理梯度為2 K/km，這個梯度取值較大。采用Tm取2 K/km和1.5 K/km兩個深度梯度，分別對于一個厚度為100 km和250 km的巖石圈板塊，令其運動速度u=1.6×10-9m/s (5 cm/a),由式(17)的黏度分布計算板底剪應力與上地幔物質流速分布，如圖3所示。

由圖3可以看出，Tm梯度取2 K/km時得到的上地幔黏度較大，在660 km深超過了1023Pa·s量級，使得地幔流動只在淺部發(fā)生，其計算出的板底剪應力較大。因此，取Tm梯度為1.5 K/km，即

Tm=2 140+1.5h.

(20)

至此，可以采用式(17)計算全球各地的板底剪應力和各板塊的板底剪切阻力矩。

將計算的全球各地板底剪應力作分布圖，與Litho1.0模型給出的LAB界面深度分布圖對比，如圖4所示。

對比圖4可以看出，對于同一個板塊，板底剪應力大小與巖石圈厚度呈高度正相關。同時，板底剪應力大小也與板塊運動速度呈正相關。AU、CP、IN、PS、PA、NZ、CO等板塊運動速度大，在巖石圈厚度相同的區(qū)域板底剪應力整體大于其余板塊。

3.3 作用于板塊的力矩與驅動力識別

依據前文所述方法分別計算了每個板塊受到的洋脊、板片負浮力、裂谷、碰撞、板底剪切5種力矩，表2列出了9個板塊所受主要力矩的計算結果。

5種力矩中，洋脊力矩和板底力矩的計算過程有較好的物理模型作為約束，其計算結果較為可靠，而負浮力力矩計算結果應當比板片凈拉力矩大一個量級且約束較差。裂谷力矩和碰撞力矩大小的計算結果只具有數量級上的參考價值。

將各種力矩與板塊運動的角動量作夾角，結果列于表3。

由表3可以識別出板塊運動的驅動力。表中所列板塊的板底力矩與板塊角動量的夾角均相當接近于180°，這與前文理論推導得到的板底力矩近似與板塊角動量方向相反的結論吻合?？傮w上，洋脊、負浮力力矩多為驅動力矩，非洲裂谷力矩阻礙了NU板塊的運動，碰撞力矩驅動了EU板塊的運動。

圖4 全球LAB深度[28]與板底剪應力分布對比Fig.4 Comparison between the global distributions of LAB depths[28] and basal shear stresses

表2 作用于板塊的各種力矩Table 2 Various types of torques acting on the plates 1025 N·m

表3 板塊作用力矩與角動量間夾角Table 3 Angles between various torques and the plate angular momentum (°)

4 板塊運動驅動機制討論

4.1 與洋脊力矩和板底力矩平衡的剩余力矩

考慮板塊力矩平衡，以上求得的同一板塊各種力矩之和應為零。如果計算與洋脊力矩和板底力矩平衡的剩余力矩，這一剩余力矩應當可以指示俯沖、裂谷、碰撞等力矩的綜合作用效果。為方便作圖，采用式(21)將力矩化為等效應力

(21)

其中：T為力矩，r為地球半徑，A為板塊面積。圖5選取了板塊中的一些代表點，給出了各板塊上與洋脊力矩、板底力矩和剩余力矩對應的等效應力分布，并區(qū)分了全球板塊邊界的類型。圖中紅線是洋中脊，綠線是轉換斷層，黑線是俯沖帶，紫線是碰撞帶，橙線是大陸裂谷，棕線是MORVEL56模型新增的4條板塊邊界。

由圖5可見，PA、AU、PS這3個板塊剩余力矩的等效應力都指向俯沖邊界方向，其剩余力矩與前文計算的負浮力力矩夾角較小，可將剩余力矩歸因為俯沖板片產生的拉力矩。剩余力矩的等效應力大小和洋脊力相當，支持了前人俯沖板片上的合力與洋脊力量級相當的觀點[18]。

CO、NZ兩個板塊的剩余力矩的等效應力指向俯沖邊界的反向。對此結果的一種解釋為前文計算高估了年輕板塊洋脊力的大小或者低估了年輕大洋巖石圈板底剪切力的大??；另一種解釋為由于東太平洋俯沖帶的俯沖形態(tài)較平緩[38]，存在的碰撞阻力大于俯沖板片拉力，使得俯沖邊界上的力沿板塊運動方向的反向，阻礙板塊運動。

對于NA板塊，洋脊力矩與板底力矩幾乎平衡，即該板塊其他邊界上作用力的合效應很小。對于SA板塊，剩余力矩可以由西側的NZ板塊俯沖帶對上覆板塊的碰撞力解釋，這一碰撞力的大小約為洋脊力的一半。對于NU板塊，剩余力矩的等效應力指向背離裂谷的方向，剩余力矩與前文計算的裂谷力矩方向基本一致，可由略小于洋脊推力的東非裂谷擴張力解釋。

對于EU板塊，剩余力矩的等效應力指向背離地中?！柴R拉雅碰撞帶的方向，剩余力矩與前文計算的碰撞力矩方向基本一致，可由與洋脊推力相當的大陸碰撞帶產生的力解釋。Warners-Ruckstuhl等[39-40]的研究認為喜馬拉雅碰撞帶上的力大小為(7.0～10.5)×1012N/m，中東地區(qū)為(1.3～2.7)×1012N/m。由此，EU板塊可由大陸碰撞力驅動。

對于OK、AM板塊，其剩余力矩可由西側的裂谷擴張和碰撞帶產生的力解釋，這種力可以視為喜馬拉雅碰撞帶產生的力經過EU板塊傳遞后施加在OK、AM板塊上，驅動了板塊的運動。

圖5 由等效應力表示的各板塊洋脊力矩(紅色箭頭)、板底力矩(綠色箭頭)和剩余力矩(藍色箭頭)比較Fig.5 Comparison between the ridge, basal and residual torques expressed by the equivalent stresses acting on the plates (The red, green and blue arrows show the ridge, basal and residual torques, respectively)

AN板塊的運動角速度很小，板底力矩的作用較小，而其四周邊界三邊為洋脊擴張，在SA、NZ板塊南部一邊為轉換斷層邊界，其剩余力矩可由轉換斷層上的碰撞力解釋。

4.2 有關古板塊角速度突然變化的討論

本文式(11)解可以合理解釋地質歷史上板塊運動的突然轉向，即板塊邊界力矩的變化導致板塊角速度變化的過程十分迅速。有關古板塊運動方向的轉變，可以討論不平衡過程前后的兩個平衡過程。用本文采用的力矩平衡方法，如果已知一個古板塊的幾何結構和絕對運動角速度，建立當時的板塊運動力矩平衡模型，比較板塊運動角速度明顯變化前后的兩個力矩平衡過程[41]，可以探究板塊角速度變化的原因。參考前人的古板塊運動模型重建結果[42-44]，對于近50 Ma來全球主要板塊運動特征的變化，有3個板塊值得討論：

1)太平洋板塊在大約43 Ma前運動方向有一次轉向，這次轉向前后的太平洋板塊運動由Hawaiian-Emperor島鏈的彎折真實地記錄下來[45]。對比轉向前后的太平洋板塊邊界可以得出，這次轉向可能與西太平洋俯沖帶的產生有關。如果建立轉向前洋脊驅動力與板底剪切阻力平衡的北向運動模型，和轉向后板塊由西向俯沖帶驅動與板底剪切阻力平衡的西向運動模型[46]，可以驗證這一猜想。

2)印度板塊在大約40 Ma前由于與歐亞板塊的碰撞，運動角速度大幅減小[41]。對比碰撞前后的板塊邊界可以看出，之前向北俯沖的俯沖邊界轉變?yōu)榇箨懪鲎矌?，俯沖拉力轉變?yōu)榕鲎沧枇?。建立碰撞前的洋中脊和俯沖帶驅動，與板底阻力平衡的快速運動模型，和碰撞后的洋中脊和俯沖帶驅動，與板底剪切阻力和碰撞阻力平衡的慢速運動模型，可以說明印度板塊角速度減慢的原因。

3)歐亞板塊在近40 Ma以來在印度、阿拉伯、地中海地區(qū)的相繼碰撞下運動速度有一次轉向，由之前的東南向運動轉變?yōu)闁|北向運動。建立碰撞前的洋脊推力驅動與板底剪切阻力平衡的模型，和碰撞后大陸碰撞力驅動與洋脊推力和板底剪切阻力平衡的模型，可以說明歐亞板塊這次轉向的原因。

5 結論

本文基于板塊的力矩平衡研究其動力學問題，提出由力矩與角動量的夾角大小識別板塊驅動力的方法，并依據最新的板塊絕對運動模型計算了板塊轉動慣量、角動量、板底剪切阻力及板塊邊界的主要作用力?；窘Y論如下：

1)板底剪切力矩與板塊角動量方向大致相反，而板塊受到的其他力矩之和則作為驅動力矩近似沿板塊角動量方向。

2)利用板塊所受力矩與板塊角動量的夾角可以識別板塊驅動力矩?？傮w上，洋脊力矩、板片負浮力力矩多為板塊運動驅動力矩。非洲裂谷力矩阻礙NU板塊的運動，碰撞力矩驅動EU板塊運動。

3)通過建立合理的物理模型計算洋脊力矩和板底力矩，可以求出與之平衡的剩余力矩。剩余力矩可以指示俯沖、裂谷、碰撞等力矩的綜合作用效果，從而解釋各板塊運動的驅動機制。

4)板底剪切力作為一個強黏性約束維持了板塊運動的穩(wěn)定性。由于板底剪切力的數量級很大，板塊邊界力矩變化時，板塊力矩不平衡的調整過程十分迅速?？梢酝ㄟ^建立變化前后兩個力矩平衡模型進行比較的方法探究古板塊運動狀態(tài)明顯變化的原因。

本文對板塊力矩的計算建立在簡單物理模型之上，計算中參數的選取具有諸多的不確定性，巖石圈密度分層也做了一定的簡化。我們將在未來的工作中繼續(xù)完善模型和數據，更深入細致地研究板塊驅動力這一重要的地球動力學課題。