1 引 言

地幔黏度存在剧烈的径向变化^{[1, 2]}和强烈的横向变化^[3].地幔横向黏度变化虽然不能对地幔对流形态起控制性作用^[4~6]，但会对对流状态和格局、相变动力学、地球内部热结构、板块动力学以及地球物理观测场等产生显著的影响^[5~15]，同时还可以在均匀边界条件下获得环型场分量^{[5, 16]}.朱涛等^[17]曾对小横向黏度变化的情形进行过深入的研究，获得的地幔浅部极型场和环型场对流形态对认识和理解现今地表构造的运动状态提供了帮助.然而，由于小横向黏度变化的假设，在朱涛等^[17]的求解过程中并没有实行交叉迭代，使得地幔对流的极型场和环型场的耦合只停留在理论意义上.针对这种情形，本文将地幔横向黏度变化提高到3 个量级，在求解过程中实施了交叉迭代，实现了地幔对流的极型场和环型场的相互耦合.依此获得的地幔浅部的极型场和环型场对流形态得到了显著的改善，能更合理地解释太平洋板块向西或西北、大洋洲向北或北北东、南美洲向西的运动状态以及北太平洋板块的右旋运动.同时揭示出东太平洋洋中脊与极型场对流的发散区对应更好，这为更深入地认识和理解现今地表构造运动的深部动力学过程提供了有益的帮助.

2 基本公式

考虑内径为R₁，外径为R₂ 充满了不可压缩的牛顿流体的同心球层，纯粹从球层底部加热，采用Boussinesq近似，并取流体的普朗特数无限大(Pr=ν/κ，ν为运动学黏度；κ 为热扩散系数.地球地幔的普朗特数约1024^[18])，球层中流体的黏度为流体平均黏度η0 加上横向黏度变化\[\tilde \eta \]，则在球坐标系(r，θ，φ)中控制变黏度流体运动的无量纲方程组见文献[4]中的式(1).将速度场u分解为极型标量场Φ 和环型标量场ψ^[19]，对动量方程(文献[4]中的式(1b))两边分别施加算子(r· Δ×)和(Λ· Δ×)，可得到关于环型场ψ 和极型场Φ 表达式

(1a)

(1b)

其中，

(1c)

(1d)

(1e)

其中，${{\dot{\varepsilon }}_{ij}}$为应变率，$\tilde{T}$ 为扰动温度，R_a为瑞利数.式(1)和热传输方程(文献[4]中的式(1c))共同组成了求解大横向黏度变化下地幔对流模型的基本方程组.可以看出，式(1a)和(1b)比小横向黏度变化情况下^[16]均多了一项1/(1+\[\tilde \eta \])，其他项均相同.在求解小横向黏度变化下的地幔对流模型时，忽略了式(1a)中因横向黏度变化(右端第二项)而产生的极型场成分^[16]，没有实施交叉迭代来求解极型和环型标量场，即极型场与环型场并没有实现真正意义上的耦合.在本文中，\[\tilde \eta \] 的变化会达到3个量级，因此在求解极型场时不能忽略式(1a)右端第二项，在求解过程中需进行交叉迭代式(1a)和(1b)，从而实现极型场与环型场真正意义上的耦合.显然，这种耦合肯定会调整极型场和环型场的对流格局，使得本文的结果比小横向黏度变化下的结果^[17]更接近真实情形，为更深入理解和认识横向黏度变化对地幔动力学的影响提供有益的帮助.

在此跟小横向黏度变化下的情形一样^[17]，研究了上边界为刚性或应力自由，下边界为应力自由等温边界的情形，其数学形式为：

上、下均为应力自由等温边界(F-F边界)

(2)

上边界为刚性、下边界为应力自由等温边界(R-F边界)

(3a)

(3b)

其中，r₁ 和r₂ 分别为球层内径R₁、外径R₂ 的无量纲值.

3 横向黏度变化形式

室内实验表明，地幔的黏度强烈依赖于温度、压力和应力(或应变率)^[20].有效黏度可以表述为

(4)

其中E^* ，V^* 分别为激活能和体积，F为常量，R是气体常量，T为温度，k是幂率指数.在文中，采用了Zhang和Christensen^[21]研究中的黏度形式，

(5)

其中H(r)是与随深度变化相关的压力、平均温度、平均应变率和矿物结构等引起的径向分层黏度.δρ为密度异常.ω 是小量，是为了防止$\dot{\varepsilon }\to $0时而导致黏度变为无穷大.常量c描述在黏度随温度线性变化和黏度随温度指数变化的假设下黏度的变化.对于一次近似，c可以根据式(6)计算

