天山地区现今地壳形变及构造应力的三维数值模拟*

2022-09-01孙明志鲁小飞张彩红刘志军

地震研究 2022年4期

孙明志，谭凯，鲁小飞，张彩红，李琦，刘志军

(中国地震局地震大地测量重点实验室，湖北武汉 430071)

0 引言

天山山系绵延于中亚腹地，全长达2 500 km，南北宽250～400 km，是世界上最大、最活跃的陆内造山带(邓起东等，2000；Vinnik，2002)。新生代以来，在印度与欧亚板块的碰撞及两侧盆地楔入背景下，天山山体复活隆升，两侧盆地凹陷，发生了强烈的地壳变形，发育了一系列断裂带，控制了一系列强震(Molnar，Tapponnier，1975；张培震等，1996)。研究天山地区的地壳形变与应力分布特征，对理解该地区构造变形的运动学机制及地震成因，分析未来潜在地震危险性有重要的科学意义。

大量学者利用GPS等大地测量资料，采用插值方法、弹性位错模型等对天山地区的构造运动特征进行了研究，得到了一系列十分有价值的研究成果(王琪等，2000；牛之俊等，2007；王晓强等，2007；杨少敏等，2008；王伟等，2014)。随着计算机技术和数值计算理论的发展，数值模拟逐渐成为了研究地学问题的重要手段。与传统研究方法相比，数值模拟体现出一定的计算优势。有限元数值模拟以较为客观的模型架构和各种地质、地震及大地测量资料为基础，可以兼顾整个区域横向及纵向上物质的非均匀性，由此可以得到更为准确、细致的构造运动特征(刘峡，2007)。

近年来不少学者对天山地区的构造变形特征进行了数值模拟研究，如Qiu等(2009)利用有限元模拟西南天山与帕米尔—西昆仑之间的汇聚趋势，得出了汇聚边界向北、向东南延伸的结论，但其模型没有考虑断裂的影响；雷显权等(2011)采用非连续接触分析方法模拟天山地区的地壳运动变形，发现地壳运动速度由西向东、由南向北逐渐减小，然而其建立的天山中段二维弹性薄片模型，不能考虑复杂介质分布情况，在空间上具有一定的局限性；王子韬等(2020)通过有限元计算了天山区域主要断裂带弹性应变能密度增加速率及应力积累速率，发现天山的断层和地震活动性主要受控于近南北向的主压应力，其建立的二维平面弹性模型将断层处理为连续变形的软弱带，不能很好地体现断裂带的滑动特征。

物理模型与地球实际情况的符合程度和精细化程度可能会帮助我们提高对地壳运动和地震活动的了解。因此，本文纳入天山地区主要活动断裂，根据天山地区的地块构造、活动断裂带分布和流变结构等资料，建立了1种弹性模型和2种黏弹性模型，共3种不同的三维块体运动断层位错模型。在GPS地壳运动观测结果的约束下，附加断层滑动速率和闭锁深度约束，使用有限元程序Pylith(Aagaard，2017)模拟天山地区的地壳运动变形，讨论岩石圈纵向分层和黏弹性效应对模拟地表速度场结果的影响，计算获得研究区地壳速度场及应力场，并分析其量值大小、分布形态及特征。

1 三维有限元模型的建立

1.1 几何模型构建

考虑到天山地区的GPS测点分布及地震活动性(图1)，将研究区选取为(75°～94°E，38°～47°N)，包含了中国境内天山及境外天山的一部分，区域面积约为1 565 124 km。大陆板块内部以块体运动为特征(邓起东等，2002)。本文参考张培震等(2003)对西域活动地块区的划分，将研究区划分为塔里木、天山和准噶尔3个活动块体。由一系列断裂耦合成的南天山山前断裂带和北天山山前断裂带为研究区内的重要活动构造带，作为划分块体的边界带。

图1 天山地区GPS测点分布与M≥5.0地震活动特征

