青藏高原及周边区域地表长期变形数值模拟

董培育　胡才博　石耀霖*

1)中国地震局地震研究所， 地震大地测量重点实验室， 武汉　430071 2)中国科学院计算地球动力学重点实验室， 北京　100049



青藏高原及周边区域地表长期变形数值模拟

董培育1，2)胡才博2)石耀霖2)*

1)中国地震局地震研究所， 地震大地测量重点实验室， 武汉430071 2)中国科学院计算地球动力学重点实验室， 北京100049

印度板块向欧亚板块俯冲挤压， 不仅令青藏高原上地壳在挤压作用下发生弹性变形和运动， 且青藏高原高温高压下的下地壳会发生柔性流动， 并对脆性的上地壳有拖曳作用， 这2种作用一起形成现今的高原运动变形场。这一动力学过程已得到GPS观测资料的证实。因此在二维平面问题中仅用上地壳在边界作用下的弹性变形解释是不够的， 还要考虑柔性下地壳流动对上地壳的拖曳作用。但是拖曳力作用的大小和方向不易确定， 故文中建立了二维平面弹性有限元模型， 利用加载等效体力来模拟下地壳流动对上地壳产生的拖曳力。以高原内部的GPS观测资料为约束， 利用试错法反演出模型中关键点的力， 其他位置上的力则用关键点上的力进行双线性插值计算。以此来反演计算出模型区域内的柔性下地壳的差异性流动对脆性上地壳产生的拖曳力(节点力的形式， 单位： N)的大小和范围， 在86°～100°E， 26°～32°N地区主要以SE向为主， 最大达到108N； 西部局部(31°～36°N， 76°～80°E)地区有较弱的W向拖曳力， 最大为107N。文中为深入研究青藏高原及周边区域的长期地表变形动力学机制提供了1个新的思路。

有限元数值模拟青藏高原拖曳力双线性插值

0　引言

印度板块长期持续向N挤压欧亚板块， 形成了青藏高原， 使其成为目前地球上最年轻、 形成过程最复杂的高原， 令其一直处于地球科学研究和讨论的中心和热点(许志琴等， 2011)。对于其隆升和运动变形机制， Tapponnier等(1982)认为高原被大型走滑断裂带切割， 在印度板块推挤下， 各块体沿着断裂带被挤出。England等(1982， 1983)利用薄板流变模型数值模拟印度板块与欧亚板块的碰撞形变， 且在一定区域中定量描述了高原地壳厚度和应变率的时空分布。傅容珊等(1999； 2000a， b)将大陆岩石层视为幂指数律控制的薄层， 上覆在黏滞性较低的软流层之上蠕变流动； 数值模拟结果表明， 青藏高原隆升的主要动力来源于印度板块对欧亚板块的碰撞挤压。许志琴等(1996)提出青藏高原腹地隆升的主要原因与大陆内部深部热驱动——地幔底辟有关。目前得到众多学者普遍认可的主要是， 上地壳为脆性物质， 下地壳为柔性物质， 且温度高， 密度低， 导电性好， 地震波速低等， 因此其黏滞系数较低， 非常柔软。 在印度板块挤压下， 下地壳容易流动， 拖曳上地壳运动， 导致高原的隆升和变形(吴功建等， 1991； 石耀霖等， 1992； 曹建玲等， 2009)。