(6)

其中T 为平均温度.令式(6)中的幂率指数k=1时，则获得随温度变化的黏度，此时为牛顿黏度.所以式(6)变为

(7)

在此，用地幔平均黏度η0 来替代分层黏度H(r)，便可获得横向黏度变化的分布

(8)

显然，对于等黏度地幔，常量c=0.0，即其最小值为0.0，而其最大值据Zhang和Christensen估计可能不超过1.0 m³/kg^[21].用地幔平均黏度η0 对式(8)进行无量纲化后可得

(9)

本文的横向黏度变化形式与朱涛等^[17]研究小横向黏度变化下的地幔对流相同.可通过式(10)由地震波速异常而获得.然后通过式(9)可获得横向黏度变化$\tilde{\eta }$，它具有三维分布形式.

假定地幔中地震波速异常对应的地幔横向不均匀结构是地幔热对流的结果，则可以将地震波速异常通过式(10)转换为地幔密度异常^[22].

(10a)

(10b)

其中，ρ 为地幔密度.v_s、v_p 分别为地震横波和纵波速度.

4 球函数展开

将极型标量场Φ、环型标量场ψ、扰动温度场珦T和横向黏度变化\[\tilde \eta \] 展成复球谐形式，分别如下

(11a)

(11b)

(11c)

(11d)

其中，Y_l^m(θ，φ)是正规化的l阶m级面谐函数，

(12)

(13)

其中，P_l^m(cosθ)是正规化的l阶m级的缔合勒让德函数.面谐函数的正交归一关系为

(14)

其中，Y_l^*m是Y_l^m的共轭.将式(1a)和(1b)展开为球谐形式为

当t≠０时，

(15a)

(15b)

当l=0，且l₁ =l₂ =l₃ =0时，

(16a)

(16b)

当l=0，且l₁，l₂，l₃ 不同时为零时，

(17a)

(17b)

其中，

积分系数H₁ ~ H₉₇ 的表达式见附录A.对于式(17)，其系数H₈₅ ~ H₉₇ 比附录A 中的减少一项，其他项均相同.

将热传输方程(文献[4]中的式(1c))展开成球谐形式为

(18)

其中，T_c 是无对流运动时的温度场，仅是径向r的函数，满足热传导方程Δ2T_c=0.积分系数Γi(i=1，2，3，4，5)的表达式见附录B.将边界条件(3)展开成球谐形式为

(19a)

(19b)

5 Galerkin方法

Galerkin方法是将连续问题转换为离散问题的一种方法.Zebib 等^[23]、叶正仁等^[24]以及朱涛和冯锐^[25]将之用来求解地幔对流模型.在此，也采用Galerkin方法来进行求解.设

(20a)

(20b)

(20c)

取基函数

(21)

则Θ_l^mj自然满足关于温度的边界条件.将式(20)和(21)代入式(15)~(19)，并利用三角函数的正交性可得

当l≠0时，

(22a)

(22b)

当t＝０，且t１ ＝t２ ＝t３ ＝０时，

(23a)

(23b)

当l=0，且l₁，l₂，l₃ 不同时为零时，

(24a)

(24b)

对于F-F边界条件，

(25)

(26a)

(26b)

对于热传输方程式(18)，

(27)

其中，

(28a)

(28b)

(28c)

由此通过式(22)~(24)和边界条件(25)或(26)求解出f_lj^m和h_lj^m，然后通过式(27)求解出Θ_lj^m，最后利用式(20)求出珦T_l^m、Φ_l^m、ψ_l^m，便可获得对流图像.

6 结果及分析

为了与小横向黏度变化下的结果对比，在文中采用了与文献[17]相同的参数，即球谐展开的最大阶数取l_max= 6，同时取j_max=k_max=6.R₁=3471km，R₂=6271km，它们分别对应于核幔边界和岩石圈底界面的位置.瑞利数R_a=106.地球地幔的地震波速异常采用SH12WM13模型求得^[25].沿经度、纬度和深度方向上采用144×72×56 的网格(相当于平面上2.5°×2.5°的网格，深度上每层厚度为50km)，利用PREM 地球模型^[26]以及地震波速异常和密度之间的转换关系(10)计算地幔密度异常，然后利用式(9)获得地幔内部的横向黏度变化.一般情况下，横向黏度变化的最大值出现在上地幔和下地幔中，达到了2~4个量级，而中地幔则为1~2个量级^[3]，因此文中在将密度异常转换为横向黏度变化时，通过选择合适的常量c，将横向黏度的变化控制在3个量级内.