模型厚度取为80 km，将纵向上分为上地壳、中地壳、下地壳和岩石圈上地幔4层介质，以构建均匀分层模型和非均匀分层模型。均匀分层模型4层介质的厚度均为20 km，非均匀分层模型每层厚度依据全球地壳模型Crust 1.0(Laske，2012)，由每个块体内的不同单元的各层厚度取均值确定。由于纳入模型的断裂的倾角均较大(冯先岳，1986)，故将其作一定简化，固定为70°。另外，研究区的面积并不太大，可以将研究区域的表面作为平面处理。将关键点经纬度坐标投影转换成笛卡尔坐标后使用Coreform Trelis软件建立几何模型并进行网格划分，网格划分时断层面格网边长设为5 km，剩余部分按照距断层面的距离以1.07的渐变梯度进行划分。均匀分层模型包含985 052个四面体单元，177 174个节点；非均匀分层模型包含1 121 486个四面体单元，201 561个节点。图2展示了三维非均匀分层模型的几何图形及网格划分。

图2 三维非均匀分层模型的几何图形及格网划分

1.2 模型约束条件与物性参数

GPS观测得到的天山地区现今地壳运动结果，可为有限元模型提供可靠的运动学边界条件约束。Wang和Shen(2020)收集和处理了多个来源的数据，并整合以往研究中的速度解，得到了统一参考框架下的中国大陆及周边区域GPS观测结果，这些数据精度较高、密度较大，还给出了中国境外天山西部的GPS速度结果。因此，本文采用这些可靠且丰富的数据作为模型的边界条件。首先，选取靠近边界且不处于断层处的GPS测点，将其GPS速度值赋给离站点最近的模型地表侧边边界节点。然后，对已获得速度值的节点进行插值，并规定边界速度不随深度改变，进而获取剩余所有侧面节点的速度值(图3)。对于模型底面，规定其在垂向上的位移为零，水平向不作约束。最终，有限元模型上表面为自由表面，侧面在水平向上由GPS观测结果约束而在垂向上自由，下表面在垂向上固定而在水平向上自由。

图3 天山地区有限元模型边界条件

考虑到断层活动在地壳运动变形过程中起到至关重要的作用(郑勇等，2007)，除边界约束外，还需要对模型中的断层进行约束。Pylith软件通过自动创建“内聚单元”实现断层两侧顶点的相对运动，从而控制断层滑动，由此可以创建包含断层的不连续体模型。参考牛之俊等(2007)和刘代芹等(2016)给出的南、北天山断裂带的闭锁深度及滑动速率结果，经过网格搜索法实验，计算地表GPS站速度模拟值与观测值间的加权均方根误差，最终确定南天山断裂带的闭锁深度为17 km，左旋走滑速率约束为2.8 mm/a，85°E以西区域挤压速率约束为8.8 mm/a，85°E以东区域挤压速率约束为3.5 mm/a；北天山断裂带的闭锁深度确定为16 km，左旋走滑速率约束为2.1 mm/a，挤压速率约束为3.3 mm/a。

Maxwell流变包含与时间无关的弹性行为和与时间有关的黏性行为。Pylith软件由介质密度、横波速度和纵波速度定义弹性，由黏滞系数定义黏性。 Crust 1.0为分辨率1°×1°的全球地壳模型，将全球划分为64 800个单元，垂向上分为9层，每层均给出P波和S波速度以及密度参数。石耀霖和曹建玲(2008)对实验室流变实验结果应用于估算岩石圈等效黏滞系数中的多种影响因素进行了讨论，并以温度和应变速率的研究成果为基础，对中国大陆地壳和上地幔等效黏滞系数做出了较为详细的估计。本文参考这些结果，确定研究区介质的物性参数。考虑到中地壳及下地壳的黏滞系数随深度增加而减小，在确定介质分界处及底面深度处的黏滞系数后，规定其余深度黏滞系数由此进行线性插值得到，以尽量实现与真实情况的吻合。表1展示了非均匀分层模型的物性参数。

表1 非均匀分层模型的物性参数

至此，通过建立的几何模型、给定的约束条件和物性参数，经有限元数值计算后，即可得到模拟结果。模型通过指定运动学上的边界条件和断层滑移速率模拟应力绝对状态下的扰动，因此忽略了重力的影响。

