东昆仑断裂的震后粘弹性松弛效应对地壳变形的影响

付誉超 刁法启 朱亚戈 陈 飞 熊 熊

1 中国地质大学(武汉)地球物理与空间信息学院湖北省地球内部多尺度成像重点实验室，武汉市鲁磨路388号，430074

近30年来，以GPS和InSAR为代表的空间大地测量技术发展迅速，实现了高精度地壳变形监测，可为深入了解不同时空尺度下的地壳变形机理提供基础资料[1]。某一时段内观测到的地壳变形是区域地球动力学过程的综合反映，不仅包含稳态的构造变形，还可能包含强震引起的瞬态同震变形和几年至几十年时间尺度的震后变形等其他局部因素[2]。Calais等[3]研究发现，蒙古地区2个8级左右地震在震后约100 a仍能造成每年数mm的地壳变形。震后变形是震后过程研究中不可或缺的约束信息，但对稳态构造运动研究(如断层稳态滑动速率和块体运动学特征)而言，震后变形则被视为“噪声”，需在地壳变形研究中予以考虑[4]。

图1 研究区构造背景Fig.1 Tectonic settings of the study area

东昆仑断裂是横跨青藏高原的大型左旋走滑断裂之一，沿东西走向展布近2 000 km，以约10 mm/a的运动速率调节南侧巴颜喀拉块体和北侧柴达木盆地的差异[5]。自20世纪初以来，东昆仑断裂带附近共发生6次7级以上地震(图1)，分别为1937年花石峡7.5级地震、1963年都兰7.1级地震、1973年玛尼7.3级地震、1997年玛尼7.5级地震、2001年昆仑山口西8.1级地震和2021年玛多7.4级地震[6-10]。其中，发生时间较近的1997年玛尼7.5级地震、2001年昆仑山口西8.1级地震和2021年玛多7.4级地震均观测到明显的震后变形[11-14]。有学者分析震后变形的力学机制发现，粘弹性松弛效应是影响长时间尺度下大范围变形的主要机制，而震后余滑主要影响短期、近场震后变形[15]。上述研究有助于深入了解青藏高原东北缘下地壳/上地幔的流变结构和震后变形过程。然而，前人通常仅针对某一时期的单个震例进行研究，缺乏整条断层长时间尺度震后效应的综合分析。

基于此，本文以东昆仑断裂为例，结合前人对区域流变结构的认识，采用较为合理的介质模型和震源参数对该区域近百年来发生的6个强震的震后松弛效应进行数值模拟。在此基础上给出该效应对区域地壳变形的影响，并对结果进行分析和讨论，从而为其他大型断裂上强震的震后变形综合研究提供理论依据。

1 模型构建和参数设置

本文采用分层介质模型，结合地质、历史地震记录等资料，考虑较为真实的地震破裂分布情况，基于前人研究得到的区域下地壳/上地幔粘滞系数，综合分析东昆仑断裂带附近区域近百年来发生的6个M≥7地震的震后粘弹性松弛效应造成的地壳变形情况及时空演化过程。

本文采用PSGRN/PSCMP软件包[16]模拟东昆仑断裂附近区域强震的震后粘弹性松弛效应。分层介质模型包含弹性的上地壳、粘弹性的下地壳和上地幔。介质模型的弹性参数取自区域地震学研究结果(图2(a)、(b))，其中下地壳和上地幔的粘弹介质采用Burgers体流变结构表示(图3)。Burgers体由一个Maxwell体和一个Kelvin体串联而成，可描述介质的不同松弛时间特征，其本构关系可表示为[17]：

(1)

式中，τ1=η1/μ1，τ2=η2/μ2，τ12=η1/μ2。

在常应力加载状态下，应变随时间的变化可表示为：

(2)

式中，应变为瞬态的弹性响应、呈指数衰减的短期粘弹响应和线性变化的长期稳态粘弹响应等3个响应的线性叠加。

