地球物理学报  2013, Vol. 56 Issue (8): 2625-2635   PDF    
青藏高原通道流模型动力环境的数值模拟
杨辉 , 滕吉文 , 皮娇龙     
中国科学院地质与地球物理研究所, 北京 100029
摘要: 中、下地壳流模型作为一种可能的动力学演化机制, 在解决诸如喜马拉雅造山带和青藏高原东缘、南缘等区域地壳中岩层的通道流或韧性剪切挤出等方面的解释给出了相应的模型和阐述.本文基于青藏高原壳、幔介质平均速度模型, 采用二维黏弹性数值模型对高原下地壳物质流动的动力学边界条件进行探讨.研究结果表明:(1)青藏高原下地壳与上地幔盖层物质作为坚硬的固态物质相接, 不具备可运动的边界条件, 难以在Moho界面处任意地域发生相互运动.壳、幔介质中需存在可供物质高速运动的边界条件, 即以上地壳底部的低速层为上滑移面, 以上地幔软流圈顶部为下滑移面, 才有可能在足够强的力系作用下促使"下地壳+岩石圈盖层"物质发生同步运移; (2)若不具备这样的初始与边界条件是难以产生深部物质运移的.因此, 青藏高原深部壳、幔物质运动不可能是普遍存在的, 只能是局部和在特异环境下才能实现.
关键词: 下地壳流      青藏高原      初始和边界条件      数值模拟      滑移面     
Numerical simulation of the geodynamical condition about the channel flow model in Tibetan plateau
YANG Hui, TENG Ji-Wen, PI Jiao-Long     
Institute of Geology and Geophysics, Chinese Academy of Sciences, Beijing 100029, China
Abstract: Being a possible geodynamical mechanism, "crustal flow" comes up with corresponding models in resolving and predicting the tectonics, metamorphism and exhumation of high terranes in the Himalayan-Tibetan orogenic system. This paper adopts 2-D viscoelastic models to discuss the geodynamical boundary condition for the deep material flow in Tibetan plateau, based on the average geophysical model of the crust and upper mantle. The result shows: (1) As non-smoothly contact mode, the solid materials of the lower crust and upper-mantle lid beneath the Tibetan plateau do not have the boundary condition to move mutually, which means that "channel flow" model cannot independently occur freely nearby Moho interface. The simultaneous transfer of the deep medium of "lower crust and upper-mantle lid" can possibly happen by the enough force source, which takes the low-velocity layers as upper detachment interface at depths 20~30 km, and the top of asthenosphere as lower detachment interface at depths 100~110 km. (2) "Channel flow" is hard to move without the former sufficient and necessary conditions. So, the deep materials flow in Tibetan plateau is not universal, and only happens in local and particular circumstance..
Key words: Lower crustal flow      Tibetan plateau      Sufficient and necessary conditions      Numerical simulation      Detachment interface     
1 引言

岩石圈具有一定的流变性质, 这是控制岩石圈构造变形的主导因素之一, 也是探讨岩石圈动力学的基础.因此, 从20世纪70年代开始, 许多学者便开始对岩石圈流变性质研究, 对大陆构造运动的动力条件及大陆在其作用下的运动和形变场进行研究, 以探讨地球深部各圈层及地表各地质单元块体的规律性运动和形变过程.但利用矿物和岩石的实验规律来定量或是半定量地外延去定义岩石圈分层结构的流变参数, 从准确性及精度上尚远远不能满足大陆动力学及板块运动学的要求, 这也使得在对不同区域基于地质、地球物理和地球化学等观测数据而得出的构造演化的动力学过程解释中出现不同的, 甚至相互矛盾的解释[1].

当前, 对岩石圈内部的流变机制认识还相当有限, 主要是基于地震波场和地热场方面的资料, 而这些地球物理量的真实测定及解释还存在很大的不确定因素, 且地球介质的属性并不全等同于物理上的流变性, 因而准确地说应当为"视流变性".显然在这种尚不能给出地球深部介质真实黏滞系数和全基于物理上的流变机制状况下, 采用各种假设条件去试图解释岩石圈动力学机制, 这会给理解岩石圈流变结构增加很大的困难, 乃至在解释方面引发很大争议.为此, 须在这一基点上对地球深部介质属性进行探索.

以青藏高原为例, 对这个由印度板块与欧亚板块因碰撞而形成的地球上最高最大的高原, 国内外已有不少人涉足于这一方面的研究, 但是对其从表层到深部的三维变形机制, 至今仍众说纷纭.近年来, 争论的焦点主要围绕着两个端元模型:一是岩石圈连续变形模型[2-3], 认为高原在演化发展的过程中, 变形从地壳到地幔是连贯且耦合的; 一是通道流(channel llow)模型[4-7], 认为变形是随深度变化的, 是由下方的软流圈流体驱动或是由上地壳拖曳运动.

可以说, 如何对研究区动力学模型进行简化, 特别是对边界条件的正确理解, 会给动力学响应过程的合理解释带来很多方面的影响.我们已清晰地认识到, 由于实验手段和方法等方面的限制, 地球科学中对物理学中流变学的借用, 与物理学中的流变学属性相比并非同一个量级, 而同时又因为当前尚缺乏一个对地球深浅部结构介质属性黏度的真实物理量的直接测定, 这不仅限制了对地球动力学和大陆动力学认识的深化, 更带来了对地球深部介质流变性量化的判断方面的困难.

因此, 本文基于青藏高原已经开展的深部地球物理探测工作所得到的一些规律性结果, 针对当前地球科学界流行的通道流模型的一些理论或假设基础, 对青藏高原下方物质流动所需要的边界条件开展了数值模拟工作, 以探讨该现象存在的一些必要条件.