随着GPS等空间大地测量技术日渐成熟， 全球有越来越多的GPS观测台站投入使用。中国在青藏高原上的GPS观测台站观测到了大量的数据， 为研究青藏高原的运动和变形提供了重要的数据支持。 如图2 所示为GPS观测值， 这些数据揭示了青藏高原水平运动的主要特征。 其SN方向上， 主要向N运动， 且从南部到北部， 运动速率逐渐衰减， 显示出SN向在缩短(Ganetal.， 2007)， 是由地壳增厚形变吸收掉了绝大部分， 另外小部分由高原内部的一系列大型走滑断裂带吸收了(任金卫等， 2003； 沈正康等， 2003； 许志琴等， 2011)。 EW方向上， 主要显示为EW向的张开， 即东部向E运动， 西部向W运动， 符合物质被挤压后向两侧移动的力学特征； 但是W向的运动相对微弱， 主要为E向的运动， 且在东南缘， 伴随有S向的运动， 显示其为顺时针方向运动。这是因为在青藏高原东部有坚硬的四川盆地阻挡， 而其东南部云南到缅甸地区恰好有1个软弱的利于岩石流动的通道， 所以物质往S流动， 青藏高原整体上形成了顺时针旋转的形态(朱守彪等， 2004； 曹建玲等， 2009)。 在青藏高原西北部为帕米尔高原， 在深部有部分物质俯冲至帕米尔高原下， 西侧物质向W挤出， 但是挤出量相对微弱(雷建设等， 2002； 张家声等， 2005； 许志琴等， 2011)。

由于在青藏高原上由GPS观测数据揭示的长期地表变形表现出如此复杂多变的特点， 仅仅用上地壳在边界作用下的弹性变形解释是不够的， 还要考虑柔软的下地壳流动拖曳上地壳运动的附加作用(曹建玲等， 2009)。那么拖曳力有多大呢？方向如何？是如何分布的？石耀霖等(2000)首先提出了在二维模拟中引入等效体力部分反映某些三维作用的方法。朱守彪等(2005)利用遗传有限单元算法反演出青藏高原受到的边界作用力、 地形扩展力以及下地壳对上地壳的拖曳力。曹建玲等(2009)指出下地壳对上地壳的拖曳力依赖于下地壳的黏滞系数， 但是黏滞系数的确切给定尚有困难。本文将建立二维平面弹性有限元模型， 利用GPS观测资料作为约束条件， 选取关键点， 利用试错法给出关键点上的力， 同时双线性插值计算其他位置上的力， 以此为基础正演模拟得到的长期地表变形与GPS观测资料吻合较好。

1　模型及边界条件

1.1区域模型

本文参考曹建玲等(2009)的研究， 主要计算青藏高原下地壳对上地壳的拖曳作用， 因此选取较简单模型， 区域范围为75°～110°E， 25°～41°N。据前人研究(Tapponnier， 2001; 张培震等， 2003； 王辉等， 2006； Ganetal.， 2007； 闻学泽等， 2011； 徐锡伟等， 2014)， 青藏高原及周边区域活动地块主要分为拉萨、 羌塘、 巴颜喀拉、 柴达木、 祁连、 川滇、 滇南、 塔里木等， 各个地块之间主要的断裂带有喜马拉雅俯冲带、 阿尔金断裂带、 祁连-海原断裂带、 昆仑山断裂带、 甘孜-玛尼-玉树断裂带、 龙门山断裂带、 鲜水河断裂带、 嘉黎断裂带等。由于一些大的断裂带周围有很多小断裂带， 我们将模型简化， 小断裂带整合为大的断裂带， 除了喜马拉雅俯冲带为20km宽外， 其他断裂带均为10km宽。由于本文为二维平面模型， 初步模型采用笛卡尔坐标系， 将研究区的经纬度坐标系投影转换到直角坐标系进行计算。模型及剖分网格如图1 所示， 采用三角形网格， 在断裂带地区加密， 共有157,649个节点， 310,439个单元。

图1　模型及网格Fig. 1　Model and grid.

表1　介质参数表

Table1　Medium parameters

杨氏模量/Pa泊松比塔里木10×10100.25柴达木9×10100.25巴颜喀拉9×10100.25羌塘9×10100.25拉萨9×10100.25雅鲁藏布江地块2×10100.22喜马拉雅俯冲带南部(印度板块)1.5×10100.22阿尔金断裂带8×10100.25祁连-海原断裂带6×10100.24昆仑山断裂带5.5×10100.23甘孜-玛尼-玉树断裂带5×10100.22龙门山断裂带5.0×10100.22龙门山北部断裂带5.5×10100.23鲜水河断裂带5.0×10100.22理塘断裂带5.5×10100.23嘉黎断裂带(西)7.5×10100.25嘉黎断裂带(东)7.0×10100.24喀喇昆仑断裂带(西)6.5×10100.24喀喇昆仑断裂带(东)8.0×10100.25喜马拉雅俯冲带(西)8×1090.22喜马拉雅俯冲带(东)1.0×10100.22