2 模拟速度场结果及其特征

为研究岩石圈流变物质黏弹性效应对地表形变的影响，本文给均匀分层几何模型分别赋予弹性参数和黏弹性参数，建立均匀分层弹性模型(M1)和均匀分层黏弹性模型(M2)。两种模型唯一的区别为弹性模型将整个域视为弹性，黏弹性模型将上地壳视为弹性，中、下地壳和上地幔视为黏弹性。黏弹性模型与弹性模型对应的弹性参数相同。另外，为非均匀分层几何模型赋予黏弹性参数，建立非均匀分层黏弹性模型(M3)，以讨论纵向介质分层对地表形变的影响。最后，根据表1设置的黏滞系数，估计出黏弹性介质的最大弛豫时间。3种模型的数值计算采用500 a为1个时间步，经过100个时间步的计算，最终得到5万年后的地壳形变计算结果。

首先，由最后1个时间步的模型计算结果计算震间速度场。然后，从模拟速度场结果中提取模型所含地表中424个GPS站的速度值，与实际观测的GPS数据进行比较，结果见表2。在3种模型中，超过95%的站点速度模拟结果与观测结果差异在4 mm/a以内，速度方位模拟结果与观测结果差异在20°以内，且都基本服从正态分布。由此可见，3种模型均能较好地模拟出地表速度场，非均匀分层模型的模拟结果略优于均匀分层的模拟结果，说明介质厚度的改变对本文模拟地表速度场的影响并不大。此外，相较于纯弹性模型，黏弹性模型得到了更优的模拟结果。下面仅对拟合地表速度场效果最优的非均匀分层黏弹性模型(以下简称M3模型)的模拟结果展开分析。

表2 3种地表GPS站模型模拟速度与观测速度对比

从使用M3模型得到的GPS站点震间速度场与实测速度场的对比图(图4a)及速度残差图(图4b)可以看出，两者符合程度较好。图4c是使用M3模型给出的地表速度场。从大小来看，天山地区形变速率呈现由南向北、由西向东逐渐减小的态势。这是因为该区域从南向北、从西向东逐渐远离碰撞板块边界，板块的推挤效应逐步减弱。研究区运动速率最大位于模型西南角，即帕米尔—西昆仑弧形断裂以南区域，达20 mm/a。在模型所含区域内，南天山汇聚速率为12～15 mm/a，中天山汇聚速率为6～9 mm/a，北天山汇聚速率为4～5 mm/a。从方向来看，由研究区西南部的NNW向逐渐向东、向北过渡到NNE及NE向，从塔里木地区到准噶尔地区速度由顺时针方向转为逆时针方向，验证了塔里木地块的顺时针旋转运动。

图4 天山地区使用M3模型模拟速度场

3 模拟应力率场结果及其特征

现今构造应力场是现今地壳变形和构造运动的动力因素，是地震活动的内在直接原因之一(邓起东等，2000)。模拟得出的应力实际上是应力的平均年增量，即应力率，它反映了现今天山地区应力场的基本特征。由1900—2020年天山地区的地震目录可知，该地区大部分地震的震源深度在10～20 km。为研究天山地区地壳内部的应力状态，本文提取地表以下15 km的应力模拟结果，计算水平最大、最小主应力率。图5为天山地区地下15 km的水平方向主应力率图和最大剪应力率云图。

