上地幔变粘度小尺度对流的数值研究
叶正仁王建
（中国科学地质与地球物理研究所, 北京 100101）
[摘要]本文基于一个二维模型，利用有限元方法，研究上地幔—岩石圈系统的变粘度小尺度对流。
考虑该系统的粘度随温度以指数形式变化。数值结果表明，当粘度随温度变化较剧烈时，由于低温高粘，系统的最上部物质不参与对流，系统发育形成一个类似于岩石圈的静止盖层。计算的表面热流、地形起伏及重力异常与对流格局有较好的相关性：高热流、上升地形、高重力异常对应于对流的上升区；反之低热流、下降地形、低重力异常与对流的下降区对应。
[关键词]小尺度对流，有限元方法，静止盖层，热流。
1 引言
研究表明，在地球内部除了与板块运动相联系的全球尺度的大规模对流外，还很可能存在小尺度的上地幔对流[1-3]。Haxby与Weisl [4]运用小尺度对流理论成功解释了观测到的太平洋和印度洋重力和大地水准面区域异常分布。近年来，随着对大陆地区岩石圈结构和动力学的研究，若干作者探讨了大陆地区上地幔小尺度对流存在的可能性及其产生的地球物理效应.[5-8]。philex
us news
高压下矿物的蠕变试验表明[9]，地幔、特别是岩石圈的粘度强烈地依赖于温度、压力。关于变粘度问题的对流研究, 自上世纪八十年代以来已有相当多的工作。由于问题的复杂性, 几乎所有的工作都是利用数值方法, 包括有限元及有限差分。Christenn等[10-12]就二维及三维的变粘度对流进行了系统的研究, 数值模拟结果表明, 粘度结构对对流的形态、格局及内部热状态有非常重要的影响。Tackley[13]及 Zhong等[14]则分别在直角坐标和球坐标下，探讨了粘度随温度变化的本构关系对全球地幔对流的影响。本文就一个二维模型，利用有限元方法，研究上地幔变粘度小尺度对流。与前人的工作比较, 其主要的区别在于所研究的区域有比较大的长宽比(本文为5), 而周知，对上地幔小尺度对流而言, 其可能的水平延伸要比深度扩展大得多。除了考虑对速度场的影响外、我们还探讨粘度随温度变化对内部温度分布、表面热流、地形变化及重力异常等地球物理场的影响。
2 物理模型和数值方法
所研究的区域为上地幔—岩石圈系统。上边界为地表；下边界为上、下地幔分界面。模型的长宽比为5。假设Boussinesq近似成立，并且内部没有热源，则对应的连续性方程、动量方程和能量方程分别为：
_________________________________________________________________________
[基金项目] 国家自然科学基金（49974022）及中国科学院知识创新工程项目 (kzcx2-112)。west是什么意思
[作者简介] 叶正仁，男，1944年生，1967年毕业于北京大学地球物理系，1981年硕士研究生毕业，现为中国科学院地质与地球物理研究所研究员，博士生导师。主要研究领域为地幔动力学、重力学及深源地震成因等。
0=⋅∇V (1a) 0=+⋅∇+∇−g σδρp (1b)
)(T k T t
T ∇⋅∇=∇⋅+∂∂V (1c)
式中为速度，p 为动力学压力，σ是偏应力张量，V T 为温度，)(1ρρδρ−=为密度异常，1ρ为参考密度。 是热扩散系数，假定为常数。在写下（1b ）式时，已考虑到系统具有非常大的Prandtle 数P k r （Pr =k
ν , 式中ν为运动粘滞系数）
，从而动量方程中的惯性项可忽略。 偏应力—速率的关系为：
)(∇+∇=V V σµ (2)
在只考虑热浮力的情况下，密度异常δρ与温度有如下关系[9]：
δρ=)(11T T −−αρ (3) 式中α是热膨胀系数，µ为粘度，T 为参考面（本文取为上表面）上温度。
1尽管在不同的压力、温度及应变率条件下，地幔—岩石圈系统有着不同的变形机制，但其中对粘度具有最主要的影响的因素是温度。本文假定粘度与温度具有如下的指数关系，即Frank-Kamenetskii 表达式[15]：
)](exp[)(11T T T −−=γµµ (4) 式中，1µ 表示上边界的粘度，γ是一个常数，它的大小控制着系统内部的粘度差异。
为便于分析及数值计算，使用如下参数对原始方程组进行无量纲化（带“′”的量为无量纲量）：
; ; ; ; ; ; ),(),(20
012p d
k p T T T T t k d t d k d y x y x ∗=′∗=+′∗∆=′∗=∗=∗′′=µµµµV V (5) 式中，为层厚，d T ∆为上地幔底部与表面的温度差。0µ为上地幔底部粘度。
得到无量纲方程组为：
0=′⋅∇′V k T Ra p ′⋅−=′⋅∇′+′∇′−σ
T T t T ′∇′=′∇′⋅′+′∂′∂′2V (6) )(∇′′+′∇′′=′V V µσ
)exp()exp()(0101T T T T ′−=′∆−=
′′θµµγµµµ
k 为方向（向上）的单位矢量。系统的运动形态由二个无量纲数数及y Ra θ值表述，这是与等粘度对流不同的，在后者情况下，系统运动形态仅由数确定。
Ra k
Td Ra 03
pgmµαρ∆=g (7) T ∆=γθ (8) 而在上下表面的粘度差异µ∆，由式)exp(θµ=∆确定。
福州成人英语学习边条件：为比较边条件的影响，取二组不同的边界条件。
（1） 所有表面为应力自由面；上表面温度01=T ，下表面温度为1700K(无量纲温度
)，侧面为绝热；
1'0=T （2） 其余同（1），除了下表面以热流输入边条件代替等温条件，即在下表面有
20的热流输入2/m mW [2]；
初始温度分布为传导温度解：T y t ′−==′1)0(。
我们利用有限元方法(FEPG 有限元程序包)解方程组(6)。具体做法是：给定初始温度场，计算粘度分布并解动量方程得到速度场，而后利用此速度场代入能量方程解下一个时间步的温度场，再利用此时间步得到的温度场解速度场，依次循环。当4)()1(10−+≤n n T T
时可认为温度达到稳定。这里)1(+n T 表示第（n+1）步、)(n T 表示第（n ）步的平均温度场。对于本文所取的参数范围，数值计算表明温度场可达到稳态。需要指出的是，在能量方程中由于对流传热项V ⋅的存在，用通常的Galerkin 有限元法将产生虚假的数值振荡，为此采用Petrov-Galerkin 方法T ∇[16]，此法能很好地抑制解的振荡而得到平滑解。计算中网格剖分数为；对速度及温度使用9节点四边形单元；而对压力4832×P ，采用4节点。双线性插值。能量方程中的时间步长遵从Courant 准则[16]。
3结果
表1列出了数值计算所用的上地幔—岩石圈系统的参数值。
表1
参数 物理含义 值
d 层厚
km 670k 热扩散系数 10
s m /26−α 热膨胀系数
5102−×T ∆ 上下界面温度差 1700
K q 进入下表面的热流
jhs2/20m mW g 重力加速度
conclusions2/8.9s m Ra Rayleigh 数 3×105
µ∆ 上下边界的粘度差异
5100−
图1为下边界为等温边条件时的上地幔—岩石圈系统的流场图。从上而下(a-f)上下表面
hm是什么东西的粘度差异依次是10。
它们共同的特征是都有良好发育的纵横比近似为1的对流环（图1a （等粘情况）即经典的Rayleigh-Banard 对流），但随着内部粘度差异的加大，上部物质的流动速度逐渐变小，直至当粘度差异等于大于10时，上部物质基本上不参与对流，系统产生一个静止的盖层（Stagnant lid ）。一些作者5
4321010,10,10,10,10,3[2,15]曾将其比拟为岩石圈。相应地，温度场随着粘度差异的变化示于图2中。该图给出对不同粘度差异（10）时的温度等值线。显而易见，当粘度差异比较小时（图2a-c ），上下边界的热边界层均发育良好，而随着粘度差异的增大（图2d-f ），上边的热边界层厚度逐渐加大。另一个值得注意的情况就是随着粘度差异的加大，“热柱”与“冷柱”的宽度增大。
5010
no other
−
图1 下边界为等温边条件时的流场图（无量纲）。 从上而下(a-f)，上、下表面
的粘度差异依次是100,101,102,103,104,105。
Fig.1 The calculated non-dimensional velocity field with isothermal lower boundary
condition. From top to bottom (a-f) the viscosity contrasts between the
upper and lower boundaries quently are 100,101,102,103,104,105.
glimp
.
图2 下边界为等温边条件时的温度分布（无量纲）。粘度差异同图1。
Fig.2 The calculated non-dimensional temperature distribution with isothermal lower boundary condition. The viscosity contrasts are the same as that in Figure 1 .
图3和图4 分别为下边界为热流边条件的情况下，对应不同粘度差异所得到的流场和温度分布。与图1及图2对比，可以发现在热流边条件情况下，同样随着内部粘度差异的增大，也存在如下的情况：一个静止的盖层、热边界层的增厚、热柱扩张。当然，细节有所不同。这表明，相对于边条件而言，系统内部的物性（主要是粘度随温度的变化）对对流形态和格局，温度分布起着控制作用。