1.2材料参数

参考前人的研究成果(Huangetal.， 2003； 王辉等， 2006； 黄建平等， 2008； 陈连旺等， 2011)， 并且考虑到各断裂带内部较软弱， 给出相对于地块较低的杨氏模量E。 在最近几十年来较活跃的， 发生过地震的断裂带相对弱一些， 而较稳定的断裂带， 稍硬一些。由于印度板块挤压欧亚板块， 青藏高原南缘抬升， 因此藏南地区软一些， 包括喜马拉雅俯冲带南部的印度板块也较软， 而北部有坚硬的塔里木盆地， 所以藏北地区硬一些。综合考虑各方面因素， 给出各个地块和断裂带的杨氏模量和泊松比， 如表1 所示。

1.3边界条件

模型采用青藏高原及周边的726个观测台站上观测得到的相对于稳定的欧亚板块的GPS水平观测数据(Ganetal.， 2007)(特别指出， 本文给定的GPS速度场及模拟计算结果均为相对于稳定的欧亚参考框架)， 如图2 中的黑线所示； 并将其均匀插值到模型边界上， 得到有限元模型的边界条件， 如图2 中的红线所示。

图2　GPS台站观测数据以及插值边界条件Fig. 2　GPS observation data and boundary condition via interpolation.

1.4初步结果

在只有GPS边界条件的情况下， 不考虑下地壳对上地壳的拖曳力作用， 利用以上有限元模型计算研究区域的长期速度场分布， 将初步计算结果与GPS观测值之间进行对比， 如图3 所示， 主要差异为：

(1)在青藏高原东部， 观测值明显大于计算值。观测值显示地表变形有强烈的SE向旋转， 而实际计算结果只能显示非常微弱的旋转， 甚至没有旋转。王辉等(2006)的青藏高原二维数值模拟结果显示， 只有观测的GPS速率边界条件， 需要在东南边施加大于观测值的S向位移， 才能够跟观测值吻合。曹建玲等(2009)利用三维流变模型模拟青藏高原GPS位移， 指出在东南缘加载应力边界条件， 结果与观测值之间的误差较小。朱守彪等(2004)对川滇地区的计算结果表明， 加载拖曳力后， 能够与GPS观测值大体吻合。

(2)在西部地区， 观测值的相对计算结果有微弱的W向偏转。朱守彪等(2005， 2006)的研究结果指出青藏高原的西部地区有W向的拖曳力。黄建平等(2008)计算出的地幔对流产生的拖曳力场显示， 在青藏高原西部地区存在有W向的拖曳力。

图3　仅加载边界条件计算结果与观测值对比Fig. 3　Comparison between observation and simulation results with only boundary condition.

观测结果和模拟计算结果对比表明， 仅仅采用GPS边界条件得到的计算结果不能够与实际观测值吻合， 因为没有考虑柔性的下地壳对上地壳的拖曳力作用。Copley等(2007)模拟了青藏高原地表的GPS速度， 指出主要是下地壳流动， 拖曳上地壳运动， 导致该区域的构造变形。曹建玲等(2009)的研究结果指出GPS绕喜马拉雅东构造结的顺时针旋转速度大于GPS观测到的地表速度， 该速度之差每年达到数毫米； 在柔性下地壳区域， 下地壳对上地壳有拖曳力作用， 大小一般在0.01～0.1MPa之间。以上结果表明， 在模拟青藏高原及其邻区长期变形时， 下地壳对上地壳的拖曳力是不可以忽略的。

2　拖曳力计算

下地壳对上地壳的拖曳力如何计算呢？采用朱守彪等(2004)文章中提出的方法， 在弹性学平衡方程中令水平剪应力分量随深度变化为常量， 那么在二维问题中， 该部分即可作为等效体力， 代表特定深度内下地壳对上地壳的拖曳剪切力作用。如式(1)所示：

(1)