下地壳/上地幔的粘滞系数较难确定，但该区域的震后变形研究可为该参数的确定提供约束(图2(c))。Diao等[18]根据1999～2015年的GPS观测数据得出该地区下地壳的长期稳态粘滞系数约为1×1018Pa·s。Diao等[15]考虑震后余滑和震后粘弹性松弛效应的综合影响后，根据2001年昆仑山口西8.1级地震震后4个月的观测资料发现，下地壳和上地幔的粘滞系数约为1×1018Pa·s。Ryder等[19]以昆仑山口西8.1级地震震后1 a的GPS资料和震后2～5 a的InSAR资料为约束，得到其下地壳瞬态和稳态粘滞系数分别为1×1018Pa·s和1×1019Pa·s。根据上述研究结果，本文设置该地区下地壳瞬态和稳态粘滞系数的下边界分别为1×1017Pa·s和1×1018Pa·s，对应的上边界分别为2.5×1018Pa·s和2.5×1019Pa·s(图2(c))。根据文献[13]的研究结果，本文将该地区上地幔瞬态和稳态粘滞系数分别设置为1×1018Pa·s和1×1019Pa·s。

VP和VS分别为地震P波和S波波速，ρ为地球介质密度，不同数字序号为不同研究得到的结果，本研究中下地壳稳态粘滞系数η2取1×1018～2.5×1019 Pa·s(浅蓝色阴影部分)，橙色五角星为该范围的平均值图2 介质模型参数Fig.2 Parameters of the media model

η1和η2分别为瞬态和稳态粘滞系数，μ1和μ2分别为2个弹簧的剪切模量图3 Burgers体结构Fig.3 Burgers body

本文根据2001年昆仑山口西8.1级地震和1997年玛尼7.5级地震的同震破裂分布深度[20]以及该地区中小地震分布深度(图2(d))，综合确定弹性上地壳的厚度[21]。大地测量数据反演结果表明，上述2个地震的同震破裂主要分布在20 km以上，且95%区域的中小地震也分布在20 km以上。上述统计结果与Ryder等[19]选取的该地区弹性上地壳的厚度(16.5 km)较为接近，因此本文将弹性上地壳厚度设置为16.5 km。本文进一步根据区域地震学成像结果[22]和震后变形研究结果[19]，将下地壳厚度设置为65 km。

此外，同震破裂模型也会影响粘弹性松弛效应。本文采用由InSAR资料反演的精细位错模型[14,20,23]对1997年玛尼7.5级地震、2001年昆仑山口西8.1级地震和2021年玛多7.4级地震进行研究，而对早期发生的3个地震采用Xiong等[24]估算的结果进行研究。

2 结果与分析

由于该地区下地壳粘滞系数存在较大的不确定性，会对结果产生显著影响，因此本文将粘滞系数设定为图2(c)中取值范围的平均值(稳态粘滞系数设置为1.3×1019Pa·s)进行模拟分析。由图4(a)可知，由于1937年花石峡7.5级地震、1963年都兰7.1级地震和1973年玛尼7.3级地震的离逝时间相对较长[24]，因此这些地震的震后粘弹性松弛造成的地表变形相对较弱(<1 mm/a)。说明对于现今地表而言，该区域的震后变形主要由离逝时间较短的2001年昆仑山口西8.1级地震和1997年玛尼7.5级地震[24]共同构成。模拟结果显示，2010～2020年震后粘弹性松弛效应使得东昆仑断裂带两侧的地表显著变形，其中最大地表变形发生在昆仑山口西8.1级地震震源区南侧，其东向变形速率可达7.6 mm/a，而断裂带北侧的西向最大地表变形速率约为6.0 mm/a。此外，震后粘弹效应在1997年玛尼7.5级地震震源区附近引起的最大地表变形速率约为6.0 mm/a。

对青海玛多7.4级地震造成的震后粘弹效应进行模拟，由图4(b)可见，2021～2030年该地区预计可产生最大约3.9 mm/a的地表变形。该结果可为相关震后变形研究提供理论依据，同时表明在进行其他研究时，应考虑地表变形的影响。

