2022年门源 MS6.9 地震前断层活动及应力状态的数值模拟

李 媛 杨周胜 庞亚瑾 梁洪宝 刘 峡

1)中国地震局地质研究所,地震动力学国家重点实验室,北京 100029 2)中国地震局第一监测中心,天津 300180 3)云南省地震局,昆明 650000

0 引言

据中国地震台网测定,2022年1月8日1时45分,青海省海北藏族自治州门源县发生MS6.9 地震,震中位于(37.77°N,101.26°E),震源深度10km,地表破裂约22km。此次地震的震中距1986年和2016年2次门源MS6.4 地震震中均约33km。3次门源MS>6.0强震的发生时间、空间距离均较近,发震断层均为冷龙岭断裂西段。与前2次逆冲型地震(姜文亮等,2017)不同,本次地震几乎为纯左旋走滑型,与现有研究显示的发震断层的运动性质一致。开展发震断层活动和应力状态研究有助于促进对地震孕震机理的认识。

针对多次门源地震,一些学者在同震、震后形变特征及应力触发等方面进行了研究(刘洋等,2019; 李振洪等,2022; 朱琳等,2022)。库仑应力结果显示,1986年门源MS6.4 地震及1927年古浪MS8.0 地震对2016年门源MS6.4 地震产生了应力加载效应(雷东宁等,2018); 2016年门源地震对2022年门源地震的发生有一定促进作用(李振洪等,2022; 朱琳等,2022)。但多期GPS应变结果(朱爽等,2022)显示,2016年门源MS6.4 地震的影响主要集中在震中断层附近非常小的区域,对2022年门源MS6.9 地震影响较小。本文认为,该局部断层及附近地震如此频繁,应力释放远高于其他区域,除与地震触发有关外,应该还与发震断层本身的运动状态、长期应力累积趋势和区域动力学背景有着必然联系。

1986—2022年门源系列地震的发震断层——冷龙岭断裂位于祁连山构造区。该区地处急剧隆升的青藏高原活动块体向内陆扩展的前缘部位,受到青藏高原块体长期的NE向推挤作用(邓起东等,2002),加之受到刚性的塔里木和阿拉善地块阻挡,新构造活动和变形强烈,断裂构造十分发育。区内大型断裂主要包括托莱山断裂、冷龙岭断裂、昌马-俄博断裂、祁连山北缘断裂、武威-天祝-庄浪河断裂等,主要展布方向为NWW向和NW向,主要运动方式为走滑、逆冲或两者的过渡类型(雷东宁等,2018)。研究表明,该区域的新构造活动以水平运动为主、垂直隆升运动为辅,两者相差1个数量级; 水平运动和变形呈明显NE向挤压缩短、顺时针旋转和向SEE挤出的构造变形特征(Wangetal.,2020; 郝明,2012,王双绪等,2013)。

关于冷龙岭断裂带晚更新世以来的活动特征,早年的地质研究曾估计其滑动速率约为19mm/a(Gaudemeretal.,1995; Lasserreetal.,2002),而近年来的研究已将其量值限定在3.35～6.4mm/a的范围内(何文贵等,2000,2010; 郭鹏等,2017)。在大地测量方面,基于高精度、长时间尺度的GPS观测数据,多位学者(甘卫军等,2005; Meade,2007; 李煜航等,2015)利用负位错模型、三维弹性球面块体模型和块体-位错模型反演了祁连块体边界主要活动断裂的长期滑动速率、亏损滑动分布及闭锁耦合变化,发现金强河断裂、毛毛山断裂和冷龙岭断裂的闭锁程度高。一些学者基于三维有限元黏弹性模型,利用单期GPS数据,在青藏高原东北缘及海原-六盘山地区开展了断层应力场状态及与地震活动性相关的数值模拟研究(杨兴悦等,2016; 石富强等,2018; 孙云强等,2018; 蒋锋云等,2021),结果显示,各个断层的应力并不均匀,强震主要发生在高应力状态的断层上,当区域应力累积到较高水平时,可在短期内经历一系列地震或地震丛集。现今祁连山构造区的木里-江仓断裂、金强河-老虎山-毛毛山断裂的剪应力积累较大。庞亚瑾等(2017)利用多期GPS数据,基于三维弹性模型研究了九寨沟MS7.0 地震前不同时段的构造应力积累情况,结果显示,现今的地壳应力积累过程有利于区域内左旋走滑型地震的发生。前人的研究为了解该区域地震过程中断层的运动、应力状态及变化特征提供了宝贵的资料和经验。但是,在确定断层不同区段的滑动速率方面,地质研究数据仍有待补充。由于在以往的数值模型中简化了断层形态,且受到研究侧重点不同等因素的制约,目前针对冷龙岭断裂在不同时间尺度上的断层运动变化及其与地震关系的研究还较少。鉴于模型的精细度将很大程度影响计算效率和结果,考虑到门源地震的发震断层在走向展布上的特殊几何形态,如果能提高模型的精细程度,对提升地震孕震过程中断层运动和应力特征的研究认识将会有较大的帮助。