考虑在上1节的有限元模型中加载拖曳力； 除此之外， 几何模型、 材料模型和边界条件均不变。由于模型是二维平面的， 拖曳力在有限元计算中以节点力或集中力的形式出现。我们采用试错法调整计算关键点上的拖曳力。关键点的选取如图4b所示， 共26个。

图4 为区域划分简图， 4个子区域组成1个大区域。如图4a中所示以1号点为中心点， 共有4个子区域， 且1号、 2号、 3号和4号， 4个点组成的是其中1个子区域。每1个子区域(由4个点组成)作为1个双线性插值区， 有4个角点的拖曳力值， 就可以插值求得整个子区域内的其他任意点上的拖曳力。同样地， 可以得到其他3个子区域上的值。这样就得到了整个大区域上的拖曳力值。参考该示意图， 我们将模型划分区域， 主要观察哪些地区的观测值和模拟值差异大， 以此作为分区依据。

图4　拖曳力区域划分图Fig. 4　Sketch of region division.a 分区简单示意图； b 真实分区图， 圆圈内的数字表示子区域号， 方框内的数字表示关键点号

将每个大区域内的每1个子区域(即4个点组成的区域)作为1个双线性插值区， 区域内各个点的节点力表示为

(2)

式(2)中的a、b、c、d、m、n、o、p由下面的方程组(3)解出

(3)

式(3)中的x1、 y1、 x2、 y2、 x3、 y3、 x4、 y4分别为4个点的坐标， Fx1、 Fx2、 Fx3、 Fx4和Fy1、 Fy2、 Fy3、 Fy4分别为根据观察和参考前人工作给出的试探节点力值(x为EW向，y为SN向)， 在试验中不断改进。然后根据式(1)即可得到各个区域内任意1个节点上的力。最终得到的区域内关键点的拖曳力值如表2 所示。插值到模型区域内各节点上的拖曳力如图5 所示。

表2　各个子区域内关键点的拖曳力值(子区域内关键点的排列顺序如图4a所示)

Table2　Drag force in different key nodes of each subregion (the order of key nodes in each subregion is shown as in Fig. 4a)

N

图5　加载至模型中的拖曳力Fig. 5　Drag force in this model.

如图6 所示， 加载拖曳力后有限元模型计算得到的速度场值与GPS观测值之间吻合较好， 只有个别地区仍存在误差。这已经是经过多次试验得到的最好的结果。 按照式(4)计算均方根误差， 式中的n是所有加载拖曳力的节点的个数，Δi是每个点的计算值和观测值之差，error是均方根误差。结果如表3 所示， 加载拖曳力之后的均方根误差明显小于无拖曳力的情况。

(4)

图6　边界条件和拖曳力共同作用的计算结果与台站观测对比Fig. 6　Comparison of observation and simulation results with boundary condition and drag force.a 观测台站上的计算结果与观测值对比； b 均匀插值分布的结果对比

表3　计算结果与观测值均方根误差对比

Table3　Comparison of MSR error between simulations and observations

EW向分量/mm·a-1SN向分量/mm·a-1矢量模/mm·a-1均匀插值无拖曳力5.335,12.434,45.773,4有拖曳力1.711,72.140,32.685,2观测台站无拖曳力4.658,63.097,95.281,9有拖曳力2.194,22.565,43.153,6

图7　计算的应变率场Fig. 7　Strain rate field by calculation.

表4　主要断裂带的走向、 主应变率方向及滑动方向

Table4　Strike and direction of the principle strain rate and slip of major faults

断裂带主压应变率方向主张应变率方向二者比较滑动方向阿尔金断裂带NE-SWNW-SE相当左旋海原断裂带NE-SWNW-SE压大于拉左旋昆仑山断裂带NE-SWNW-SE相当左旋甘孜玉树断裂带NE-SWNW-SE相当左旋鲜水河断裂带NWW-SEENNE-SSW相当左旋龙门山断裂带NWW-SEENNE-SSW相当右旋(逆掩不明显)嘉黎断裂带NNE-SSWNWW-SEE压稍小于拉右旋喀喇昆仑-雅鲁藏布江断裂带近SN近EW压远大于拉右旋喜马拉雅断裂带近SN近EW压远大于拉逆掩