2 问题的提出

自Kuhn[8]于1979年提出了通道流的原始概念以来, 三十年间, 通道流的模型得到了迅速的发展.通道流模型认为以青藏高原为例的一些造山带区域, 相对于较为刚性的上地壳和岩石圈上地幔, 处于中间层的中、下地壳作为强度较低的塑性流变层, 具有更强的流动性, 因而在高原厚地壳的重力差异驱动或是上地幔的拖曳作用下易发生流动而形成通道流[9].最初, 通道流的模型运用于对造山带腹地存在高阶变质岩的地带进行构造解释, Gruic[10]首先将通道流模型(Poiseuille和Couette流)力学机制用于喜马拉雅造山带腹地的高喜马拉雅结晶岩系(Greater HimalayanSequence, GHS)的韧性挤出构造演化模式中.高喜马拉雅结晶岩系位于北部的藏南拆沉系(South Tibetan Detachmentsystem, STD)及南部的主中央断裂带(Main Central Thrust, MCT)之间, 出露大范围的前震旦纪变质岩, 一般由高级变质岩组成, 属于印度板块北缘结晶基底的组分, 剪切变形的应变场特征剧烈.藏南拆沉系为北倾的低角度正断层系, 而主中央断裂带为北倾的逆冲断裂带, 两条断裂系近于平行, 具有相反的剪切特征, 活动特征具有明显的同期性, 约发生于22 Ma左右, 在中中新世中止[11].通道流模型作为其中一种可能的演化方式, 认为位于藏南拆沉系和主中央断裂系之间北倾的高喜马拉雅结晶岩系是中地壳物质向南挤出的条带区域, 可以说在解释高喜马拉雅结晶岩系的构造演化和应变模式方面得到了很好的应用.

随后, NelSOn[12]在对藏南构造演化研究时首先明确提出了"channel low"一词, 通道流模型得到了进一步的发展.期间, 青藏高原众多深部地球物理探测工作, 如重力与地表地形、地热、大地电磁和地震测深等资料解释对通道流模型的发展提供了帮助, 同时, 还开展了相应的数值模拟工作[7, 13-17].这些工作的展开将通道流模型运用到了更多的构造背景[18], 如:软流圈倒流(asthenospheric counterllow); 下地壳通道流(lower crustal channels); 造山带壳内通道流(intra-crustal channels); 俯冲带通道流(subduction channels); 盐构造通道流(salt tectonics)等, 并且, 其应用区域已不仅仅局限于喜马拉雅造山带和青藏高原周缘, 还应用到如希腊的西奈山(Hellenides)[19]的中地壳韧性变形挤出模式和阿巴拉契亚内山麓(Appalachian InnerPiedmont)[20]的通道流模式等这类古老造山带的构造演化动力学模型中.

由此可见, 通道流模型的适用与否, 不仅关系到青藏高原隆升的机制问题, 还牵涉到全球大陆动力学系统的合理解释.

在通道流模型的两种端元模型中, Couette流是由顶部或底部刚性较强层的运动所导致的剪切力系牵连效应使中间的通道流层发生流动, 如重力垮塌引起上部刚性上地壳的横向运动, 地幔流动对下地壳的拖拽运动, 以及岩石圈向下俯冲而通过层间耦合驱动软弱的中、下地壳的流动, 其通道流中速度最大的区域发生在顶部或底部; 而Poisemlle流的上、下边界是稳定的, 由两侧的压力梯度差导致通道流层发生运动, 如存在地壳密度差异、地壳厚度差异或是地貌高度差异引起的横向压力梯度驱动壳内软弱层的流动, 其最大流动速度位于通道流层的中部.显然, 通道流的存在必须是低黏度的物质被高黏度的物质所包围, 形成物理意义上的通道.而由于深部地壳真实介质属性难以直接测定, 因而低黏度层的特异黏滞系数主要由数值模拟工作者来界定, 认为只要该层设定为黏滞系数很低的物质组成, 便可产生通道流效应.然而, 必须清晰地认识到, 黏滞系数的正确给定和理想模型的设定是具有一定的局限性, 即必须从实际出发加以约束, 方能提取一个合理的动力学模型.

通道流模型认为低黏度区域的特异黏滞系数往往是由部分熔融或含水流体所致, 而地球内部各圈层的黏滞系数可以通过多要素约束下的边界条件进行理论计算, 亦可通过高温高压实验、地震学、地磁学、重力学和大地测量学等多种途径获得.然而, 用不同手段和方法所给出的数值结果之间其差异明显, 误差范围在最佳情况下亦可达2~3个数量级, 而用固体潮方法得到的黏滞系数差别更高, 甚至可达到1011Pa·s不等(孙和平, 个人交流, 2010).

已开展的数值模拟工作亦存在同样的问题. Clark和Royden[6]的二维模型很好地重现了通道流的演化, 但通道流层内物质的黏滞系数设置过低.同样, Shen等[21]所用三维牛顿黏性流体模拟软弱的下地壳底部时, 对岩石圈地幔的强度设置过高, 这些都难以被深部地球物理工作所完全支持.可见, 复杂的数值模型具有很强的特殊性, 其中的因果关系亦往往不是轻易就能一目了然辨析而出的.

因此, 这样的模拟所得结果仅可供讨论问题或推断及猜测认识时的一个参考, 若以其作为定量的模型或解决深部物质的逼近运移轨迹则尚难以令人信服.

3 青藏高原深部壳、幔物质流展的边界条件和数值模型的构建