为更直观地展示东昆仑断裂附近强震震后松弛效应的影响，将模拟的震后位移场投影到平行于断层走向的剖面上，并将其与震间变形模型预测的地表变形速率进行比较。图4(c)、(d)分别将速度场投影到A-A′和B-B′剖面上，其中剖面A-A′ 和B-B′处的断层滑动速率(10 mm/a)取自地质学和大地测量学研究结果[4,18]，断层闭锁深度(16.5 km)则取自Ryder等[19]的研究结果。

由图4可见，剖面A-A′上变形速率分量主要反映的是1997年玛尼7.5级地震和2001年昆仑山口西8.1级地震的震后变形结果。由图4(c)可见，2005～2010年和2010～2020年震后粘弹性松弛效应对剖面A-A′造成的最大影响分别约为27.2 mm/a和8.4 mm/a，均高于震间变形模型预测的变形速率。对玛多7.4级地震而言，2021～2025年和2025～2030年粘弹性松弛效应对剖面B-B′造成的最大影响分别约为3.4 mm/a和2.6 mm/a，其影响不可忽略。

模拟结果显示，距断层较近区域(≤100 km)的震后松弛效应的影响与震间变形模型预测的结果大致相同(图4(c)、4(d))，说明若仅根据观测到的地表速度场来研究断层的滑移率而不考虑其震后松弛效应，必然会使研究结果产生较大偏差。

3 讨 论

3.1 粘滞性参数对研究结果的影响

3.1.1 下地壳粘滞系数的影响

研究表明[26]，震后粘弹性松弛效应响应时间主要受粘滞系数约束。本文以2001年昆仑山口西8.1级地震为例，选取震源区附近一点(92.53°E, 36.33°N)，讨论不同下地壳粘滞系数的震后地表水平位移场时空变化特征(图5)。由图可见，下地壳粘滞系数变化对粘弹性松弛效应的弛豫时间影响显著。粘滞系数越小，早期地表变形速率越快，粘弹性松弛效应衰减也越快；粘滞系数越大，应力松弛越慢，地表变形速率越慢，震后粘弹性松弛效应衰减越慢，持续时间也更长。

基于此，为研究粘滞系数不确定性对模拟结果的影响，对前文选取的下地壳粘滞系数范围进行讨论(图2(c))。由图6(a)可见，2005～2010年震后松弛效应在断层附近区域造成的地壳变形速率为9.6～44.4 mm/a；由图6(b)可见，2010～2020年变形速率迅速衰减至mm量级，但仍可达3.6～9.0 mm/a；由图6(c)可见，2021～2025年东昆仑断裂附近系列强震的震后松弛效应造成的地表变形速率可达6.2～11.7 mm/a；由图6(d)可见，2025～2030年该效应引起的周边地表变形速率预计约为5.9～7.6 mm/a。上述分析表明，当下地壳粘滞系数在取值范围内变化时，以1997年玛尼7.5级地震、2001年昆仑山口西8.1级地震和2021年玛多7.4级地震等为代表的东昆仑断裂带附近强震的震后松弛效应，在震后数10 a内对震源区附近地表造成的影响仍可达每年数mm。

(a)和(b)表示前5次地震的影响；(c)和(d)表示所有6次地震的影响图6 下地壳粘滞系数对地表变形的影响Fig.6 Effect of lower crust viscosity on surface deformation

由图7(a)可见，2005～2010年前5个地震的震后粘弹性松弛效应对剖面A-A′造成的影响可达7.0 ～43.0 mm/a，即便是距断层超过100 km的区域，断层两侧的变形速率仍大于震间变形模型的模拟结果；由图7(b)可见，2010～2020年断层附近A-A′剖面变形速率可达3.4～9.4 mm/a；由图7(c)可见，对于发生时间较近的2021年玛多7.4级地震，2021～2025年震后松弛效应对剖面B-B′造成的影响约为2.3～9.5 mm/a；由图7(d)可见，2025～2030年震后松驰效应对剖面B-B′ 造成的影响预计会减小至1.7～6.2 mm/a，其影响不可忽略。

(a)和(b)表示前5次地震的影响；(c)和(d)表示玛多7.4级地震的影响；灰色阴影区域表示模拟得到的跨断层变形速率范围，黑色曲线和红色曲线分别表示模拟结果的平均值和根据震间变形模型得到的结果[25]图7 下地壳粘滞系数对跨断层滑动速率的影响Fig.7 Effect of lower crust viscosity on slip rate across the faults