计算出的应变率场如图7 所示， 与前人文中的应变率场基本一致， 具有可比性(沈正康等， 2003； Ganetal.， 2007)。图7 显示喜马拉雅处于SN向挤压状态， 对应逆掩断层； 藏南地区EW向拉应变率明显大于SN向压应变率， 与该区SN走向的正断层吻合； 高原内部的数条走滑断层， 其右旋或者左旋滑动性质， 计算与观测吻合。主要断裂带的详细情况如表4 所示。

3　结论与讨论

模拟计算依赖于材料参数和边界条件， 需要参考更多的资料对材料参数的差异性分布进行约束。本文依靠GPS数据插值得到边界条件， 而研究区域的GPS台站数不足， 尤其在青藏高原中部地区台站较少， 插值得到的数据精度会受到很大的影响。模型西北部即阿尔金断裂带和喀喇昆仑断裂带交界的区域构造非常复杂， 在帕米尔高原还发生着双向俯冲(许志琴等， 2011)， 还需要进一步细化研究。

另外， 地球曲率的影响不容忽视， 石耀霖等(2006)指出在利用位移计算应变时， 如果采用平面坐标， 尤其在高纬度地区， 会产生较大的系统误差。由于本文为初步模型， 暂且考虑直角坐标系。同时， 虽然本文采用的二维平面模型能够将问题简化， 但是不能考虑真实的三维结构， 如下地壳黏滞系数， 不能与垂向的GPS实测值进行对比等。因此在深入计算中， 需考虑三维球坐标， 以及更加精细的物质结构和边界条件等。

本文采用加载等效体力的方式来模拟柔性下地壳对脆性上地壳的拖曳剪切力作用。选取关键点， 给出试探值， 利用试错法， 给出最优结果， 同时在其他位置点上利用双线性插值计算其等效体力。最终得到区域内的上地壳下部受到的拖曳力的大小和方向， 以此为基础模拟得到的速度场结果能够与GPS观测有较好的吻合， 误差较小， 为深入研究青藏高原及其邻区地表变形提供了新的思路。

曹建玲， 石耀霖， 张怀， 等. 2009. 青藏高原GPS位移绕喜马拉雅东构造顺时针旋转成因的数值模拟 [J]. 科学通报， 54(8): 1398—1410.

CAO Jian-ling， SHI Yao-lin， ZHANG Huai，etal. 2009. Numerical simulation of GPS observed clockwise rotation around the eastern Himalayan syntax in the Tibetan plateau [J]. Chinese Science Bulletin， 54(8): 1398—1410(in Chinese).

陈连旺， 詹自敏， 叶际阳， 等. 2011. 流变特性对青藏高原构造变形影响的数值模拟 [J]. 大地测量与地球动力学， 31(3): 8—14.

CHEN Lian-wang， ZHAN Zi-min， YE Ji-yang，etal. 2011. Numerically modeling the influence of rheological properties on tectonic deformation of Tibet Plateau [J]. Journal of Geodesy and Geodynamics， 31(3): 8—14(in Chinese).

傅容珊， 黄建华， 徐耀民. 2000a. 印度与欧亚板块碰撞的数值模拟和现代中国大陆形变 [J]. 地震学报， 22(1): 1—7.

FU Rong-shan， HUANG Jian-hua， XU Yao-min. 2000a. Numerical simulation of the Inida-Asia collision and deformation of China Mainland [J]. Acta Seismologica Sinica， 22(1): 1—7(in Chinese).

傅容珊， 李力刚， 黄建华， 等. 1999. 青藏高原隆升过程的三阶段模式 [J]. 地球物理学报， 42(5): 609— 616.

FU Rong-shan， LI Li-gang， HUANG Jian-hua，etal. 1999. Three-step model of the Qinghai-Xizang Plateau uplift [J]. Chinese Journal of Geophysics， 42(5): 609— 616(in Chinese).

傅容珊， 徐耀民， 黄建华， 等. 2000b. 青藏高原挤压隆升过程的数值模拟 [J]. 地球物理学报， 43(3): 346—355.