近四十年来, 在青藏高原腹地已取得了大量关于地壳和上地幔分层速度、密度和电性结构的结果, 尽管在高原的北部、南部、东部和中部有着一定的异同, 但是在总体上仍存在着一些基本的共同规律.在提取青藏高原的壳、幔介质平均模型时, 我们对人工源地震所得的精细速度结构结果给予了更大的权重, 而以天然地震和大地电磁测深所得结果为辅助约束, 构建了平均模型.从而得到青藏高原地壳与上地幔分层结构的基本模型为(滕吉文, 2011年在中国科学院地质与地球物理研究所学术报告):a)上地壳, 厚约20±5 km, P波速度为5.9±0.1 km/s; b)地壳低速层, 厚约8±2km, P波速度为5. 7±0. 1km/s; c)中、下地壳厚度为40±10km, P波速度为6.7±0. 1km/s; d)上地幔盖层厚度为30±10km, P波速度为8.1±0. 1km/s; e)上地幔软流层顶部埋藏深度为110±10km, P波速度为7. 4±0. 1km/s.

基于青藏高原深部的分层地震波速度结构, 现设青藏高原岩石圈由三层介质组成, 即上地壳、下地壳与岩石圈地幔, 岩石圈下方为软流圈.模型取二维, 深度为150 km, 水平向宽度为300 km.为了便于讨论和明晰起见, 分别考虑到上地壳与下地壳之间的低速层以及软流圈顶面的滑移面, 构建了三个模型进行比较:

模型一:速度结构考虑了两个低速层.a)上地壳厚度为23km, VP=5. 9±0. 1km/s; b)下方23~30 km深度处为低速层, VP=5. 7±0. 1km/s; c) 30~70 km处为下地壳, VP=6. 7±0. 1km/s; d)岩石圈地幔, 即上地幔盖层位于70~97 km深处, VP=8. 1±0. 1km/s; e)软流圈顶部滑移面, 位于深部97~103 km之间, VP=7. 8±0. 1km/s; f)103~150 km之间为软流圈, VP=7. 4±0. 1km/s.

模型二:与模型一相比, 仅考虑上地壳与下地壳之间的低速层, 不考虑软流圈顶部的滑移面, 其分层结构共计5层.a)上地壳厚度为23km, VP=5. 9±0. 1km/s; b)下方23~30 km深度处为低速层, VP=5. 7±0. 1km/s; c)30~70 km处为下地壳, VP=6. 7±0. 1km/s; d)岩石圈地幔, 即上地幔盖层位于70~100 km深处, VP=8. 1±0. 1km/s; e)岩石圈下方, 即100~150 km之间为软流圈, VP=7. 4±0. 1km/s.

模型三:与模型一相比, 两个低速层均不予考虑, 只考虑上地壳、下地壳、岩石圈地幔和软流圈四层.a)上地壳厚度为30km, Vp=5. 9±0. 1km/s; b) 30~70 km处为下地壳, Vp=6. 7±0. 1km/s; c)岩石圈地幔, 即上地幔盖层位于70~100 km深处, Vp=8. 1±0. 1km/s; d)100~150 km之间为软流圈, Vp=7.4±0.1km/s.

为便于计算, 三个模型的边界均为直立或是水平.以上述几何模型和地震波速度为基础, 本项工作数值模型的构建及计算处理均采用有限元软件Ansys, 模型采用了8节点Plane 183的四边形平面单元, 网格密度在上、下之间的低速层及软流圈顶部的滑移面均加密为0. 5km×0. 5 km, 其余层位均设置为1. 0km×1. 0 km, 整个有限元模型单元数为49229, 节点数为49629.模型二和模型三在同一地层中的网格密度与之相同, 但单元网格数及节点数均有所减少.

本文所选取的150 km深度范围内的地球深部介质, 在地质时间尺度下可以简化为黏弹性的本构关系来模拟其力学行为.因而在单元本构关系的构建中, 各圈层均采用了Maxwell体的黏弹性本构关系, 但为了表征不同深度处的圈层结构不同的演化过程, 本文对各圈层介质采用了不同的物理-力学参数.模型一中所采用的力学参数如表 1所示.其中, 弹性模量由密度与弹性波速度之间的关系计算而得, 泊松比均取为0.25.

表 1 模型一中的力学参数 Table 1 The mechanical parameters of materials in the 1st model

在地球动力学数值模拟研究中, 对等效黏滞系数进行合理的估计, 是取得可靠结果的基础. Royden等[13]采用牛顿流体的方法研究了高原东部上、下地壳流变性质差异和地表变形的关系时, 假定地壳浅部(h < 55km)的黏度为1021Pa·s, 地壳深部黏度随深度指数衰减.Shen等[21]利用幂律流体模拟印度的推挤造成40 Ma后青藏高原的抬升高度, 认为上地壳黏度在0. 5×1021~3×1021Pa·s之间.Flesch等[22]利用偏应力和应变速率资料计算了青藏高原及周边地区垂向平均的等效黏滞系数, 结果认为在青藏高原的等效黏滞系数在0.5×1022~5×1022Pa·s之间.

参考先前的研究结果, 黏度结构的选取中, 模型选择青藏高原数值模拟中常用的"三明治"结构, 即具有力学性质强的上地壳和岩石圈地幔, 黏度取1.0×1022 Pa·s, 低速层和滑移面作为软弱层, 黏度取1.0×1020Pa·s.而下地壳和软流圈的力学性质比上地壳和岩石圈地幔弱, 比低速层和滑移面强, 黏度取1.0×1021Pa·s.模型二和模型三中, 相同地层的力学参数选取与模型一相同.