从图5a可以看出，研究区大部分区域的主压应力率大于主拉应力率。天山两侧褶皱带的主压应力率大于天山内部主压应力率，且呈现与山脉走向相正交的特征，与Chen等(2005)的研究结论相似。主压应力率较大的地区集中分布在天山褶皱带上，主压应力率最大的区域为南天山与帕米尔高原相接区，达4.780 kPa/a，此区域也是发生过几次强震的地震活跃区，其未来地震危险性值得重点关注。天山西段(75°～84°E)平均主压应力率为1.896 kPa/a，天山东段(85°～94°E)平均主压应力率为1.478 kPa/a。另外，塔里木地块和准噶尔地块作为较稳定块体，其应力率大致均匀分布且明显小于天山的应力率，平均主压应力率分别为0.471 kPa/a、0.235 kPa/a。图5b中的五角星代表发生于1900—2020年、震源深度介于10～20 km、≥5地震，由此可以看出，最大剪应力率大的地区发震次数多，地震活动性强，如南、北天山西段，最大剪应力率较大，地震频发；最大剪应力率小的地区地震少发，如南、北天山东段，塔里木和准噶尔，最大剪应力率小，发震的次数少。

图5 天山地区地下15 km截面模拟应力率场

研究区主压应力方向大多近N-S向，绝大部分最大、最小主应力接近水平。南天山从东到西主压应力总体方向为NNW向，主压应力轴集中在NNW15°～NNE10°，与李杰等(2012)根据GPS结果计算的南天山主压应变轴方向相符。北天山从东到西主压应力方向从NNE向过渡为NNW向，主压应力轴集中在NNW16°～NNE26°。

4 讨论

本文使用3种模型对地表速度场进行拟合，结果表明，使用M3模型计算得到的天山地区震间速度场与GPS观测结果有较好的一致性，说明本文建立的有限元模型是合理的，模拟结果是可靠的。此外，非均匀分层模型的模拟结果优于均匀分层模型模拟结果，并且黏弹性模型提高了对震间地壳形变的拟合程度，这与前人(王辉等，2007；Li，2015)得出的实验结果相一致。王辉等(2007)建立了川滇地区均匀分层和不均匀分层模型，探讨莫霍面起伏对模拟速度场结果的影响，结果表明不均匀分层模型的结果较均匀分层模型的结果有局部改善，但总体改进不多。Li等(2015)使用有限元模型来研究黏弹性对震间形变的控制作用，结果证实岩石圈流变物质的黏弹性效应对震间地表形变场有很大的影响。

具体来看，天山地区速度场从南向北、由西到东逐渐减小，运动速率最大达20 mm/a，位于模型西南角。在模型所含区域内，南天山汇聚速率为12～15 mm/a，中天山汇聚速率为6～9 mm/a，北天山汇聚速率为4～5 mm/a，这与牛之俊等(2007)利用GPS数据得出的天山地壳缩短速率相近。塔里木地块速度场方向顺时针旋转，准噶尔地块速度场方向逆时针旋转。本文模拟得到的天山地区现今构造应力场以近南北向水平挤压应力为主，与龙海英等(2007)利用震源机制解给出的解释一致。研究区主压应力率从南到北呈现不均匀分布特征，主要分布在天山南北两侧褶皱带，最大值处于南天山与帕米尔高原及西昆仑交界地带。主压应力率场成近N-S向分布。天山西段(75°～84°E)平均主压应力率较大。塔里木地块和准噶尔地块作为较稳定块体，其应力率大致均匀分布且明显小于天山的应力率。

本文模拟得到的研究区地壳形变速度场和应力分布结果显示，挤压速率在天山内部较为均匀，而在跨越断裂带时快速减小，表明断裂带起到了吸收形变的重要作用，最大剪应力结果也表明在这些区域积累了大量能量，这些现象与该区域的地震活动分布呈现出很好的一致性。因此，南、北天山山前断裂带的蠕滑对天山地区现今形变场和应力场起到一定的调整作用，天山地区地壳的缩短变形主要发生在天山南北两侧褶皱带，这影响了天山的地震活动性分布。

欧亚板块与印度板块之间的碰撞挤压是引起天山现今构造变形的直接原因。本文的模型受代表帕米尔高原推挤和塔里木顺时针旋转的运动学边界条件的约束，两者共同给天山地区提供了西强东弱的推挤作用力。本文模拟得到了与前人研究结果一致的由西到东的差异性构造运动特征。因此，天山地区的现今构造变形特征是帕米尔高原向北推挤和塔里木块体顺时针旋转运动共同作用的结果，但两者谁起到主导作用，还需进一步探讨。