3.1.2 下地壳介质非均质性、流变模型及上地幔粘滞系数选取

研究发现，东昆仑断裂带南北两侧下地壳粘滞系数存在明显差异[21,27]，说明建立更为精细的三维流变结构有限元模型可进一步深入理解震后变形过程。然而，目前的观测数据对该区域下地壳粘滞系数的约束能力较差，不同研究给出的该地区下地壳粘滞系数结果之间的差异(10～100倍)甚至远高于断裂带两侧粘滞系数之间的差异(约10倍)[28]。因此，本文采用横向均匀的介质模型进行震后粘弹性松弛模拟研究。

相比于Maxwell体，Burgers体需要考虑更多参数，并且能同时兼顾震后瞬态和长期的变形模拟[29-30]。本文在采用Burgers体进行模拟时，根据前人研究结果，将瞬态粘滞系数设定为稳态粘滞系数的1/10[13]。模拟结果表明，由于该区域的震后松弛主要由下地壳的粘性流动引起[19]，上地幔的贡献相对较小，且现有数据对该区域上地幔粘滞系数的约束能力较弱，研究相对较少，因此本文将上地幔粘滞系数设为固定值。

3.2 震后粘弹性松弛效应对断层滑动速率估算的影响

Bell等[31]分析发现，1997年玛尼7.5级地震前玛尼断裂周边的地表变形速率约为1 mm/a，在玛尼7.5级地震发生10 a后，玛尼断裂附近的地壳变形速率可达8 mm/a，远大于震前的变形速率。由图4(a)可见，2005～2010年震后松弛效应在玛尼断裂附近造成的地表变形速率最大约为6.0 mm/a，说明震后松弛效应或许可以解释1997年玛尼7.5级地震附近现今地壳变形速率明显大于震前变形速率这一现象。Li等[32]研究发现，2009～2015年昆仑山口西8.1级地震震源区域震后地表变形速率相较于震前增加约40～80 mm/a；Liu等[33]研究发现，2007～2015年昆仑山口西8.1级地震的震后变形速率约为10 mm/a。由图6(a)可见，昆仑山口西8.1级地震震后松弛效应在震后10 a内造成的地表平均变形速率约为9.6～44.4 mm/a，与前人研究结果基本一致。

Gan等[34]利用跨东昆仑断裂GPS观测数据和弹性半空间模型得到东昆仑断裂中段的滑动速率约为20 mm/a，而地质学研究结果[35]与其他大地测量研究结果[5]中的滑动速率分别约为8～12 mm/a和10～15 mm/a。Wang等[36]与Li等[32]将地质学和大地测量研究得到的滑移速率结果差异归因于不同的空间位置，Gan等[34]推断该差异可能是由震后粘弹性松弛效应所致。

基于此，本文根据Zheng等[37]发表的GPS速度场，采用2001年昆仑山口西8.1级地震近场GPS台站数据(图1中红色三角形)构建垂直东昆仑断裂中段的速度剖面，并使用一维模型搜索该段断层的滑动速率。已有研究表明，垂直于断层走向的高密度GPS观测数据可为断层滑移率估算提供可靠约束[38]。本文采用前人基于弹性位错理论提出的震间变形模型进行反演[25]，该模型中平行于断层走向的台站变形速率V(x)可用断层闭锁深度D、滑动速率V和台站距断层的距离x来表示：

(3)

图8 断层滑移速率反演结果对比Fig.8 Comparison of inverted results of fault slip rates