三个模型在计算过程中均选择平面应变模式, 求解器为软件自主选择的稀疏矩阵求解器, 该求解器的计算以直接消元法为基础, 病态矩阵的存在不会引起求解的困难, 因而计算过程具有良好的收敛性, 但对内存及硬盘的配置要求较高.

张培震[23]认为高原东侧相对四川盆地其地壳厚度、地貌高差和垂直形变上均存在横向压力梯度, 但没发现上地壳和岩石圈地幔相对于下地壳的运动, 因而通道流模型更可能是Poiseuille模式, 而非Couette流.Klemperer[9]认为藏南可能是Poiseuille流, 而藏北可能是Poiseuille流与Couette流的混合模式.为此, 本文中的初始模型设置为Poiseuille流的动力模式, 即在左边界施加由高原重力加载随深度变化的静岩压力, 在模型的计算中考虑了重力项所带来的影响, 平均静岩压力简化取为3倍静水压力值.边界条件的选取中, 模型底部施加竖直向约束, 使其水平向能自由移动.模型的右边界施加水平方向的约束, 在竖直方向上自由移动.图 1所示为模型一的卡通图模型, 图中从上到下六种不同颜色的单元层分别代表着上地壳、低速层、下地壳、岩石圈地幔、软流圈顶部滑移面和软流圈六个地层.左边界的箭头代表了静岩压力的方向, 箭头长短表征了其数值大小随深度而变化.

图 1 模型一结构示意图 Fig. 1 Structure sketch of the 1st numerical model

三个模型总的求解时间均为2.5万年, 远大于岩石圈和软流圈介质的黏弹性松弛时间(约数十至数千年), 这就使模型最终计算结果均能进入黏性效应.求解中选择的时间步为1500步, 平均每一步的运算时间步长为16.7a, 小于相应黏度条件下的应力松弛时间, 以确保计算结果的可靠.

4 数值模拟结果对比 4.1 滑脱面的效应

图 2为三种模型所计算得到的位移场分布图, 图 2d为模型一在0.2Ma后的位移场分布图.从图 2a图 2b图 2c三个模型的对比可以看出, 在左边界动力源附近位移值最大, 向右边界方向逐渐减小.同一水平位置所对应的深度上, 模型一的下地壳和岩石圈地幔位移量均比上层的上地壳和下层的软流圈位移量大, 其变化界面正是地壳间低速层和软流圈顶部的滑脱面.同时, 与没有低速层和滑移面的模型三相比, 仅仅2.5万年之后, 模型一的位移场形态便有了与模型三差异明显的变化, 而且相同位置处位移量比模型三增量超过10%.模型二只具有地壳间的低速层, 因而其变形形态及大小介于模型一和模型三之间.可见, 地壳间低速层和软流圈顶部滑移面的共同存在, 对下地壳和岩石圈地幔相对上下层面的快速有序移动起到了明显的增强作用.

图 2 位移场分布图 (a)模型一; (b)模型二; (c)模型三; (d)0. 2 Ma后模型一. Fig. 2 Distribution of displacement field (a) inmodel 1;(b) inmodel2;(c) inmodel3;(d) inmodel 1after 0. 2 Ma.

同时, 从图 2d可见, 0.2Ma之后, 模型一下地壳及岩石圈地幔相对上下层位的位移变化更为明显, 但岩石圈地幔黏度大, 力学强度高, 抵抗变形的能力比下地壳强, 因而位移量弱于下地壳.而上地壳在下地壳快速移动所导致的拖曳作用之下, 其位移情况呈现出与岩石圈地幔近乎同步的变形趋势, 这或许也是导致在这种边界条件下的深部地震波资料解释中SKS横波分裂所表征的各向异性现象不明显的一个因素.

从三个模型里各自取距离左边界100 km处的下地壳和岩石圈地幔中两个单元节点, 观测它们的水平位移量随时间的变化情况, 如图 3所示.图 3a中的节点为下地壳层中距地表深度50 km处, 图 3b中的节点则距地表深度为95 km处, 位于岩石圈地幔地层中部.

图 3 所选节点水平位移与时间关系图 (a)下地壳中; (b)岩石圈地幔中. Fig. 3 Relation graph between the horizontal displacement and the time of the selected nodes (a) in the lower crust; (b) in the upper mantle.

图 3a中, 万年时间尺度之后, 比较三个模型的下地壳水平位移量, 模型三明显小于模型一和模型二, 而模型一和模型二之间下地壳变形的差别一直较小, 表明下地壳顶部低速层对下地壳位移响应的影响效果明显.至于图 3b岩石圈地幔的水平位移量, 三个模型的差别显著地小于它们之间下地壳水平位移的差别量, 这或许是与岩石圈地幔强度较大, 较难变形相关.但模型一和模型二的岩石圈地幔水平位移差别则较它们之间的下地壳水平位移差别明显, 表明模型二中由于缺乏软流圈顶部滑移面的作用, 导致其岩石圈地幔位移能力较模型一为弱.

图 4为三个模型右边界的等效应力值与深度变化的对比, 图上等效应力值的转折处清晰地对应各圈层的界面位置.各圈层内部的等效应力值比较稳定, 而圈层间界面处应力变化明显.可以看到, 软弱层作为滑脱面, 它与上、下圈层界面处是应力差最大的区域, 提供了各界面连接处最大的应力差, 易在此区域形成集中且数值较大的剪切力差, 而软弱层抗剪切能力差, 在有破裂或空隙等情况下, 容易发生剪切破坏现象, 增加了通道内地层物质在压力下有序运移的机率和可能性.下地壳作为力学性质软弱的圈层, 在三个模型中均体现为岩石圈中等效应力强度值最小的区域, 而上地壳和岩石圈地幔均是应力强度值较大的区域, 这与青藏高原地震震源深度基本全在0~20 km的上地壳和70~110 km的岩石圈地幔相符合.

