岩石圈弯曲的黏弹性松弛及其对确定弹性岩石圈厚度的影响*

邵翠法，王世民

(中国科学院大学地球与行星科学学院 中国科学院计算地球动力学重点实验室， 北京 100049)

地壳和地幔中的固体岩石具有显著的黏弹性性质，即在低温和快速变化的载荷作用下变形由可逆的弹性应变占主导，而在高温和缓慢变化的载荷作用下变形以不可逆的黏性流动为主。

弹性岩石圈的定义以及弹性岩石圈厚度的确定是地球动力学研究的中心课题。事实上，弹性岩石圈并没有一个截然的下界面。弹性岩石圈的一种定义方法是借助于岩石的黏弹性松弛时间(黏度与弹性模量之比)。如定义弹性岩石圈内松弛时间大于一千万年[1]，相当于干橄榄石温度低于675 ℃、湿橄榄石温度低于475 ℃。考虑到一般以1 300 ℃等温线定义热岩石圈与软流圈界面[1]，弹性岩石圈应看作是热岩石圈的一个组成部分。为论述方便起见，下文中一律将弹性岩石圈以下的热岩石圈部分称作黏性岩石圈。显然，整个热岩石圈都具有黏弹性性质(虽然在弹性岩石圈中以弹性变形为主，而在黏性岩石圈中以黏性变形为主)，且弹性岩石圈向黏性岩石圈的过渡是主要由温度连续变化导致的渐变而非突变。根据松弛时间定义的弹性岩石圈厚度在原理上无懈可击，但由于地下岩石的黏度无法精确测量，并不能实际应用于具体区域的弹性岩石圈厚度确定。

在地球动力学研究中，弹性岩石圈厚度通常由弹性板弯曲分析反推，即通过以弹性薄板在静岩恢复力(浮力)作用下弯曲的解析解[1]拟合实测岩石圈弯曲导致的地形(如岛链、海沟的外缘隆起位置)来反演岩石圈的抗弯刚度进而确定弹性岩石圈等效厚度。这样得到的弹性岩石圈厚度与根据松弛时间定义的弹性岩石圈厚度是否一致并没有经过严格的检验。

严格来说，弹性岩石圈的厚度应当是载荷作用时间的函数。在地震波传播的时间尺度上，整个地壳和地幔都表现出弹性性质，而在地质时间尺度上，弹性变形只在近地表的弹性岩石圈内占主导，黏性流动则主导弹性岩石圈以下地幔岩石的变形(如地幔对流)。岩石圈对各种不同时间尺度载荷的响应可以反映相应时间尺度下的弹性岩石圈厚度。如对震后变形[2-7]、冰后回弹[8-10]、俯冲过程[11-14]和岩石圈弯曲[15-19]等问题的黏弹性模拟可以约束相应载荷作用时间下的弹性岩石圈厚度。

描述岩石圈弯曲的物理模型经历了从简单到复杂的过程。单层弹性薄板加无黏流体地基的模型[20-28]是描述岩石圈弯曲变形最简单且应用最广泛的模型。通过对比模型预测和火山、沉积物等载荷作用下地表形态的实际观测数据，可以研究弹性岩石圈厚度与板块年龄及载荷作用时间的依赖关系[29-31]。基于这种模型，Watts等[32]曾通过海底地形反演计算全球范围内9 758处弹性岩石圈的厚度。

为克服弹性板弯曲模型结果与时间无关的缺陷，之后又发展出了均匀黏度黏弹性薄板模型[33-36]。 Walcott[33]采用均匀黏度黏弹性薄板模型研究大洋和大陆岩石圈的弹性厚度; Lambeck 和 Nakiboglu[37]给出多层弹性板、均匀黏度黏弹性板以及弹性层叠加黏弹性层模型的解析解，并计算海洋岩石圈的有效弯曲刚度。为了更好地反映黏度随深度的变化，文献[18, 38]采用多层黏弹性模型，并解释弹性厚度大范围变化的观测数据。