由于上、下地幔密度异常最大值是不同的，因此在将密度异常转换成横向黏度变化时，需选择不同的常量c.在文中，对于上地幔，利用SH12WM13模型^[24]求得的正、负密度异常绝对值的最大值均出现在约100km 深度处，分别约为83.11kg/m³ 和58.63kg/m³，用式(9)可获得常量c分别约为0.083m³/kg和0.118m³/kg;而对于下地幔，正、负密度异常绝对值的最大值分别出现在约2550km和2700km 深度处，约为38.10kg/m³ 和55.63kg/m³，则常量c分别约为0.181m³/kg和0.124m³/kg.图 1显示了200km、300km、400km 深度处横向黏度变化的分布图像.图 1揭示出，大陆下的横向黏度变化基本为正值，而大洋下，除了200km 深度处太平洋中部的横向黏度变化为正值外，其他区域基本为负值；在上地幔，横向黏度变化剧烈程度随深度而降低，且变化最剧烈的区域也随深度而不同：对大陆而言，200km 深度处横向黏度变化最剧烈的区域位于北美大陆下，最大横向黏度变化约309.0(无量纲值);300km 深度处出现在欧亚大陆下，最大横向黏度变化约137.5;400km 深度处则出现在欧亚、澳洲及非洲北部大陆下，最大横向黏度变化约27.1.对海洋而言，200km、300km 和400km 深度处均出现在东太平洋区域，最大横向黏度变化绝对值分别约240.7、48.6 和34.0.文中提供的不同深度的横向黏度变化(图 1)的这些特征是由地幔中密度异常的分布来决定的(见式(9)).

6.1 数值方法及程序验证

在数值离散过程中，对一阶和三阶偏微分采用向前差分格式，对二阶和四阶则采用向后差分格式.对非线性方程组采用Levenberg-Marquardt 方法^[27]进行求解.两种方式验证了计算程序的可靠性^[4].

6.2 极型场对流图像

图 2显示了深度为200km 处R-F 和F-F 边界下的极型场对流图像与地表构造的关系.与文献[17]中的图 3对比发现，本文结果相对于小横向黏度变化下的结果具有较好的改善.主要体现在太平洋板块、大洋洲和南美洲下的对流运动方向以及东太平洋洋中脊处对流状态.本文结果明显揭示出太平洋板块向西或西北、大洋洲向北或北北东、南美洲向西的运动状态以及对流发散区与东太平洋洋中脊的良好对应.而小横向黏度变化下，两种边界下的极型场对流图像难以明显揭示出现今地表板块的运动状态(图 3).对于欧亚大陆、非洲大陆以及北美洲下的对流运动状态，没有明显的改善.

图 2揭示出，对流运动在碰撞带(如喜马拉雅-阿尔卑斯碰撞带)和俯冲带(如西太平洋俯冲带)表现出汇聚流，而在洋中脊(如东太平洋、印度洋、大西洋洋中脊)则表现为发散流.对流速度较大的区域主要集中在西太平洋俯冲带和洋中脊附近.洋中脊如东太平洋和印度洋洋中脊两侧对流运动速率存在明显差异.这些结果与HS2-NUVEL1 模型获得的地表板块运动状态(图 3)比较一致.与小横向黏度下的结果一样，本文结果(图 1，特别是图 1a)也揭示出在非洲大陆的东北部下存在一个对流运动的发散中心，其北或北北东向对流速率明显较大，这或许是东非大裂谷形成的深部动力学原因.

300km 深度处R-F边界下极型场对流图像(图 4a)与200km 深度处的(图 2a)非常相似，也与地表大型构造带具有良好的对应关系.但400km 深度处R-F边界下极型场对流图像(图 4b)与200km(图 2a)和300km(图 4a)的特征存在明显的差异，主要体现在东太平洋、印度洋洋中脊处和非洲大陆下的对流发散区的消失以及印度洋板块下东北部对流运动速率的明显减小.这些现象可能表明与欧亚大陆作用的印度洋板块的主体为400km 以上的部分；产生东太平洋洋中脊和东非大裂谷的地幔热柱主要存在于400km 以浅.

从图 2和图 4看出，本文计算的极型场对流运动的最大速率已经明显大于现今板块运动的总速率(图 3)，这可能是在计算中引入的约束不够的缘故.