基于前人工作的基础,本文利用地震前1991—2015期和2017—2021期的GPS速度场,综合考虑各级地质单元、地壳的分层结构、介质的横纵向不均匀性、复杂的断层系统及断层形态等多种因素,建立了祁连山构造区精细的三维有限元黏弹性模型,计算分析了祁连山构造区在长期的构造运动环境下区域应力积累的基本格局,区域内断层的长期滑动速率、应力累积速率,以及这些量值在门源MS6.9 地震前约5a时段内的变化特征。同时,结合跨断层短水准、震源机制解等资料,探讨门源MS6.9 地震前近震区强震发生的动力学背景和孕震环境。

图1 祁连山断裂带构造背景图

1 有限元数值模型构建

1.1 区域及活动断裂带的三维模型

为充分利用研究区内的GPS观测数据约束模型,同时考虑块体运动差异对断层运动的影响和模型计算的边界效应等,将模型扩大为青藏高原东北部的整个区域。如图2a 所示,三维模型考虑了地壳的分层结构、介质的横纵向不均匀性因素,涵盖了多个活动地块,东至鄂尔多斯地块腹地,东南至四川盆地,西北至阿拉善地块腹地,中南部为祁连山地块。各地块边界以中国大陆活动地块划分方案(张培震等,2003)为依据。模型厚度为100km,分为上、中、下地壳和地幔层,地壳分层结构主要依据全球地壳模型Crust1.0(Laskeetal.,2013)。由于莫霍面起伏对模型结果具有显著影响,且Crust1.0数据的精度较低(1°×1°),模型的壳幔分层界面采用精度较高的中国海陆1︰500万莫霍面深度数据(郝天珧等,2014),地壳厚度自西南的祁连山地块向东北的阿拉善-鄂尔多斯地块逐步递减。

图2 三维有限元模型的边界条件和拟合残差结果

同时,模型考虑了复杂的断层系统,包含30余条大型断裂,其中祁连山断裂带主要包括党河南山断裂、阴凹槽断裂、托莱山断裂、冷龙岭断裂、昌马-俄博断裂、祁连山北缘断裂、榆木山断裂、龙首山断裂、金强河断裂、武威-天祝断裂、庄浪河断裂等。断层形态主要依据中国断层图(邓起东,2007)详细刻画。如表1所示,大部分断层的倾角为60°～80°或更为陡立,以走滑运动为主。参考前人的研究(杨兴悦等,2016; 庞亚瑾等,2017; 石富强等,2018),不考虑断层介质横向不均匀性的影响,同时为兼顾计算收敛问题,将断层统一简化为宽约3km的直立弱化带,埋深至上地壳底部深度。

表1 模型物性参数设置

模型采用ANSYS软件中具有 2 次位移模式的Solid187六面体单元进行网格划分。断层单元的边长尺寸为3km,其他单元控制在5km以内。模型共有379633个单元,71080个节点。网格划分后的三维模型见图2a,其中图框部分为本文重点讨论的祁连山构造区。规定模型的xyz坐标系为:x轴向E,y轴向N,z轴向上。

1.2 模型物性参数的选取

