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

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

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

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

为研究岩石圈流变物质黏弹性效应对地表形变的影响,本文给均匀分层几何模型分别赋予弹性参数和黏弹性参数,建立均匀分层弹性模型(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站模型模拟速度与观测速度对比
Tab.2 Comparison between simulation results and observation results of GPS stations

从使用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模型模拟速度场
Fig.4 Simulation results of velocity field by using M3 model

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

从图5a可以看出,研究区大部分区域的主压应力率大于主拉应力率。天山两侧褶皱带的主压应力率大于天山内部主压应力率,且呈现与山脉走向相正交的特征,与Chen等(2005)的研究结论相似。主压应力率较大的地区集中分布在天山

图5 天山地区地下15 km截面模拟应力率场
Fig.5 Results of simulated stress rate field with -15 km cross section

褶皱带上,主压应力率最大的区域为南天山与帕米尔高原相接区,达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、M≥5地震,由此可以看出,最大剪应力率大的地区发震次数多,地震活动性强,如南、北天山西段,最大剪应力率较大,地震频发; 最大剪应力率小的地区地震少发,如南、北天山东段,塔里木和准噶尔,最大剪应力率小,发震的次数少。

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

本文使用3种模型对地表速度场进行拟合,结果表明,使用M3模型计算得到的天山地区震间速度场与GPS观测结果有较好的一致性,说明本文建立的有限元模型是合理的,模拟结果是可靠的。此外,非均匀分层模型的模拟结果优于均匀分层模型模拟结果,并且黏弹性模型提高了对震间地壳形变的拟合程度,这与前人(王辉等,2007; Li et al,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)平均主压应力率较大。塔里木地块和准噶尔地块作为较稳定块体,其应力率大致均匀分布且明显小于天山的应力率。

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

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

本文通过建立3种不同的三维有限元模型,以GPS观测结果为边界约束,对断层附加滑移速率及闭锁深度进行约束,获取了天山地区的形变场及应力场,所得结论如下:

(1)模型纵向分层的不同对模拟地表速度场结果的影响并不大,而考虑黏弹性效应可以提高对震间地壳形变的拟合程度。

(2)主压应力率最大的地区是南天山与帕米尔高原及西昆仑交界地带,南、北天山西段的主压应力率也较大,这些地区未来的地震危险性值得重点关注。天山地区主压应力场方向近N-S向,呈现出以近南北向水平挤压应力为主的特征。

(3)南、北天山山前断裂带的蠕滑对天山地区现今形变场和应力场起到一定的调整作用。天山地区地壳的缩短变形主要发生在天山南北两侧褶皱带,这影响了天山地区的地震活动性分布。