6.3 环型场对流图像

图 5显示了200km 深度处R-F 和F-F 边界条件下的环型场对流图像与地表构造的关系.与小横向黏度变化下的环型场对流图像(文献[17]中的图 5)相比，图 5主要改善之处为北太平洋下的环型场对流的旋转中心的位置更偏北以及南美洲和东太平洋边界处的旋转状态.通过与图 5 中显示的现今板块运动的径向涡度(radialvorticity)的等值线^[28]进行对比，发现这种改善使得本文获得的环型场对流比小横向黏度变化下的更加合理.

图 5揭示出总体上北半球处于右旋、南半球处于左旋循环对流状态，并且在赤道附近区域存在自东向西的环型场对流条带，在条带的两侧对流速率总体上存在差异：北半球对流速率相对较大且主要存在于环赤道低纬度区域，而南半球相对较小，这种差速流动或许是导致环赤道大型剪切系统^[29]及影响其强震活动的深部动力学原因.

对比图 2和图 5可以看出，计算的环型场对流速率远远小于极型场对流速率，这与现今板块运动的极型场与环型场能量几乎相等的事实还存在显著差异.同时使得地幔对流总速度场与极型场速度分布特征几乎相同(图 6).

7 结论与讨论

本文的结果相对于小横向黏度变化下的结果具有显著的改善.通过与地表板块运动进行对比分析后，发现本文的结果更加合理.这可能要归功于两个方面.首先，本文将横向黏度变化的量级显著提高，达到了3个量级.横向黏度变化大小的改变对计算结果进行了合理的调节；其次，由于横向黏度变化量级的提高，使得在求解过程须经过交叉迭代，这就使得极型场和环型场达到了真正的耦合，从而结果更趋合理.不过，计算获得的极型场对流最大速率明显大于、环型场对流最大速率明显小于现今地表板块运动的最大速率；某些区域如欧亚大陆和非洲大陆下的对流状态与现今地表板块的运动状态有较大差异.

自20世纪80年代以来，很多研究者已经对环型场的产生机制及其地球物理意义做过深入的研究.主要从两个方面入手.首先，利用运动板块激发地幔的环型对流.在这类模型中，一般要求板块内部为刚性，而边界则为变形区域(黏度比内部低得多)，并且作用于板块的总应力或力矩为零.研究结果表明，能获得能量几乎相等的极、环型场和良好的地表速度场以及大地水准面异常的拟合，但是环型场能量随着深度会快速衰减，主要存在于上地幔中，而极型场却在整个层中相对稳定(除两个边界处有所衰减外)^[30~33].然而，由于板块边界处应力奇点的存在导致施加于旋转板块上的力矩趋于无穷会导致模型不自恰；其次，利用地幔黏度横向变化产生环型对流.地幔黏度横向变化的引入，即使在均匀边界下，也能够自然地获得环型场分量.计算的环型场对流速率远远小于极型场对流速率^{[7, 16]}，本文的结果与此一致，但与板块运动的极型场分量和环型场分量功率几乎相等具有明显差异.由此，建立同时包含了板块自身运动和地幔黏度横向变化的对流模型可能对解释地表观测能提供更大的帮助.

附录A

(A1)

(A2)

(A3)

(A4)

(A5)

(A6)

(A7)

(A8)

(A9)

(A10)

(A11)

(A12)

(A13)

(A14)

(A15)

(A16)

(A17)

(A18)

(A19)

(A20)

(A21)

(A22)

(A23)

(A24)

(A25)

(A26)

(A27)

(A28)

(A29)

(A30)

(A31)

(A32)

(A33)

(A34)

(A35)

(A36)

(A37)

(A38)

(A39)

(A40)

(A41)

(A42)

(A43)

(A44)

(A45)

(A46)

(A47)

(A48)

(A49)

(A50)

(A51)

(A52)

(A53)

(A54)

(A55)

(A56)

(A57)

(A58)

(A59)

(A60)

(A61)

(A62)

(A63)

(A64)

(A65)

(A66)

(A67)

(A68)

(A69)

(A70)

(A71)

(A72)

(A73)

(A74)

(A75)

(A76)

(A77)

(A78)

(A79)

(A80)

(A81)

(A82)

(A83)

(A84)

(A85)

(A86)

(A87)

(A88)

(A89)

(A90)

(A91)

(A92)

(A93)

(A94)

(A95)

(A96)

(A97)

附录B

(B1)

(B2)

(B3)

(B4)

(B5)