模型中,上地壳采用脆性结构,其他分层采用黏弹性体,介质参数包括杨氏模量、泊松比和黏滞系数。在假定介质弹性各向同性的前提下,基于波速结构结果(Stolketal.,2013)和Crust1.0(Laskeetal.,2013)岩石密度结果,依据式(1)(张东宁等,1989)计算不同地块的杨氏模量和泊松比。结果显示,各块体间的弹性性质在横向上的差距较小,与前人的计算结果相符(Xiaoetal.,2015; 杨兴悦等,2016; 祝爱玉等,2019)。各地块在上地壳、中地壳、下地壳及上地幔的杨氏模量、泊松比和密度设置如表1所示。参考前人对断层的设置方法,将断层的杨氏模量设置为较周围介质偏小约1个量级,即1×1010Pa,泊松比较周围介质高0.03(Xiaoetal.,2015; 杨兴悦等,2016; 祝爱玉等,2019)。关于介质黏滞系数,相关研究表明,模型中不同地块的岩石圈流变性质在横向和纵向上均有差异,且在中国西部低黏滞系数的下地壳普遍存在,下地壳的黏滞系数比岩石圈地幔一般低1～2个数量级,尤其是青藏高原存在着强度最低的下地壳(石耀霖等,2008; 陈连旺等,2011; 孙玉军等,2013)。本文采用祝爱玉等(2019)和孙玉军等(2013)给出的黏滞系数结果。

(1)

其中,ρ为岩石密度,VP、VS为横波速度及纵波速度,ϑ为泊松比,E为杨氏模量。

1.3 模型加载和计算

为探索不同时间尺度形变特征,本文采用门源MS6.9 地震前1991—2015期(Zhengetal.,2017)和2017—2021期的GPS速度场作为模型边界加载条件。其中,2017—2021期GPS水平运动速率场的数据处理和解算由最新的GAMIT/GLOBK软件(2)http:∥geoweb.mit.edu/gg/。完成。在该时间段内青藏高原东北缘发生了九寨沟MS7.0 地震,在速度场中去掉了该地震的同震影响。

选取位于模型边界附近的GPS测站,将其观测结果通过双次样条插值至模型侧面作为水平向运动约束。注意在GPS站点的选取过程中,需要剔除与周边测点运动显著不协调的测点,以降低模型边界附近的断层活动等局部运动信息对模拟运算的影响。参考前人研究(刘峡等,2010; 石富强等,2018; 万永魁等,2021),设定水平速度不随深度改变。另外,模型的上表面为自由面,下表面水平向自由、垂向固定。

本文采用静力学分析,岩石圈视为线弹性体和理想黏弹性体串联的三维黏弹性模型,其变形服从Maxwell体的本构关系(Liuetal.,2003)。Maxwell体在短时间内表现为弹性体,随着外力作用时间的增加,储存在弹性体内的势能逐渐消失在黏性体中,随后应力松弛。弛豫时间一般定义为黏滞系数与杨氏模量的比值。当流变时间远大于弛豫时间时,应力得到松弛并趋于稳定(Caoetal.,2009)。本文模型的弛豫时间范围为114～1655a。本文将加载时间设定为30000a。第1步,施加1991—2015期GPS速度场约束,开始采用的步长为100a,持续加载至29800a后将步长改为1a,继续加载195a,此步计算的目的是得到一个反映较长时间尺度、相对稳定的运动场和应力场; 第2步,在第1步的基础上继续施加2017—2021期GPS速度场约束,加载时间为5a,步长为1a,其计算目的是获取2022年门源MS6.9 地震发生之前约5a内研究区运动场和应力场的变化。在计算过程中可以观察到,当加载至约10000a时,模型地壳和地幔应力已趋于平稳。计算采用Newton Raphson迭代算法进行求解,将计算划分为一系列载荷增量步,在每个步长内自动划分出多个子步长; 采用残余力和位移的范数控制收敛,收敛精度为0.1%。需要说明的是,本文计算显示,长期稳态构造加载下的应力场平均累积速率约为1kPa/a。而该区域历史大地震对区域应力加卸载的影响量级在10～100kPa(万永革等,2007; 邓园浩等,2018; 汤大委等,2023),相比于万年尺度的构造应力场而言,量值非常小,故本文在加载中忽略地震事件的影响。

2 数值模拟结果分析

2.1 数值模拟计算的合理性分析

图2b为1991—2015、2017—20212期GPS观测值和模型表面边界的加载。如图2c 所示,模型范围内198个GPS测点的观测值与模拟值的残差方向较为散乱,呈无序状态,且残差值大小不一,说明其中不包含显著的整体运动信息。对模型内GPS测点的N向和E向速率残差值进行统计分析,结果显示(图2d,e),2期残差值基本服从正态分布,残差<2mm的测站分别约占总测站数的90%和81%,其中残差<1mm的测站分别占比约76%和68%。综上,本文认为2期的边界约束和模拟结果比较合理。