图 4 模型右边界等效应力值 Fig. 4 Sketch of the equivalent effective stress at the right boundary

需要注意的是, 在远离力源作用的远端边界处, 模型二和模型三在Moho面以下的等效应力值趋势基本一致, 而模型一在岩石圈深度以下的软流圈深度范围内, 其等效应力值亦与模型二和模型三基本一致.可见, 通道流效应需要有足够强的初始力源作为保障, 而不像物理学意义上的流体一样, 由于抗剪切能力差, 较小的力源便能发生流动.并且, 需要考虑其作用范围, 即在初始力源作用下, 通道流效应存在着一定的效应范围, 当水平方向位移尺度增加时, 通道流效应逐渐减弱.

4.2 重力及黏度的效应

前人的研究工作也表明, 通道流中软弱层的位移量取决于流体通道的形态、岩石黏度和部分熔融的程度, 以及板块边界的位移速率.在通道壁不平行的情况下, 如果流体黏度足够低, 差应力可以引起通道内高速率的物质浮力折返流动[24].显见, 黏度的选取成为通道流形成的一个重要因素[17].

Ryder等[25]利用Mw7. 6级玛尼地震后4年变形观测得到下地壳黏度为4×1018Pa·s. Bürgmann和Dreson[26]的综述文章里面认为根据震后蠕变变形而估计的黏滞系数可以比1019Pa·s更小.应该注意到, 大地震震后断层附近百公里范围内的变形应变速率在数年内可维持在10-13s-1, 而陆内稳定板块地质时间尺度范围内长期变形的应变速率通常小于或等于10-15s-1, 两者应变率相差可达两个数量级, 同时作用时间也存在数量级方面的差别, 所推算出来的等效黏滞系数亦可以相差20~30倍.

在上文的计算中, 模型黏度参数在可能的取值范围内取值时相互差别较小.为了探讨通道流模型在差别更大的黏度范围内的响应情况, 我们采用前文的几何模型和弹性参数, 重新构建了相对应的模型四、模型五和模型六, 在上文所述前人工作的黏度范围内取值.与模型一相对应, 模型四的力学参数如表 2所示.

表 2 模型四的力学参数 Table 2 The mechanical parameters of materials in the 4 th model

模型依然采用常见的"三明治"结构, 具有力学性质较强的上地壳和岩石圈地幔盖层.软弱的下地壳黏度比用震后变形所得到的黏度约大20倍左右, 低速层和滑移面的黏度比下地壳低一阶, 取值范围均在前人研究的范围内[17].单元的选取和网格的划分方式以及变形时间及时间步的大小不变, 底界面及右边界的位移约束条件同样不变.但在左边界动力学边界条件的施加中, 为了减小各层之间由于密度差所带来的影响, 抵消重力效应, 模型左端除了施加静岩压力外(数值上为三倍静水), 深部每一圈层均增加与上方所有圈层由于密度差所带来的随深度变化的压力差值.这样, 压力从数值上更大地减小了重力对边界力抵消所带来的影响, 并确保了模型的力源约束属于Poiseuille流方式.

图 5为第二组模型在2. 5万年之后的位移场计算结果.图 5a中, 由于模型四地壳间低速层和岩石圈底部滑脱面的存在, 使得下地壳与上地壳之间以及岩石圈地幔与软流圈之间显示位移解耦的现象, 下地壳和软流圈位移明显大于上地壳和岩石圈地幔, 而且越靠近左边界力源差别越显著, 其差异向右边界方向逐渐减小.图 5b中, 模型五由于岩石圈底部与软流圈顶部没有了滑脱面, 在模型中部开始, 岩石圈地幔与软流圈之间的位移已逐渐一致.而低速层与岩石圈地幔之间的位移差也逐渐向右递减, 在模型中部下地壳、岩石圈地幔和软流圈位移趋于一致.而壳内低速层的存在也确保了从下地壳开始, 深部圈层物质的运移比上地壳大.图 5c右端界面处, 四个圈层的位移趋于相等.模型六上地壳的位移明显比模型五小, 在其牵拉作用下, 下地壳顶部和底部位移差别较大, 下地壳顶部Moho面位置成为地壳中位移速度最大的地方, 而这难以从介质属性和机理上解释.

图 5 位移场计算结果 (a)模型四;(b)模型五;(c)模型六. Fig. 5 Displacement field (a) in model 4; (b) in model 5; (c) in model 6.

软弱层的存在, 从力学性质上降低了研究区的等效力学强度, 增大了内部位移及应变变形, 有利于软弱层之间物质的有序流动.

5 结论和讨论

(1)在青藏高原南、北向挤压与东、西方向拉张的构造背景下, 青藏高原腹地和四周地壳与上地幔介质属性和深部物质的运移轨迹清晰地表明:深部物质属性在纵向和横向存在明显不均匀性, 构造展现出南北分区、东西分带的复杂深、浅格局, 并且局部存在有构造变异.在高原深、浅部介质与结构的不均匀性和各向异性作用下, 形成了一系列大型走滑断裂, 促使深部壳、幔物质以这类走滑断裂为边界向东(主体)、西(部分)方向流展, 故又导致了浅部(尤其是高原西部)一系列张性断裂的形成.