FU Rong-shan， XU Yao-min， HUANG Jian-hua，etal. 2000b. Numerical simulation of the compression uplift of the Qinghai-Xizang plateau [J]. Chinese Journal of Geophysics， 43(3): 346—355(in Chinese).

黄建平， 傅容珊， 郑勇， 等. 2008. 地幔对流拖曳力对中国大陆岩石层变形的影响 [J]. 地球物理学报， 51(4): 1048—1057.

HUANG Jian-ping， FU Rong-shan， ZHENG Yong，etal. 2008. The influence of mantle convection to the lithosphere deformation of China mainland [J]. Chinese Journal of Geophysics， 51(4): 1048—1057(in Chinese).

雷建设， 周蕙兰， 赵大鹏. 2002. 帕米尔及邻区地壳上地幔P波三维速度结构的研究 [J]. 地球物理学报， 45(6): 802—811.

LEI Jian-she， ZHOU Hui-lan， ZHAO Da-peng. 2002. 3-D velocity structure of P-wave in the crust and upper-mantle beneath Pamir and adjacent region [J]. Chinese Journal of Geophysics， 45(6): 802—811(in Chinese).

任金卫， 马宗晋. 2003. 东亚地区现代地壳运动特征与构造变形 [J]. 地学前缘， 10(特刊): 58— 65.

REN Jin-wei， MA Zong-jin. 2003. Crustal movement and tectonic deformation of eastern Asia [J]. Earth Science Frontiers， 10: 58— 65(in Chinese).

沈正康， 王敏， 甘卫军， 等. 2003. 中国大陆先进构造应变率场及其动力学成因研究 [J]. 地学前缘， 10(特刊): 93—100.

SHEN Zheng-kang， WANG Min， GAN Wei-jun，etal. 2003. Contemporary tectonic strain rate field of Chinese continent and its geodynamic implications [J]. Earth Science Frontiers， 10: 93—100(in Chinese).

石耀霖， Assumpcao M. 2000. 巴西构造应力场的遗传有限单元法反演 [J]. 地球物理学报， 2000， 43(2): 166—174.

SHI Yao-lin， Assumpcao M. 2000. Genetic algorithm-finite element inversion of stress field of Brazil [J]. Chinese Journal of Geophysics， 43(2): 166—174(in Chinese).

石耀霖, 朱守彪. 2006. 用GPS位移资料计算应变方法的讨论 [J]. 大地测量与地球动力学, 26(1) : 1—8.

SHI Yao-lin, ZHU Shou-biao. 2006. Discussion on method of calculating strain with GPS displacement data [J]. Journal of Geodesy and Geodynamics, 26(1): 1—8(in Chinese).

石耀霖， 朱元清， 沈显杰. 1992. 青藏高原构造热演化的主要控制因素 [J]. 地球物理学报， 35(6): 710—720.

SHI Yao-lin， ZHU Yuan-qing， SHEN Xian-jie. 1992. Tectonic processes and thermal evolution of the Qinghai-Xizang(Tibetan)Plateau [J]. Chinese Journal of Geophysics， 35(6): 710—720(in Chinese).

王辉， 张国民， 石耀霖， 等. 2006. 青藏活动地块区运动与变形特征的数值模拟 [J]. 大地测量与地球动力学， 26(2): 15—23.

WANG Hui， ZHANG Guo-min， SHI Yao-lin，etal. 2006. Numerical simulation of movement and deformation of Qinghai-Tibet Plateau [J]. Journal of Geodesy and Geodynamics， 26(2): 15—23(in Chinese).

闻学泽， 杜方， 张培震. 2011. 巴颜喀拉块体北和东边界大地震序列的关联性与2008年汶川地震 [J]. 地球物理学报， 54(3): 706—716.

WEN Xue-ze， DU Fang， ZHANG Pei-zhen. 2011. Correlation of major earthquake sequences on the northern and eastern boundaries of the Bayan Har block， and its relation to the 2008 Wenchuan earthquake [J]. Chinese Journal of Geophysics， 54(3): 706—716(in Chinese).