除解析方法外，数值模拟方法也越来越多地应用于岩石圈弯曲问题研究。如，Zhong等[10]采用三维球坐标有限元模型研究不可压缩地球介质在自重作用下的冰后回弹问题；Zhong 和 Watts[17]采用三维直角坐标有限元模型计算岩石圈在表面载荷作用下的黏弹性变形和应力，其中考虑了复杂的非线性黏弹性流变效应，但同样假设材料为不可压缩的，与地震波速给出的岩石圈泊松比远低于0.5的观测事实明显不符。

本文采用黏度随深度连续变化的黏弹性本构模型统一描述弹性岩石圈与黏性岩石圈，并建立有限元模型求解在表面垂向载荷(对海山或岛链载荷等的一般化处理)作用下的岩石圈弯曲问题。利用经典弹性板弯曲解析解，我们将本研究黏弹性模拟计算得到的外缘隆起位置随载荷作用时间的变化转换成弹性岩石圈厚度随时间的变化，并与由松弛时间定义的弹性岩石圈厚度随时间的变化进行定量对比。因此，本研究不但实现了弹性岩石圈向黏性岩石圈的逐渐自然过渡，而且能够对弹性岩石圈两种传统定义(即松弛时间定义和弹性板弯曲定义)的一致性加以严格检验。

1 弹性岩石圈的两种传统定义及其联系

弹性岩石圈的传统定义包括弹性板弯曲定义和松弛时间定义两种[1]。

在弹性板弯曲定义中，将弹性岩石圈之下的地幔介质近似为应力完全松弛因而处于静水压状态的流体，而弹性岩石圈由弹性薄板近似，如图1所示[1]。

图1 大洋岩石圈受外载荷作用而发生弯曲Fig.1 The oceanic lithosphere deflected by an applied load

依据弹性薄板弯曲理论，利用实测岩石圈在表面载荷作用下发生弯曲导致的地形(如岛链、海沟的外缘隆起位置)能够反演岩石圈的抗弯刚度进而确定弹性岩石圈等效厚度[39]。有效弹性厚度Te与薄板抗弯刚度D之间的定量关系为

(1)

式中：弹性参数E和ν分别为杨氏模量和泊松比。

根据松弛时间定义[1]的弹性岩石圈实质上是要求弹性岩石圈内的黏弹性松弛时间大于载荷作用时间，即载荷作用在弹性岩石圈中尚未松弛，而在下伏黏性岩石圈中已经松弛。根据Turcotte和Schubert[1]的分析，初始应力σ0松弛一半需要的时间为

(2)

假定杨氏模量E=70 GPa，松弛时间τr=107a，初始应力σ0=100 MPa，对于干的橄榄石Ea=523 kJ·mol-1，C1=4.2×105MPa-3·s-1，由公式(2)给出弹性岩石圈底界相当于干橄榄石温度675 ℃，而对于湿的橄榄石Ea=398 kJ·mol-1，C1=5.5×105MPa-3·s-1，由公式(2)给出弹性岩石圈底界相当于475 ℃等温线。

以上两种弹性岩石圈的定义之间没有直接联系。基于松弛时间定义的弹性岩石圈理论上合理，但因黏度难以直接测量无法广泛应用，而基于弹性板弯曲定义的弹性岩石圈则广泛应用于实际弹性岩石圈厚度的确定。

2 有限元模型

2.1 物理模型

本研究采用的物理模型如图2所示。在线状分布的地表垂向载荷(海山或岛链载等)作用下的岩石圈弯曲问题近似模拟为平面应变问题。岩石圈底部采用弹簧型边界条件模拟静水恢复力。上表面加载区域指定为垂向载荷线性变化，其余区域指定为自由边界。左侧边界指定为对称边界，而右侧边界位移固定。

图2 物理模型Fig.2 Physical model