由式(3)可得到一条关于x的反正切函数曲线，即震间变形曲线(图8(c)中绿色和黑色虚线)。需要说明的是，本文采用的是Zheng等[37]发表的速度场数据，但并未考虑昆仑山口西8.1级地震震后粘弹性松弛效应的影响，因此分2种情况展开研究：1)情况1，直接采用Zheng等[37]发表的速度场进行断层滑移速率反演；2)情况2，首先在Zheng等[37]发表的速度场基础上扣除该时段内模拟的震后粘弹性松弛效应，以此获得震间稳态的地表速度场，然后进行断层滑动速率反演。初步测试发现，地表变形数据对断层闭锁深度D并不敏感，因此对断层闭锁深度D取固定值。由于该断层闭锁深度D与区域弹性上地壳层厚度[19]具有较好的一致性，因此在参考昆仑山口西8.1级地震同震破裂深度分布结果[20]时，将断层闭锁深度D固定为16.5 km。各台站原始GPS观测数据和扣除震后粘弹性松弛效应后的地表变形数据如图8(c)中蓝色和红色误差棒所示。由图8(a)、(b)可见，情况1得到的东昆仑断裂中段的滑动速率V约为13 mm/a，情况2得到的滑动速率V约为10 mm/a。需要说明的是，由于数据限制，本文未考虑东昆仑断裂南侧块体内部次级断层的作用，但这些次级断层的滑动可能也会对反演结果产生一定的影响[38]。

由图8可见，情况2反演得到的东昆仑断裂中段的滑移速率结果与地质学研究结果[35]具有更好的一致性，而情况1反演的结果则会高估约30%的断层滑移速率。Li等[32]利用2001年昆仑山口西8.1级地震震后GPS观测数据(2014～2017年)反演得到的东昆仑断裂中段的滑移速率(17.4±1.2 mm/a)远高于使用震前GPS数据(1991～2001年)反演结果(7.2±2.2 mm/a)；而Wang等[36]采用1999～2014年GPS数据，剔除该时段内强震的同震和震后变形影响后得到的该段断层滑移速率(10.0±0.3 mm/a)与本文结果十分吻合，说明在稳态断层滑移速率研究中，需要考虑并扣除强震震后粘弹性松弛效应。

由图4(c)可见，对于现今地表而言，震后松弛效应对昆仑山口西8.1级地震破裂带附近地壳变形的最大影响可达8.4 mm/a，与震间积累的变形速率影响相当。因此，在利用现今大地测量观测资料研究断层滑动速率时，震后松弛效应的长期影响不容忽视。

3.3 不同震后变形机制的影响

震后变形过程可能是震后断层余滑、下地壳/上地幔的粘弹性松弛和震后孔隙流体调整的共同作用结果[39]。Liu等[33]认为，2001年昆仑山口西8.1级地震震后短期变形机制应该包括余滑和粘弹性松弛效应；贺鹏超等[27]认为，忽略震后余滑的影响会低估介质的粘滞系数；Zhao等[28]研究表明，余滑效应会影响下地壳瞬态粘滞系数。由于震后断层余滑和孔隙流体调整可能仅作用在震后很短时间内，且影响范围较小，因此在考虑震后长期影响时，需要将震后粘弹性松弛效应视为主要因素。

4 结 语

本文以东昆仑断裂为例，构建粘弹变形模型，并设定合理的粘滞性结构和地震位错模型，模拟近百年来该断裂附近6个M≥7地震的震后粘弹性松弛效应，并分析该效应对稳态地壳变形的影响。结果表明：

1)东昆仑断裂附近历史强震的震后粘弹性松弛效应在2010～2020年对地壳造成的最大变形速率约为7.6 mm/a，而该效应造成的最大跨断层变形速率约为8.4 mm/a，大于震间模型估计的断层两侧地表变形速率。

2)2021年玛多7.4级地震震后粘弹效应预计将在2025～2030年造成最大跨断层变形速率约为2.6 mm/a，若忽略震后粘弹性松弛效应，将会高估震间断层滑动速率。

3)考虑震后粘弹性松弛效应得到的东昆仑断裂中段滑动速率与地质学研究结果具有更好的一致性，而不考虑震后粘弹性松弛效应则会高估30%的断层滑移率。

致谢:2021年玛多7.4级地震震源机制解和区域历史地震数据分别取自美国地质调查局USGS(https:∥earthquake.usgs.gov/earthquakes/map)和中国地震台网中心国家地震科学数据中心(http:∥data.earthquake.cn/index.html)；文中大部分图件由GMT软件绘制，在此一并表示感谢。