吴功建， 高锐， 余钦范. 1991. 青藏高原亚东—格尔木地学断面综合地球物理调查与研究 [J]. 地球物理学报， 34(5): 552—560.

WU Gong-jian， GAO Rui， YU Qin-fan，etal. 1991. Integrated investigations of the Qinghai-Tibet Plateau along the Yadong-Golmud geoscience transect [J]. Chinese Journal of Geophysics， 34(5): 552—562(in Chinese).

徐锡伟， 程佳， 许冲， 等. 2014. 青藏高原块体运动模型与地震活动主体地区讨论： 鲁甸和景谷地震的启示 [J]. 地震地质， 36(4): 1116—1134.

XU Xi-wei， CHENG Jia， XU Cong，etal. 2014. Discussion on block kinematic model and future themed areas for earthquake occurrence in the Tibetan plateau: Inspiration from the Ludian and Jinggu earthquakes [J]. Seismology and Geology， 36(4): 1116—1134(in Chinese).

许志琴， 姜枚， 杨经绥. 1996. 青藏高原北部隆升的深部构造物理作用： 以 “格尔木—唐古拉山”地质及地球物理综合剖面为例 [J]. 地质学报， 70(3): 195—206.

XU Zhi-qin， JIANG Mei， YANG Jing-sui. 1996. Tectonophysical process at depth for the uplift of the northern part of the Qinghai-Tibet Plateau: Illustrated by geological and geophysical comprehensive profile from Golmud to the Tanggula Mountains， Qinghai Province [J].Acta Geologica Sinica， 70(3): 195—206(in Chinese).

许志琴， 杨经绥， 李海兵， 等. 2011. 印度-亚洲碰撞大地构造 [J]. 地质学报， 85(1): 1—33.

XU Zhi-qin， YANG Jing-sui， LI Hai-bing，etal. 2011. On the tectonics of the India-Asia collision [J]. Acta Geologica Sinica， 85(1): 1—33(in Chinese).

张家声， 单新建， 李建华， 等. 2005. 帕米尔地区现今大陆深俯冲： 地震构造和动力学解释 [J]. 岩石学报， 26(4): 1215—1241.

ZHANG Jia-sheng， SHAN Xin-jian， LI Jian-hua，etal. 2005. Recent deep subducting of continental crust in Pamier: An interpretation on seismotectonics and geodynamics [J]. Acta Petrologica Sinica， 26(4): 1215—1241(in Chinese).

张培震， 邓起东， 张国民， 等. 2003. 中国大陆的强震活动与活动地块 [J]. 中国科学(D辑)， 33(增刊): 12—20.

ZHANG Pei-zhen， DENG Qi-dong， ZHANG Guo-min，etal. 2003. Strong earthquake activity and active block in China mainland [J]. Science in China(Ser D)， 33: 12—20(in Chinese).

朱守彪， 石耀霖. 2004. 用遗传有限单元法反演川滇下地壳流动对上地壳的拖曳作用 [J]. 地球物理学报， 47(2): 232—239.

ZHU Shou-biao， SHI Yao-lin. 2004. Genetic algorithm-finite element inversion of drag forces exerted by the lower crust on the upper crust in the Sichuan-Yunnan area [J]. Chinese Journal of Geophysical， 47(2)： 232—239(in Chinese).

朱守彪， 石耀霖. 2005. 青藏高原地形扩展力以及下地壳对上地壳拖曳力的遗传有限单元法反演 [J]. 北京大学学报(自然科学版)， 41(2): 225—234.

ZHU Shou-biao， SHI Yao-lin. 2005. Genetic algorithm-finite element inversion of topographic spreading forces and drag forces of lower crust to upper crust in Tibetan plateau [J]. Acta Scientiarum Naturalium Universitatis Pekinensis， 41(2): 225—234(in Chinese).

朱守彪， 石耀霖. 2006. 中国大陆及邻区构造应力场成因的研究 [J]. 中国科学(D辑)， 36(12): 1077—1083.