我们用Maxwell模型黏弹性体代表岩石圈。假设弹性变形体积可压缩，但弹性模量和泊松比不随深度变化；黏度则指定为随深度连续变化的情形[1]

(3)

式中：C为指前系数、Ea为活化能、T为绝对温度、R为普适气体常数。因为岩石圈所受压力较低，忽略压力对黏度的影响。本文考虑的黏度变化是与Karato 和 Wu[40]的结论相符的。他们指出位错蠕变在高温(例如大洋中脊)时是重要的，而低温时扩散蠕变占主导。由于海山或岛链载荷作用的时间尺度较大(>1 Ma)，在这个时间尺度内岩石圈有相当部分的浅部低温区域经历了应力松弛，因此扩散蠕变是一个合理的假设[18]。设软流圈温度为Tm=1 300 ℃，指前系数C可以用参考软流圈黏度ηref表示

(4)

因而式(3)可以重新写为

(5)

只要给定温度随深度的变化以及活化能Ea与参考软流圈黏度ηref就可以用式(5)确定黏度随深度的变化。Watts和Zhong[18]以及Watts等[16]曾用这样的扩散蠕变黏度公式结合板块冷却模型[41]和分层均匀黏弹性结构研究岩石圈弯曲和等效弹性厚度问题。他们推荐的参数组合是Ea=120 kJ/mol，ηref=1020Pa·s。本文采用同样的参数设定，以方便模型及结果比较。

从式(5)可以看出温度是控制黏度的关键因素。我们采用板块冷却模型[41]计算岩石圈内温度分布

T(y)=T0+(T1-T0)

(6)

式中：y是深度，ts为大洋岩石圈年龄，yL0是热岩石圈最大厚度，κ是热扩散系数，Tm和T0分别是软流圈和海底温度。模型参数取值见表1。

表1 模型参数Table 1 Modeling parameters

对年龄给定为80 Ma的大洋岩石圈，得到的温度和黏度随深度变化分别由图3(a)和3(b)给出。

图3 模拟采用的岩石圈温度和黏度剖面Fig.3 Adopted profiles of temperature and viscosity in the modeled lithosphere

假设岩石圈上表面受到不随时间变化的垂向载荷作用是关于y轴对称的，为减小计算量，只考虑x>0部分的岩石圈弯曲即可。垂向载荷函数指定为

(7)

由于模拟问题是线性的，计算位移与应力均与施加的载荷成正比，而黏弹性松弛过程及外缘隆起的时间演化均与施加载荷的绝对大小无关。

岩石圈下表面受到软流圈静水恢复力与岩石圈下表面垂向位移成正比，即

q(x)=kv(x),

(8)

这里k=(ρm-ρw)g，ρm和ρw分别是岩石圈之下软流圈和岩石圈之上海水的密度，g是重力加速度，v(x) 是岩石圈下表面的垂向位移。

2.2 数学模型

岩石圈弯曲过程中，应力张量σij满足准静态平衡方程

σij,j=0，

(9)

位移与小应变之间满足几何方程

(10)

式中：ui为i方向上的位移，而应力与应变由体积可压缩的三维Maxwell体本构方程相联系，即

(11)

式中：λ,μ为弹性拉梅常数；η为黏度(见式(5))。利用平面应变εzz=0的条件，三维问题可简化为二维问题。

2.3 数值模型

用有限元软件COMSOL[42]数值计算上文定义的岩石圈弯曲问题。模拟的黏弹性岩石圈为2 000 km×81 km的长方形区域。除精度检验外，本文结果一律基于参考模型计算。在参考模型中，采用5阶三角形单元对计算域进行剖分，共有126 520个单元，总自由度数为8 864 562。在区域的右半部分最大单元尺寸小于2 km，而x

图4 部分模型域的计算网格剖分Fig.4 Computational mesh in part of the modeling domain

2.4 计算精度检验

为了验证数值模拟的计算精度，我们比较了不同网格尺寸、单元阶次和时间步长的模拟结果。