此外,表2 给出了本文第1载荷步计算获取的各断层的运动性质与速率的模拟结果,作为对比,也列出了前人的地质考察结果。模拟得到的断层性质与前人的地质研究结果相符,模拟给出的大部分断层的滑动速率与地质结果相近,进一步印证了模拟计算的合理性。

表2 祁连山构造区主要活动断裂几何形态和活动特征表

2.2 1991—2015期区域应力场模拟结果揭示强震在长时间尺度上的孕育背景

应力场可直观地反映地壳受力状态及地震孕育的背景。图3a 为利用较长时间跨度的1991—2015期GPS速度约束,计算给出的研究区上地壳10～15km范围的最大、最小主应力和最大剪应力累积率结果,为应力的平均年变化率,代表了构造应力场长期的变化趋势。根据该图,应力场自西向东逐渐由NNE向挤压和NWW向拉张,顺时针旋转为NE向挤压和NW向拉张,主压应力方向大致与断层走向垂直,与中国地壳应力场(谢富仁等,2004)中的最大主应力方向相符。此外,本文收集整理了祁连山地区1970年1月—2021年9月利用HASH和CAP方法进行解算和精定位的,共619个3级以上地震震源机制解结果,基于断层面解的阻尼应力张量反演法(Hardebecketal.,2006),获得水平构造应力场,见图3b。本文的应力模拟结果为上地壳网格单元的应力结果,深度范围为10～15km,而震源机制解涉及的深度则更宽泛,尽管2种应力结果在某些局部存在差别,但总体上大致符合。值得关注的是,模拟结果显示在门源MS6.9 震源区及附近,应力场处于NE向挤压和NW向剪切共同作用的动力学背景之下,为主压应力场由NNE向向NWW向转换的过渡地区,与GPS应变率场研究结果(Zhengetal.,2013; Wangetal.,2020)一致。可见,现今应力场的长期的、整体的积累趋势,促进了NWW向断层的走滑和挤压运动,有助于托莱山断裂、冷龙岭断裂等左旋走滑断层上走滑型和逆冲型地震的孕育。

图3 研究区上地壳主应力累积率场分布

图3a显示,研究区自NW向SE最大剪切应力累积率逐渐减小,西部地区最大达2.4kPa/a,东南地区为0.6～1.2kPa/a。GPS研究给出的1999—2017年最大剪应变率结果(朱爽等,2022)同样显示了西强东弱的特征,与本文剪应力累积率相符合。此外,在托莱山和冷龙岭断裂交会部位,即门源MS6.9 震源区附近的剪应力累积率较周围大,最大为2.0kPa/a,其东西两侧约为1kPa/a,揭示门源MS6.9 震源区受到强烈的剪切作用,更容易发生地震。

2.3 1991—2015期断层长期平均运动与应力场模拟结果

本文模型利用连续、具有一定宽度、黏滞系数和杨氏模量都显著低于周围稳定地块的薄层单元,模拟断层的运动学特征和释放能量的力学性质,与实际断层的非连续错动滑移仍有区别。尽管如此,从模拟结果看,断层两侧的位移仍然有显著差异(图4)。参考已有的处理方法(祝爱玉等,2015),本文的断层挤压和走滑速率分量分别是断层两侧网格节点垂直于断层走向和沿断层走向的相对位移速率,相应地,断层应力累积率为垂直于断层走向和沿断层走向的法向正(挤压)应力和剪切应力累积率。

图4 断层长期滑动分量与应力积累速率模拟结果

图4a、b为利用较长时间跨度的1991—2015期GPS速度场作为约束,计算给出的断层长期平均滑动速率。结果显示,区内断裂活动总体表现为:大部分NWW向断裂受到左旋走滑兼挤压作用,走滑速率为0～5mm/a,挤压速率为0～3mm/a; NNW向断裂表现为右旋走滑兼挤压,如庄浪河断裂和榆木山断裂东段,与地质考察结果较为符合(国家地震局地质研究所等,1993; 侯康明等,1999; 袁道阳等,2002)。断层应力累积率结果(图4c,d)显示,力的性质与运动性质相匹配,大部分NWW向断层既受到挤压正应力又受到左旋剪切应力,NNW向断层受右旋剪应力。上述结果与杨兴悦等(2016)、戴华光等(1996)给出的有限元结果及李煜航等(2015)根据GPS计算的主要大断层的滑动速率一致,与该区域晚新生代以来处于长期挤压与剪切共同作用的动力学背景相符。