青藏高原的通道流模型中[4-7, 9, 13-18, 27], 通常认为通道流的发生层位, 在拉萨块体南部是中地壳, 北部是中、下地壳, 而在青藏高原北部的羌塘和松潘一甘孜块体则发生在中地壳或者是下地壳.虽然中、下地壳在长时间地质年代下呈黏弹性, 低Q值, 物理学意义上存在流动性, 但青藏高原地壳底界的反射波能量强, 体现出其Moho面上部的下地壳属于高地震波速的物质.这种高地震波速的物质, 在地震波场的短时间内, 力学性质可体现为弹性.在固定外力的作用下, 随着时间的增加, 岩石的应变不断增加, 但应变速率逐渐减慢, 最终达到应变以恒定的速率增长, 呈现稳态蠕变现象.这时候, 在应变率特别低的情况下, 应力可以相当低, 但不会为零.而当深部岩石由于外作用力突然降低为零而保持恒定的形变时, 岩石圈内的应力会随时间的增加而逐渐减小, 呈现松弛现象.这类呈黏性流动的蠕变和松弛现象, 说明深部高温岩石在一些方面具有物理学上高黏度流体的性质, 但另一方面, 在低于屈服应力的情况下, 岩石的形变达到一定程度后将保持不变.因而岩石圈的深部岩石还是呈现为固体.

在下地壳和上地幔盖层均为地震波速度随深度增加而递增的高速固态物质的情况下, 即下地壳物质和上地幔盖层之间为固态物质的接触方式[9, 28-33], 这种边界限定了高原下地壳不可能像物理学意义上一样任意流动.

同时, 高原深部整个地壳属于同步增厚, Moho界面被诸多深大断裂所切割, 分别构成了喜马拉雅地体、拉萨地体、羌塘地体、巴颜喀拉地体和昆仑地体, 而且Moho面起伏变化剧烈, 极为凹凸不平, 它不仅不可能构成一个下地壳流动运移的滑移面, 而应体现为阻止下地壳物质沿其界面运动的阻碍面.另外, 若中、下地壳相对于上地壳和上地幔盖层强行运动, 则施加于中地壳或下地壳介质上的力必将足够大, 而在其强迫运移的进程中必当会产生一系列新的物理-力学响应现象, 如高强度的摩擦生热效应, 沿上、下地壳之间的界面和下地壳与上地幔的交界面处亦会有一系列的同源、同深和一定强度的破裂效应以及一系列小地震的发生这种效应会影响这两个交界面上、下及周边介质发生变形和地震波场及速度的变化.然而, 至今均未见这些必然现象的产生[1, 34-35].同时, 通道流模型中, 中、下地壳以塑性变形的方式发生流动, 流动速率随深度增大, 在下地壳处达到最大, 而Moho界面处流动速率最大或相对很大的分布情况, 难以从介质属性和机理上(即高速介质=低黏度介质)得到解释.

由于不具备中、下地壳能够运移的上、下边界条件, 故中、下地壳难以有序运移, 这是任何在给定具有相当误差范围黏滞系数条件下的数学模拟和推断所不能解释的.显然, 青藏高原深部下地壳物质流动, 必须具备一定的物理-力学边界条件和可行的运移轨迹, 否则难以构成该区壳、幔物质的有序流展.

(2)在青藏高原的地壳中, 上地壳底部普遍存在一低速层(低阻层), 它与上地壳为软-硬衔接, 故可解耦而构成壳、幔物质运移的上滑移面, 即提供了一个物质可以移动的自由边界.

下地壳与上地幔盖层物质亦为两坚硬的固态物质相接, 且Moho界面被诸多深大断裂所切割、形成了一个犬牙交错的不规则界面(带), 故不具备可运动的边界条件, 难以发生相对运动.而青藏高原岩石圈底界只能是上地幔低速层的顶部, 在板块构造中一个十分核心的问题就是岩石圈板块漂曳在软流圈上运动, 这亦与岩石圈漂曳在其上运动的边界层相一致, 因而下地壳与上地幔盖层物质只能以上地幔软流层顶部为下滑移面, 才有可能在力系作用下促使其运动, 下滑移面即为软流圈的顶部.

青藏高原上、中地壳之间存在一低速层(深度为20±5km), 而上地幔软流层埋深亦相对较浅(深度为110±10km), 并向高原边缘地域减薄、抬升, 这两个滑脱层面是构成高原内部壳、幔物质能有序运移的唯一边界条件.这种条件下, 不仅下地壳及岩石圈地幔的固态物质较易运移, 而且在上、下边界上不会产生强烈的破裂、变形和热事件.在研究岩石圈深部流变性质时, 多种方法的比较研究是不可缺少的.但在比较中要特别注意的是要了解各种方法估算的前提条件, 明确各种方法的特点和弱点, 以及给出的等效黏滞系数所隐含的使用条件.因此, 这种边界条件的设定满足于青藏高原具体地区的实际观测结果和该边界条件在运动中所导致的效应, 同时不必拘泥于在相当误差条件下对黏度的估计.

而上述的数值模型, 从位移场和等效应力场等的分布和随时间演化情况等几个方面的数值模拟工作证实了软弱层的存在很大程度地增强了高原内部壳、幔物质有序运移的动力学机制.