在网格单元尺寸的检验中，在参考模型网格的基础上，将最大网格尺寸加倍，单元数从参考模型的126 520减少到31 906，自由度从8 864 562减少到2 237 512。在单元离散阶次的检验中，将单元阶次从参考模型的5阶减小到4阶，自由度从8 864 562减少到5 826 450。在计算时间步长的检验中，将参考模型每一阶段采用的时间步长都增加到2.5倍，4阶段的最大步长依次为250、2 500、25 000和250 000 a。

以上3种检验模型的结果与参考模型比较由图5示出。图中对比的是上表面t=100 Ma时刻x=0处的垂向位移。网格变疏、离散阶次减小、时间步长增大的结果在图4中依次由加号、叉号、星号标记对应的线表示。结果显示4条曲线几乎完全重合，网格变疏、离散阶次减小、时间步长增大的3种检验模型相对于参考模型相对误差最大值的量级分别为10-6、10-8和10-9，证明本文得到的结果对计算网格、单元阶次及时间步长的选取不敏感，具有足够的计算精度。

图5 数值模拟计算精度检验Fig.5 Accuracy tests for numerical modeling results

3 模拟结果与讨论

3.1 模拟结果

本节给出的模拟结果均基于参考模型。岩石圈上表面的垂向位移随着载荷加载时间的变化如图6所示。图6(a)与6(b)分别绘出t=0+～9 Ma与t=9～100 Ma两个时间段的垂向位移随x坐标的变化。图6(a)中的t=0+解是程序给出的瞬时弹性位移解(图8中t=0+结果表示瞬时弹性应力)。结果表明，前一时间段外缘隆起向内(加载一侧)运动，而后一时间段外缘隆起向外移动，即随着载荷作用时间的增加，外缘隆起与载荷中心处的距离先减小后增大。

图6 岩石圈上表面的垂直方向位移随载荷作用时间的变化Fig.6 Vertical displacement of the upper surface of the lithosphere as a function of loading time

上表面载荷中心处(x=0)的垂向位移，以及上表面x=250 km处的垂向位移与水平位移随时间的连续变化如图7所示。图7除清楚表明载荷中心与外缘隆起相反的垂向运动外，也清楚地反映了外缘隆起区域存在先内后外的水平运动特征。

图7 计算位移的时间演化Fig.7 Evolution of the computed displacement

数值模拟得到的黏弹性应力松弛过程由图8给出。图中绘出的是x=0处水平和垂向正应力σxx,σyy以及x=250 km处的水平方向正应力σxx和剪应力τxy在不同时刻随深度的变化。图8表明，岩石圈中应力的黏弹性松弛是一个持续的动力学过程。在加载之后，岩石圈整个厚度均产生弹性响应，水平正应力随深度近乎线性分布，与弹性薄板假设很接近。随着载荷作用的持续，岩石圈深部的应力逐渐松弛，而浅部的正应力持续增加，岩石圈的弹性厚度因而持续减小。外缘隆起下方区域的剪应力呈现先增加后减小的趋势。这样的剪应力演化过程与外缘隆起区域先内后外的水平运动特征是一致的。

图8 岩石圈黏弹性应力松弛过程Fig.8 Processes of viscoelastic stress relaxation in the lithosphere

3.2 讨论

Turcotte和Schubert[1]给出利用岩石圈弹性弯曲产生的外缘隆起位置反推出弹性岩石圈厚度的经典公式

(12)

由外缘隆起的位置xb可计算弹性岩石圈厚度Te；Te随xb的单调变化图像如图9(a)所示。图9(b)给出本研究计算得到的外缘隆起位置随时间的演化规律。图9(b)表明，外缘隆起与载荷中心的距离先是持续减小，直到t=9 Ma达到最小值242 km。而后，外缘隆起与载荷中心的距离又随时间持续单调递增直到1 000 Ma。