空间上,断层运动和应力呈现不均匀分布。研究区内断层的运动趋势为:自西向东挤压速率逐渐减小,由3.0mm/a减小为0.2mm/a; 自西向东走滑速率逐渐增大,由0.9mm/a增加为4.6mm/a。这揭示了祁连山构造带在青藏高原东北缘构造变形中的转换和调整作用,且尤其以托莱山断裂—冷龙岭断裂的断层运动和应力分段性特征最为明显:从托莱山断裂西段至冷龙岭断裂东段,挤压速率从2.5mm/a逐渐减少为0.2mm/a,而走滑速率从1mm/a逐渐增大至4.6mm/a; 剪切应力积累年变化率为0.5～3.8kPa/a,挤压正应力积累年变化率为0.2～1.0kPa/a。从图4b 中可以看出,相对周边断层,托莱山断裂和冷龙岭断裂左旋走滑运动较强,但在冷龙岭断裂西段的几何拐折部位,挤压运动分量和走滑分量较两侧断层段落更弱,而断层的法向挤压应力和左旋剪应力积累率较大,表明运动受阻。

2.4 2017—2021年模拟断层运动与应力增量场的结果分析

在2017—2021时段GPS速度场的约束下,计算得到了2022年门源MS6.9 地震发生之前约5a内的断层滑动速率和应力累积率分布。为了更加直观地反映断层运动的变化情况,将2017—2021期的模拟结果与1991—2015期结果相减,得到断层滑动速率和应力累积率的变化量,即增量(图5)。整体来看,挤压、走滑速率增量的绝对值范围分别为0～1.5mm/a和0～2.5mm/a。大部分断层段的变化量并不大,在1mm范围内,局部断层段变化较大,如榆木山断裂东段、党河南山断裂中段和冷龙岭断裂西段。断层的增量分布表现为显著的分段差异和不均匀性:一方面,走滑速率增量(图5b)显示,以101°E为界,西部大多数断层的走滑速率增量为正值,东部大部分断层为负值,表明在2022年门源MS6.9 地震前,西部断层左旋走滑运动有所增强,而东部则有所减弱,尤其是冷龙岭断裂西段减弱较为显著,量值约为1mm/a,而此段落的剪应力积累速率增量(图5d)为正值,即断层的剪切应力积累加速; 另一方面,图5a,c显示,挤压速率和正应力累积率(挤压速率和挤压应力为负值)增量为正的断层分段主要包括党河南山断裂、阴凹槽断裂、托莱山断裂东段、冷龙岭断裂、祁连山北缘断裂中段及NNW走向的断层,表明近期这些断层的挤压运动有所减弱,压应力积累减缓。尤其在托莱山断层东段及冷龙岭断层西段邻接的部位,近期断层的滑动速率和压应力累积率显著减弱,这一变化有利于断层错动。

图5 震前约5a较长期的断层滑动速率与应力累积率变化结果

3 讨论

3.1 发震断层长趋势运动和应力累积特点与门源 MS6.9 地震的关系

根据本文利用1991—2015期GPS结果给出的区域应力场、断层长期的运动和应力累积分布,整体而言,托莱山和冷龙岭断裂受剪切作用力较强(图3),断层左旋走滑速率达4.6mm/a,与野外地质考察显示的这2条走滑型断裂具有4.0～5.38mm/a和3.36～4.62mm/a的水平滑动速率(何文贵等,2000,2010)差别不大。但在托莱山断裂和冷龙岭断裂邻接几何拐折部位及附近区域,挤压速率和走滑速率比相邻的断层分段更弱(图4a,b),同时该部位的断层法向挤压应力累积率和左旋剪应力累积率较大(图4c,d),这与该断层段可能处于闭锁状态(叶茂盛等,2018; Guoetal.,2022; 刘雷等,2023)有关。在闭锁阶段,断层运动速率很小而能量快速积累,直到发生一次或若干次较快滑动,将积累的能量迅速释放掉。可以推断,这样的断层区段更容易发生地震。