由此, 青藏高原的岩石圈结构, 在这样的壳、幔分层速度结构为边界条件下, 在两陆-陆板块长期碰撞-挤压作用下, 且高原下地壳与上地幔盖层在其北部和东部均被周边高速、高电阻率和高密度物质所圈闭, 阻隔它们"溢出", 故使地壳短缩增厚.待高原在垂直方向增厚达一定高度后, 深部物质在两陆-陆板块强烈挤压的力系以及高大地形重力势能差的共同作用下, 壳、幔物质运移方式在平面为向东、西方向上伸展, 即在张应力作用下"突破重围", 沿大型走滑断裂系所围通道进行侧向运移.这种运移方式在深部体现为以上、下地壳间的低速层(低阻层)底界为上滑移面, 以上地幔软流层顶部为下滑移面, 岩石圈在地壳低速层与上地壳解耦条件下, 并在驱动力系作用下, 实现"下地壳+上地幔盖层"高速物质在这上、下两个滑移面之间, 沿岩石圈底部, 即软流圈顶部呈同步有序地运动, 而不是下地壳流的力学机制[1, 34-35].

当然, 并非在任何地域的深部都会任意产生这样的物质运移, 只能是局部的和在特异环境下才能实现.因为它必须满足于: 1)壳、幔介质中存在可供高速物质运动的边界条件, 即存在软弱层, 使上、中地壳中的低速层与上地壳解耦, 以及上地幔盖层与软流圈解耦; 2)被移动的介质必须有足够强的力系作用, 即南一北向碰撞、挤压和东一西向拉张; 3)存在主体向东的合理的运移轨迹, 即高原东部大型走滑断层为水平面上深部物质运移的通道边界, 深部存在"下地壳+上地幔盖层"同步运移.

致谢

感谢两位评审专家对本稿提出的宝贵意见和建议.本工作在开展过程中得到了中国科学院青藏高原研究所何建坤研究员和北京大学蔡永恩教授的认真指导, 在此深表感谢.感谢中国南车戚墅堰机车有限公司计算中心陈晓百高工为本工作提供Ansys计算平台, 感谢中国地震局中国台网中心张雪梅副研究员对本工作的支持.