结合图9(a)和9(b)，即将本文模拟得到的外缘隆起位置变化由公式(12)转换成弹性岩石圈厚度，得到图10所示的Te随加载时间的演化曲线。由于外缘隆起区域先内后外的水平运动导致图10中Te先随加载时间减小而后又随时间增加。显然，弹性厚度随时间增加阶段的结果是违反物理规律的。因此，可以得到结论：前人广泛采用的利用外缘隆起位置反演弹性岩石圈厚度的做法具有严重缺陷，在载荷作用时间超过9 Ma以后会严重高估弹性岩石圈厚度。

事实上，即使在载荷作用时间小于9 Ma的时间内，利用外缘隆起位置反演得到的弹性岩石圈厚度与根据松弛时间定义的弹性厚度也有明显的区别。根据Turcotte和Schubert[1]的公式

(13)

我们计算了本文所用黏度剖面对应的深度yEL随松弛时间的变化关系，并绘制在图10中。如果用松弛时间大于载荷作用时间作为弹性岩石圈的定义，显然松弛时间曲线对应于弹性厚度的另一种定义。

图10清楚地表明，两种定义给出的弹性厚度仅在0.2～5.0 Ma的时间范围内较为符合(差异不超过0.64 km， 且在t=0.30 Ma和t=4.5 Ma时两线重合)。在载荷作用时间小于0.2 Ma时，两种定义有显著差别但总体趋势一致。在t=0.001 Ma时两条曲线差距最大(4.4 km)，接着差异慢慢减小到零(t=0.001 8 Ma)，随后逐渐增大到3.64 km(t=0.008 Ma)后又逐渐减小。在t>5 Ma时，两种定义给出的弹性岩石圈厚度差异随时间快速增大。由于松弛时间定义的弹性厚度yEL是直接由模型给定的光滑单调变化的黏度剖面计算得到的，图10中两条曲线间差距的复杂交替变化完全由岩石圈黏弹性松弛导致的外缘隆起移动特征造成。由图9(b)可以看出，外缘隆起位置随加载时间对数的变化先后经历了斜率增加到接近常斜率、而后斜率又逐渐减小为零、最后斜率反向并逐渐增加的复杂演变。正是岩石圈黏弹性松弛的复杂性导致弹性岩石圈两种定义给出的厚度曲线出现多次交叉的现象。

图10 弹性岩石圈两种传统定义给出的厚度对比Fig.10 Comparison between the thickness values calculated based on two classical definitions of elastic lithosphere

4 结论

本文利用变黏度可压缩黏弹性有限元模型定量研究岩石圈弯曲的黏弹性松弛问题。模拟结果显示，外缘隆起位置具有先向弯曲中心运动后向弯曲外侧移动的特点，反向运动发生在载荷作用约9 Ma后。由于经典弹性板弯曲理论给出外缘隆起与弹性岩石圈厚度的单调正相关关系，据此反演的弹性岩石圈厚度不能应用于载荷作用超过9 Ma情形。否则，会出现弹性岩石圈厚度随加载时间增加的荒谬结果，违反岩石圈黏弹性松弛的物理规律。进一步与基于松弛时间的弹性岩石圈厚度定义进行对比，本文模拟结果表明，由基于弹性板弯曲解利用外缘隆起位置约束的弹性岩石圈厚度仅在载荷持续0.2～5.0 Ma的时间段内能够给出较为合理的结果，误差不超过0.64 km。在载荷作用时间小于0.2 Ma时，基于两种传统定义给出的弹性岩石圈厚度总体趋势一致，但差异最大可达4.4 km。而在载荷作用时间大于5 Ma时，两种定义给出的弹性厚度差异随时间快速增大。两种传统定义给出的弹性岩石圈厚度曲线出现多次交叉的现象，明确揭示岩石圈弯曲的黏弹性松弛是一个复杂的动力学过程，因此在确定弹性岩石圈厚度的地球动力学相关研究中应当充分考虑黏弹性松弛的影响。