中国大陆板内地震是在统一的、大的构造应力场作用下发生的(Tapponnieretal.,1977),主要受发震断层段的构造几何形态、介质物性、断裂活动习性和应力积累程度等孕震环境控制(陶玮等,2011; 张培震等,2013)。鉴于祁连山地块内部地壳结构和介质物性参数横向差异不大(Laskeetal.,2013; 郝天珧等,2014),推断引起冷龙岭断裂西段滑动和应力积累减弱的因素主要与该断层分段走向所呈现的特殊形态紧密相关。过去的研究多关注块体边界断裂,且大多对其形态进行简化。而1986—2022年门源系列地震发生在冷龙岭断裂西段走向呈现拐折形态的区域,具有特殊性,采用简化的断层模型难以解释其孕震过程的特征性变化。本文构建的连续三维黏弹性有限元模型(图2a 所示)考虑了各级地质单元、壳幔分层、地壳的横纵向不均匀性和介质的黏弹性等多种因素,同时构建了复杂的断层系统,细化了断层的几何形态,所获得的区域形变场、应力场和断层运动,应是周边块体相对运动、岩石层的黏弹结构和断层几何形态和空间展布等共同作用的结果。与已有模型的结果相比(杨兴悦等,2016; 石富强等,2018; 孙云强等,2018),本文模拟的断层运动及其应力分布具有显著的分段性特征,体现了断层的几何形态对断层活动的重要影响。

除2022年门源MS6.9 地震外,许多强震发生在断层拐折附近(Andrews,1994)。例如,San Andreas大断裂分别于1690年、1630年和1857年发生了3次MS>7.0地震,震中均位于断层拐折附近; 1976年唐山MS7.8、2013年芦山MS7.0 及1927年古浪MS8.0 地震同样为发生在拐折断层上的事件(Schwartzetal.,1990; Aochietal.,2002; 张竹琪等,2013)。Andrews(1994)指出断层的拐折部位对断层的滑动具有抑制作用; 薛霆虓等(2009)研究表明,几何弯曲的断层导致了应力的集中,断层拐折附近剪应力显著增加且多发育强震。针对拐折断层强震多发的位置,通过物理实验和数值模拟研究表明:地震并未发生在拐折显著的位置,而是发生在拐折附近或拐折两侧的平直断层处,且拐点起到类似阀门的作用(马瑾等,2012)。而1986—2022年3次门源MS>6.0地震的震中正位于拐折断层的两端靠近平直断层分段的位置。冷龙岭断裂带西段特殊的形态,很大程度上导致了该局部断层长趋势的运动亏损,与两侧的断层段运动存在显著不协调,促使应力在该断层段快速积累。因此,冷龙岭断裂带西段具有强震多发的孕震背景和构造条件。

3.2 发震断层震前中短期运动和受力特点与2022年门源 MS6.9 地震的关系

图6 为2022年门源MS6.9 地震前,祁连山构造区跨断层短水准观测的时序曲线,所用数据经过去周期和4～6阶多项式趋势拟合等处理。曲线趋势向下代表所跨断层出现挤压、逆冲运动,向上则为拉张、正断运动。将跨断层短水准表现出的张压运动特征,与2017—2021期相对于1991—2015期挤压速率增量的模拟结果(图5a)进行对比,发现两者具有较为一致的变化趋势。例如,祁连山北缘断裂东段上的扁都口、石灰窑口和冷龙岭断裂东段的白土庄跨断层场地短水准测段在2017年后均表现为挤压运动偏离线性和趋势减缓的特征; 大野口测段2018年以来为拉张增强或挤压减弱的运动趋势,与此同时,模拟结果(图5a)显示这些跨断层场地所处的断层分段,其挤压速率增量均为正值,表现为与跨断层水准曲线相符的挤压减弱变化。再如,金强河断裂上的乌鞘岭、松山跨断层短水准在2018年以来呈挤压增强的运动趋势,模拟结果(图5a)也表明该断层的挤压速率增量为负,近期挤压运动显著增加。又如,昌马-俄博断裂东段的童子坝和龙首山断裂东段的滑石口井短水准自2009年以来挤压运动速率几乎未发生改变,而模拟结果(图5a)显示相应断层分段的挤压速率增量接近零值。上述跨断层形变实测数据可与模型计算结果相互印证。