ZHU Shou-biao， SHI Yao-lin. 2006. Causes of tectonic stress field in China and adjacent areas [J]. Science in China(Ser D)， 36(12): 1077—1083(in Chinese).

Copley A， Mckenzie D. 2007. Models of crustal flow in the Inida-Asia collision zone [J]. Geophysical Journal International， 169(2): 683— 698.

England P， Mckenzie D. 1982. A thin viscous sheet model for continental deformation [J]. Geophysical Journal International， 70(2): 295—321.

England P， Mckenzie D. 1983. Correction to： A thin viscous sheet model for continental deformation [J]. Geophysical Journal International， 73(2): 523—532.

Gan W J， Zhang P Z， Shen Z K，etal. 2007. Present-day crustal motion within the Tibetan plateau inferred from GPS measurements [J]. Journal of Geophysical Research： Solid Earth(1978—2012)， 112(B8).

Huang Z X， Su W， Peng Y J，etal. 2003. Rayleigh wave tomography of China and adjacent regions [J]. Journal of Geophysical Research: Solid Earth(1978—2012)， 108(B2).

Tapponnier P. 2001. Oblique stepwise rise and growth of the Tibet plateau [J]. Science， 294(5547): 1671—1677.

Tapponnier P， Peltzer G， Dain A Y，etal. 1982. Propagating extrusion tectonics in Asia: New insights from simple experiments with plasticine [J]. Geology， 10(2): 611— 616.

Abstract

The subduction of the Indian plate underneath Eurasian plate results not only in deformation and movement of the elastic upper crust， but also flow of the ductile lower crust in the high temperature and high pressure which drags the brittle upper crust to move at the same time. These two actions work together producing the present movement and deformation field in Tibetan plateau. The dynamics progress has been verified by GPS observation data. Therefore， in a two-dimension plain model， only the elastic deformation with the boundary action at the upper crust cannot explain the deformation well， the drag force acted on the base of upper crust by the drag of ductile flow of the lower crust also need to be considered. However， it’s hard to figure out the magnitude and direction of the drag force. Thus， we established a two-dimension plain elastic finite element model， with the equivalent-body force approach to simulate the drag force. With the internal GPS observation data of Tibetan plateau as constraint condition， we calculated inversely the drag force of key nodes in the model with trial method， and the other nodes in the model with bilinear interpolation method. Finally， we got the drag forces(nodal forces， unit: N)caused by the difference flow of ductile lower crust dragging the brittle upper crust， which are distributed mainly in the region of 86°～100°E and 26°～32°N， the direction is east and south， and the maximum reaches to 1e8N; in some areas in the western part of the study region at 31°～36°N and 76°～80°E， the direction is west， and the maximum reaches to 1e7N. All these work provides a new thought for further research on long-term dynamic mechanism of surface deformation in Tibetan plateau and its surrounding area.

NUMERICAL SIMULATION OF LONG-TERM DEFORMATION OF TIBETAN PLATEAU AND SURROUNDING AREA

DONG Pei-yu1，2)HU Cai-bo2)SHI Yao-lin2)

1)InstituteofSeismology，ChinaEarthquakeAdministration，Wuhan430071，China2)KeyLaboratoryofComputationalGeodynamics，CAS，UniversityofChineseAcademyofSciences，Beijing100049，China

finite element method， numerical simulation， Tibetan plateau， drag forces， bilinear interpolation

2015-01-19收稿， 2015-12-10改回。

中国地震局地震研究所所长基金项目(IS201526237)、 国家科技支撑项目“地震预报实用技术”(2012BAK19B035)、 深部探测数据集成与共享服务(201511028)与国家自然科学基金(41474085， 41274027)共同资助。
*

石耀霖， 男， 教授， E-mail: shiyl@ucas.ac.cn。

P315.2

A

0253-4967(2016)02-410-13

董培育， 女， 1987年生， 2015年毕业于中国科学院大学地球科学学院固体地球物理学专业， 获博士学位， 助理研究员， 主要研究方向为地球动力学问题的数值模拟， E-mail: dongpeiyu97@163.com。

doi:10.3969/j.issn.0253- 4967.2016.02.014