参考文献
[1] 滕吉文, 杨辉, 张雪梅. 中国地球动力学研究的方向和任务. 岩石学报 , 2010, 26(11): 3159–3176. Teng J W, Yang H, Zhang X M. Development direction and task of the geodynamical research in China. Acta Petrologica Sinica (in Chinese) , 2010, 26(11): 3159-3176.
[2] England P C, Houseman G. Finite strain calculations of continental deformation: 2Comparison with the India-Asia collision. J. Geophys. Res. , 1986, 91(B3): 3664-3676. DOI:10.1029/JB091iB03p03664
[3] Flesch L M, Holt W E, Silver P G, et al. Constraining the extent of crust-mantle coupling in central Asia using GPS, geologic, and shear wave splitting data. Earth Planet. Sci. Lett. , 2005, 238(1-2): 248-268. DOI:10.1016/j.epsl.2005.06.023
[4] Zhao W L, Morgan W J P. Injection of Indian crust into Tibetan lower crust: a two-dimensional finite element model study. Tectonics , 1987, 6(4): 489-504. DOI:10.1029/TC006i004p00489
[5] Bird P. Lateral extrusion of lower crust from under high topography in the isostatic limit. J. Geophys. Res. , 1991, 96(B6): 10275-10286. DOI:10.1029/91JB00370
[6] Clark M K, Royden L H. Topographic ooze: building the eastern margin of Tibet by lower crustal flow. Geology , 2000, 28(8): 703-706. DOI:10.1130/0091-7613(2000)28<703:TOBTEM>2.0.CO;2
[7] Beaumont C, Jamieson R A, Nguyen M H, et al. Crustal channel flows: 1. Numerical models with applications to the tectonics of the Himalayan-Tibetan orogen. J. Geophys. Res. , 2004, 109: B06406. DOI:10.1029/2003JB002809
[8] Kuhn T S. The Structure of Scientific Revolutions. Chicago: University of Chicago Press, 1979 .
[9] Klemperer S L. Crustal flow in Tibet: geophysical evidence for the physical state of Tibetan lithosphere, and inferred patterns of active flow.//Law R D, Searle P, Godin L eds. Channel Flow, Ductile Extrusion and Exhumation in Continental Collision Zones. Geological Society, London, Special Publications, 2006, 268: 39-70.
[10] Grujic D, Casey M, Davidson C, et al. Ductile extrusion of the Higher Himalayan Crystalline in Bhutan: evidence from quartz microfabrics. Tectonophysics , 1996, 260(1-3): 21-43. DOI:10.1016/0040-1951(96)00074-1
[11] Searle M P, Rex A J. Thermal model for the Zanskar Himalaya. J. Metamorphic Geol. , 1989, 7(1): 127-134. DOI:10.1111/jmg.1989.7.issue-1
[12] Nelson K D, Zhao W J, Brown L D, et al. Partially molten middle crust beneath southern Tibet: synthesis of project INDEPTH results. Science , 1996, 274(5293): 1684-1688. DOI:10.1126/science.274.5293.1684
[13] Royden L H, Burchfiel B C, King R W, et al. Surface deformation and lower crustal flow in eastern Tibet. Science , 1997, 276(5313): 788-790. DOI:10.1126/science.276.5313.788
[14] Royden L H. Coupling and decoupling of crust and mantle in convergent orogens: implications for strain partitioning in the crust. J. Geophys. Res. , 1996, 101(B8): 17679-17705. DOI:10.1029/96JB00951
[15] Beaumont C, Jamieson R A, Nguyen M H, et al. Himalayan tectonics explained by extrusion of a low-viscosity crustal channel coupled to focused surface denudation. Nature , 2001, 414(6865): 738-742. DOI:10.1038/414738a
[16] Clark M K, Bush J W M, Royden L H. Dynamic topography produced by lower crustal flow against rheological strength heterogeneities bordering the Tibetan Plateau. Geophys. J. Int. , 2005, 162(2): 575-590. DOI:10.1111/gji.2005.162.issue-2
[17] 王晓芳, 何建坤. 下地壳管道流动与青藏高原东缘大尺度构造地貌关系. 中国科学(D辑) , 2012, 55(8): 1338–1390. Wang X F, He J K. Channel flow of the lower crust and its relation to large-scale tectonic geomorphology of the eastern Tibetan Plateau. Science in China (Series D) (in Chinese) , 2012, 55(8): 1338-1390.
[18] Godin L, Grujic D, Law R D, et al. Channel flow, ductile extrusion and exhumation in continental collision zones: an introduction.//Law R D, Searle P, Godin L eds. Channel Flow, Ductile Extrusion and Exhumation in Continental Collision Zones. Geological Society, London, Special Publications, 2006, 268: 1-23.
[19] Xypolias P, Kokkalas S. Heterogeneous ductile deformation along a mid-crustal extruding shear zone: an example from the External Hellenides (Greece).//Law R D, Searle P, Godin L eds. Channel Flow, Ductile Extrusion and Exhumation in Continental Collision Zones. Geological Society, London, Special Publications, 2006, 268: 497-516.
[20] Hatcher R D Jr, Merschat A J. The Appalachian Inner Piedmont: an exhumed strike-parallel, tectonically forced orogenic channel.//Law R D, Searle P, Godin L eds. Channel Flow, Ductile Extrusion and Exhumation in Continental Collision Zones. Geological Society, London, Special Publications, 2006, 268: 517-541.
[21] Shen F, Royden L H, Burchfiel B C. Large-scale crustal deformation of the Tibetan Plateau. J. Geophys. Res. , 2001, 106(B4): 6793-6816. DOI:10.1029/2000JB900389
[22] Flesch L M, Haines A J, Holt W E. Dynamics of the India-Eurasia collision zone. J. Geophys. Res. , 2001, 106(B8): 16435-16460. DOI:10.1029/2001JB000208
[23] 张培震. 青藏高原东缘川西地区的现今构造变形、应变分配与深部动力过程. 中国科学(D辑) , 2008, 38(9): 1041–1056. Zhang P Z. The present day's deformation, strain distribution and deep dynamic process on the Western Sichuan, Eastern Tibetan Plateau. Science in China (Series D) (in Chinese) , 2008, 38(9): 1041-1056.
[24] Gerya T V, Stöckhert B, Perchuk A L. Exhumation of high-pressure metamorphic rocks in a subduction channel: a numerical simulation. Tectonics , 2002, 21(6): 1056. DOI:10.1029/2002TC001406
[25] Ryder I, Parsons B, Wright T J, et al. Post-seismic motion following the 1997 Manyi (Tibet) earthquake: InSAR observations and modelling. Geophys. J. Int. , 2007, 169(3): 1009-1027. DOI:10.1111/gji.2007.169.issue-3
[26] Bürgmann R, Dreson G. Rheology of the lower crust and upper mantle: evidence from rock mechanics, geodesy, and field observations. Annu. Rev. Earth Planet. Sci. , 2008, 36(1): 531-567. DOI:10.1146/annurev.earth.36.031207.124326
[27] Grujic D. Channel flow and continental collision tectonics: an overview.//Law R D, Searle P, Godin L eds. Channel Flow, Ductile Extrusion and Exhumation in Continental Collision Zones. Geological Society, London, Special Publications, 2006, 268: 25-37.
[28] Chen W P, Yang Z H. Earthquakes beneath the Himalayas and Tibet: evidence for strong lithospheric mantle. Science , 2004, 304(5679): 1949-1952. DOI:10.1126/science.1097324
[29] Fan G W, Lay T. Characteristics of Lg attenuation in the Tibetan Plateau. J. Geophys. Res. , 2002, 107: B10, 2256. DOI:10.1029/2001JB000804
[30] Hacker B R, Gnos E, Ratschbacher L, et al. Hot and dry deep crustal xenoliths from Tibet. Science , 2000, 287(5462): 2463-2466. DOI:10.1126/science.287.5462.2463
[31] Kusznir N J, Matthews D H. Deep seismic reflections and the deformational mechanics of the continental lithosphere. J. Petrol. , 1988, 29(1): 63-87.
[32] Ratschbacher L, Frisch W, Liu G H, et al. Distributed deformation in southern and western Tibet during and after the India-Asia collision. J. Geophys. Res. , 1994, 99(B10): 19917-19945. DOI:10.1029/94JB00932
[33] Wernicke B. The fluid crustal layer and its implications for continental dynamics.//Salisbury M H, Fountain D M eds. Exposed Cross-Sections of the Continental Crust, NATO ASI Series, 1990, 317: 509-544.
[34] 滕吉文, 白登海, 杨辉, 等. 2008汶川Ms8.0地震发生的深层过程和动力学响应. 地球物理学报 , 2008, 51(5): 1385–1402. Teng J W, Bai D H, Yang H, et al. Deep processes and dynamic responses associated with the Wenchuan Ms8.0 earthquake of 2008. Chinese J. Geophys. (in Chinese) , 2008, 51(5): 1385-1402.
[35] 滕吉文.地球内部物质与能量交换和动力过程.//21世纪100个交叉科学难题.北京:科学出版社, 2005: 327-338. Teng J W. The exchange of substance and energy and dynamic process in the earth interior.//The 100 Crossing-Subject Problems in 21st Century (in Chinese). Beijing: Science Press, 2005: 327-338.