图6 震中距150km以内的异常测段时序曲线

2017—2021期相对于1991—2015期的断层运动速率增量结果(图5a,b)显示,2022年门源MS6.9 震前中短期,近震区断层的挤压运动和走滑运动有所减缓,尤其是托莱山断裂东段和冷龙岭断裂带西段,其走滑速率由原来的2.0mm/a减少至1.0mm/a,挤压速率由原来较弱的运动速率1.5mm/a减少至0.3mm/a,说明断层挤压和剪切变形在震前中短期进一步受阻,可能与地震前冷龙岭断裂存在闭锁现象(Guoetal.,2022)有关。利用GPS剖面分析得到的断层滑动速率和区域应变场动态结果(朱爽等,2022)也显示,冷龙岭断裂带的左旋走滑速率和区域剪切应变率在2017年后有所减弱。震前中短期运动趋势有减缓的现象,在1976年唐山MS7.8、1989年美国LomaprietaMS7.1、2013年岷县漳县MS6.6、2008年汶川MS8.0 等地震发生前也曾出现(Lindeetal.,1992; 杜方等,2009; 牛安福等,2020)。此外,2017—2021期相对于1991—2015期断层剪应力累积率的增量结果(图5d)显示,冷龙岭断裂西段的剪应力积累速率进一步增强,挤压正应力积累速率有所减缓,剪应力、正应力的变化有利于左旋走滑型地震的发生。上述现象说明震源区在震前中长期已处于高应力积累水平,随着地震临近,发震断层面法向挤压应力缓慢增高,剪切摩擦力同步增大,导致断层运动“迟滞”(李媛等,2023),促使运动趋势偏离线性和减缓。

不仅如此,2017—2021期相对于1991—2015期的增量结果也显示,在与冷龙岭断裂西段弯曲构造相似的六盘山东麓断裂的弯曲部位,断层挤压和走滑运动剪切变形在近期略有减少,运动进一步受阻,可能预示着六盘山东麓断裂的弯曲部位附近的强震危险性有所增加。

4 结论

本文利用1991—2015期和2017—2021期GPS速度场,通过建立青藏高原东北部三维有限元模型,计算了地壳形变和应力积累的长趋势特征及其在门源震前中短期的变化。结合震源机制解、跨断层短水准等结果,综合分析了门源MS6.9 地震近震区强震多发的孕震背景和震前形变特征,对其机理进行讨论。得到如下结论:

(1)在青藏高原NE向的推挤作用下,现今门源地震近震区受到NE-SW向挤压和NW-SE向拉张的构造应力场作用,有助于促进该区域内NWW向断层的左旋走滑和挤压运动,利于左旋走滑断层上地震的孕育。

(2)冷龙岭断裂西段具有强震多发的孕震背景和构造条件:1)门源MS6.9 地震震中局部区域最大剪切应力累积率比周围区域更大; 2)冷龙岭断裂西段呈几何拐折形态,其特殊的构造导致不同的断层分段之间应力积累的差异性和运动不协调,与两侧断层分段相比,该几何拐折部位的应力积累速率较高,但走滑运动在空间上受到了一定阻碍,断层速率偏低,运动存在亏损; 3)自2017年以来至2022年门源MS6.9 震前,冷龙岭西段断层剪切应力积累速率加快,断层运动尤其是走滑运动显著减弱,断层运动进一步“迟滞”。上述断层运动和应力状态有利于冷龙岭断裂西段上地震的孕育和发生。综上所述,冷龙岭断裂西段具有高应力积累状态和断层运动亏损的孕震背景,判断该断裂分段未来依旧是强震频发的危险区段。

此外需要指出的是,受计算量、收敛问题等多方面的限制,本文模型仅是在几何形态(即走向)上对断层进行了详细处理,不考虑倾角,将断层简化为直立断层; 同时,模型边界约束沿壳幔深度方向并无变化。这可能导致断层面倾滑运动分量的模拟值偏小,因此仍有待进一步改进。

致谢跨断层数据来源于中国地震局第二监测中心和甘肃省地震局; 本研究使用的有限元软件为中国地震局地质研究所提供的正版Ansys商用有限元软件; 马胜利研究员对文章给予了多次指导; 审稿专家为本文提出了宝贵建议,使得文章更加完善。在此一并表示感谢!