基于断层倾角分段变化的玛多地震发震断层构造应力场演化数值模拟分析

吴立新， 卢菁琛， 毛文飞*， 胡俊， 周子龙， 李志伟， 齐源， 姚汝冰

1 中南大学地球科学与信息物理学院， 长沙 410083 2 中南大学地灾感知认知预知实验室， 长沙 410083 3 中南大学资源与安全工程学院， 长沙 410083

0 引言

2021年5月22日2时4分，青海省果洛藏族自治州玛多县发生了MS7.4强震，震中位于东昆仑断裂带以南约70 km(34.59°N，98.34°E)、鲜水河断裂的北侧，震源深度17 km(中国地震台网中心，2021).震后大量研究人员开展了野外地质调查研究(潘家伟等，2021；李智敏等，2021)、地球物理反演(王未来等，2021；尹欣欣等，2021；詹艳等，2021；杨君妍等，2021)、卫星遥感分析(祝爱玉等，2021；李志才等，2021；刘雷等, 2021; Qi et al., 2021；Wang et al.,2022).现场地质调查资料显示，玛多地震形成了包括张剪裂隙、挤压鼓包等多种破坏类型的复杂地表破裂带，为典型左行走滑位错性质.该破裂带主要沿东昆仑断裂带南部的江错断裂分布，整体呈N105°E走向，全长约151 km，最大水平位错量量约2.9 m，位于江错断裂西段末端(潘家伟等，2021).

玛多地震发生在青藏高原中部巴颜喀拉块体上，该区域是中国大陆的地震活跃区，其周缘孕育了包括东昆仑断裂、鲜水河断裂、甘孜—玉树断裂、马尔盖查卡—若拉岗日断裂、龙门山断裂和阿尔金断裂在内的多条活跃断裂(张培震等，2003；邓起东等，2002，2010).在印度板块的挤压作用下，巴颜喀拉块体呈NE向运动，且其滑动速率自西向东逐渐减小(Kirby et al.，2007；李陈侠等，2011；李煜航等，2015).近二十多年来，在巴颜喀拉地块周缘已经发生了9次7.0级以上的大地震，如2008年四川汶川MS8.0地震(徐锡伟等，2008b)、2010年青海玉树MS7.1地震(陈立春等，2010).

玛多地震的发震断层为江错断裂，是一条与达日断裂、玛多—甘德断裂和东昆仑断裂近平行展布的深大断裂(潘家伟等，2021，王未来等，2021)，总长约700 km，整体走向为NWW，断层倾向以NE为主，倾角近垂直，为左旋走滑性质(王未来等，2021).从活动断裂划分上来看，江错断裂属于东昆仑断裂南侧的次级断裂(邓起东等，2002；徐锡伟等，2008a)，玛多地震震中大体位于地表破裂带的中段，且在破裂带两端均出现分枝状破裂现象(潘家伟等，2021).

据已有研究，地震序列间大多具有触发和延缓效应(程佳等，2011；杨兴悦等，2013；李平恩等，2019；刘雷等，2021)；而本次玛多地震沿发震断裂带两端形成了分叉式地表破裂现象，这在青藏高原的走滑地震中首次显现(潘家伟等，2021)；而且，断层倾角的空间展布变化相较前序地震更为复杂，呈现出一种凹凸交替的分段变化形态(徐志国等，2021；王未来等，2021).因此，有必要详细分析玛多地震发震断层上构造应力的演化特征及未来破裂发震的危险性.

目前针对区域构造应力场的演化模拟研究中，陈连旺等(2008)模拟了在川滇地区主要活跃断裂上发生的强震对潜在强震库仑破裂应力的加卸载作用，发现活跃断裂同震时引起的破裂库仑应力变化大多呈增加趋势.杨兴悦等(2013)采用单元降刚法模拟了自1900年以来巴颜喀拉块体及周缘地区发生的MS7.0以上的地震，计算了该地区先前地震对后续地震的触发效应，结果表明：历史地震的诱发效应并不是完全连续的，而是具有一定的时间滞后性.吴萍萍等(2014)针对1893年以来强活跃性的左行走滑断裂——鲜水河断裂带发生的MS6.7以上地震，计算了各次地震对于后续地震的同震库仑破裂应力的变化，以及库仑应力的变化对后续地震的触发作用.李平恩等(2019)以GNSS水平运动速率观测值和最大主压应力方向测量值为约束条件，获得了构造背景的初始应力场；进而模拟自1900年以来巴颜喀拉块体周缘MS7.0级以上强震序列，研究了强震与库仑破裂应力变化的关系；结果表明：巴颜喀拉活跃块体上的强震对后续地震的发生既有触发效应也有延缓作用.而针对2021MS7.4玛多地震，祝爱玉等(2021)基于InSAR地表形变观测结果开展了区域应力场数值模拟计算，分析认为发震断层西南端和东北端应力场变化有利于未来余震的发生，且玛多地震的发生对部分临近断裂带上的应力分布具有显著影响，增大了其未来发生滑动的可能性.冯淦等(2021)基于地震破裂模型和均匀弹性半空间模型，计算和分析了玛多地震对其周围地区应力场和位移场的影响，发现玛多地震的发生，增大了发震断层西段末端以及东段末端西侧区域的库仑破裂应力.

鉴于模型复杂度会在很大程度上影响数值模拟计算效率以及结果的收敛性，且从尺度效应上来说适当简化模型并不会对模拟结果产生过大的影响，过去在研究块体边界或某一断裂带上发生的地震时，大多采用简化模型——将断层处理为单一倾角断层或直立断层.但对于2021年玛多地震而言，震后余震丛发，断层空间的分段与倾角变化十分复杂，若依旧采用简化模型进行模拟，则难以从应力演化的角度来解释其余震位置、地表破裂的分段性特征(潘家伟等，2021；王未来等，2021).因此，本文基于玛多地震发震断层面结构产状的分段特征(王未来等，2021)，精细化重构了发震断层面模型，以此对玛多地震发震断层面上构造应力场的时空演化过程进行了数值模拟分析.

1 三维走滑断层有限元模型

1.1 模型构建

为模拟玛多MS7.4左行走滑地震的应力演变过程，根据已有的地质资料以及震后中国地震局台网提供的信息，构建了33°51′N—35°30′N，96°41′E—100°33′E区域范围内的有限元模型.模型覆盖范围内的构造断裂包括：达日断裂、江错断裂、甘德南缘断裂、西藏大沟—昌麻河断裂，如图1所示.

震后基于余震序列的精确定位研究发现(王未来等，2021；徐志国等，2021)，江错断层(发震断层)并不是整体上单调南倾或北倾，其走向和倾角均呈现分段变化的特征，故该断层在空间上呈现出凹凸体形态.因此，在构建江错断层模型时，根据地震序列重定位所获得的各分段走向及倾角结果(王未来等，2021)，对江错断层的走向和倾角进行了细化处理.如图2a所示，江错断层的走向可分为13段，由西向东依次为AB段、BC段、…、MN段，各分段走向介于280°～298°之间，而断层整体平均走向为285°(王未来等，2021).相应地，发震断层各分段具有不同的倾向和倾角(图2b)，断层中间段区域(F→K)整体为北倾，倾向为11°～27°，倾角为76.6°～86.4°，东、西两段(A→F和K→N)则为南倾，倾向为190°～208°，倾角为76.6°～89.3°；为保证不同分段间的连续性，在构建几何模型时相邻分段连接处的断层倾角进行过渡式缓变处理.考虑到本文主要是针对发震断层(江错断层)的模拟分析，图2c所示区域内除发震断层外，其余各断层均简化为直立断层.基于上述考虑，构建如图2c所示有限元网格模型.模型中各断层几何参数如表1所示.模型中X轴为东西向，Y轴为南北向，Z轴为深度变化方向；模型WE长233 km，NS长100 km，深度20 km；经自适应网格划分后的节点数为135886个，单元数为37191个.采用有限元计算软件ANSYS对数值模型进行计算.

图2 三维有限元模型及断裂分布特征(a) 发震断层各分段平均走向，正北走向为0°； (b) 发震断层各分段平均倾向及倾角，如NE12°∠87.2°表示该段断层倾向为北东12°，倾角为87.2°； (c) 有限元模型及断裂分布(红色为发震断裂，黑色为非发震断裂，虚线框范围为本次精细化倾角区域). Fig.2 Three dimentional finite element model and faults distribution(a) Average strike of each fault segmentation, where north is 0°; (b) The average dip and dip angle of each fault segmentation,where NE12°∠87.2° illustrates that the fault segmentation dips to northeast (12°) with a dip angle 87.2°; (c) Finite element model and faults distribution (red represents seismogenic fault, black stands for non-seismogenic fault, and the dashed frame is sophisticated dip angle zone).

表1 模型几何及介质参数Table 1 Geometric and physical properties of model

1.2 本构关系及材料参数

从构造应力长期演变来看，地壳介质应当是黏性的，体现出在构造运动长期流变环境下的应力积累特性；从短期来看地壳介质又是弹性的.为了反映这种长期为黏性、短期为弹性的变形特征，本研究采用了黏弹性本构关系.根据研究区域的岩性分析结果(刘大明等，2018；白国典等，2018)，选取了模型主体介质参数的密度为2700 kg·m-3，弹性模量为1800 MPa，泊松比为0.28.同时参考已有围绕巴颜喀拉块体地震模拟中所采用的物性参数(李平恩等，2019)，选取介质的黏滞系数为1023Pa·s.考虑到断层介质的性质要弱于周围活动地块介质，故参照前人研究方法将断层考虑为具有一定厚度的软弱结构带，即主体介质的杨氏模量参数乘以弱化系数，以体现断层与周围块体间的差异，在本研究中所采用的弱化系数为1/3(杨兴悦等，2013；李平恩等，2019).因江错断裂南侧的断层活动程度远小于江错断层，故本研究将南侧甘德南缘断裂和达日断裂的强度进行适度提升、减弱其活跃性，以便进行相似性模拟.各断层的介质参数具体见表1.

1.3 边界条件及模型合理性验证

本文采用GNSS年均速度场(Zheng et al., 2017)作为有限元模型的边界约束，对GNSS观测值进行插值从而获得边界上的速度场.由于模型的实际地理范围为33°51′N—35°30′N，96°41′E—100°33′E，区内的GNSS观测站点数量有限(仅有9个)，若直接通过这些站点获得的运动速度值进行内插，会影响插值结果的可靠性，故将测站的范围适当放大至32°N—37°N，96°E—101°E(均在巴颜喀拉块体内部).通过对空间范围放大后的区域内测站数据(共35个测站)进行内插，获得原区域边界上的年均运动速度值，并与加载时间(10万年)相乘，从而获得模型边界上总的位移加载量.由于在深度方向上GNSS观测值的差异暂无法获知，因此本研究设定：在模型的不同深度下运动速度保持一致.此外，考虑到研究区域东边界受龙门山断裂带(近南北向)影响，其运动速率远小于其他三个边界，故将模型东边界视为固定边界.同时模型底部水平向自由，垂向固定，上表面为自由边界.

为了检验模型参数的合理性，在上述模型几何参数、物性参数及边界约束条件下进行加载模拟，对比分析模拟得到的速度场与GNSS观测获得的速度场.结果如图3所示，靠近模型四周边界上的模拟值和观测值在大小和方向上偏离程度均较小；而模型中部区域，二者在大小上虽然在局部区域存在一定差异(受断层阻挡导致其运动速率减小)，但从整体的运动方向上来看，模拟结果很好地体现了受印度板块推挤，巴颜喀拉块体内部区域呈自西向东运动的特征.

图3 GNSS观测值与数值模拟值的对比Fig.3 Comprison between GNSS observations and numerical simulations

2 模拟结果

2.1 库仑破裂应力变化速率

根据库仑破裂准则，当岩石趋近破裂时库仑破裂应力σCFS为

σCFS=τ+μ(σn-p)，

(1)

式中μ为摩擦系数，τ为剪切应力(沿断层面方向)，σn为正应力(压应力为正)，p为孔隙压力.地震模拟中通常采用库仑破裂应力的变化来表征应力的积累速率(Harris，1998)，当μ为常数时，可得库仑破裂应力的变化ΔσCFS为

ΔσCFS=Δτ+μ(Δσn-Δp)，

(2)

对于各向同性均匀介质，孔隙压力对摩擦系数的影响可用视摩擦系数μ′=μ(1-B)来表示，其中B为Skempton系数(取决于岩石膨胀系数和流体所占比例)，取值范围为0～1(Rice，1992)，则

ΔσCFS=Δτ+μ′Δσn.

(3)

视摩擦系数μ′取值通常为0～0.8，0～0.4为大变形，0.4～0.8为小变形，平均值约为0.4(Freed，2005).视摩擦系数取值大小对于库仑破裂应力的极性分布影响较小，仅对应力大小有一定的影响(万永革等，2000)，故在本研究采用平均视摩擦系数0.4进行计算.考虑到研究区构造运动的时间相似性(李平恩等，2019)，将获取的三维模型边界处运动速度场乘以10万年的时间尺度，从而获得模拟计算所需总位移加载量；为了便于分析应力演化过程，整个计算过程分为40个时间步进行位移加载.

由图4可知，江错断层面上库仑破裂应力年变化速率较大的位置有3处，分别为断层西段末端的中上部区域、震源附近和断层东段的末端区域(图4).A、B、C三处的年库仑破裂应力变化速率均为正值，即库仑破裂应力呈增加趋势，且各处最大库仑破裂应力积累速率分别为717.5 Pa·a-1、388.6 Pa·a-1、2691.2 Pa·a-1，分别是江错断层面上整体平均库仑应力积累速率(19.6 Pa·a-1)的36.6倍、19.8倍、137.3倍.

图4 江错断层面上库仑破裂应力的年变化速率Fig.4 Coulomb stress changing velocity per year on the Jiangcuo fault

从区域的差异性来看，震源以东断层面的平均库仑破裂应力积累速率为39.3 Pa·a-1，显著大于震源以西断层面的平均应力积累速率14.6 Pa·a-1.西段A区靠近鄂陵湖，深度范围为0～10 km；其中包含有库仑破裂应力变化速率为负值的子区域(图4中蓝色区域)，变化速率最大为-269.3 Pa·a-1，表明该负值子区域在构造运动作用下应力逐渐释放.中段震源附近(B区域)的库仑破裂应力积累速率呈现以震源为中心向四周辐射状减小的趋势，其中震源位置的库仑破裂应力积累速率为64.8 Pa·a-1.断层东段末端C区是整个江错断层上应力积累速率最大的区域，应力积累集中分布在0～20 km深度范围内，同时在8 km、12 km和20 km深度处存在3个积累速率大于1000 Pa·a-1的子区域；在4 km、10 km和7 km深度处，则有3个库仑破裂应力变化速率为负值的环状区域，其最值为-244.7 Pa·a-1.本研究在震源附近的模拟计算结果与李平恩等(2019)计算所得的库仑破裂应力积累速率(约为50 Pa·a-1)基本一致.

前人的地震模拟研究表明，库仑破裂应力积累速率快与断裂带历史强震复发周期短具有一致性，较低水平的应力积累速率对应强震发生的应力积累周期更长(刘雷等，2021).本次玛多地震所在江错断层面上A、B、C处的库仑破裂应力积累速率处于高水平状态，即这些区域更易形成应力集中区，意味着这些区域发生地震(或强震)的可能性要远高于其他区域，这或与玛多主震、余震的空间分布特征具有一定的相关性.

2.2 江错断层面上的应力演化

(4)

式中σx、σy、σz、τxy、τyz、τzx为应力张量分量.等效应力增大会加快介质失效而促进孕震过程，等效应力减小会阻碍介质屈服而延缓孕震过程.模拟得到的江错断层面上的应力演化结果如图5所示.

由图5可知，江错断层面上的应力演化过程为：首先在江错断层面的东段末端4～20 km深度范围内形成最值为4 MPa的应力集中区域(图5a)；随后，在江错断层面西段末端(鄂陵湖南侧区域)0～8 km深度范围出现了约为2 MPa的倾斜带状应力集中区(图5b)；之后，在震源附近出现了应力集中现象(图5c，呈向东、西两侧扩展的趋势)；此后，这3处的应力进一步积累、扩展和强化.当数值模拟的加载时间结束时，江错断层东段末端、江错断层西段末端和玛多地震主震震源处的最大等效应力分别为296 MPa、109 MPa和40 MPa.

图5 江错断层面上的应力演化结果(a) 加载时间为2500年； (b) 加载时间为7700年； (c) 加载时间为1.3万年； (d) 加载时间为2.1万年； (e) 加载时间为5.4万年； (f) 对应的加载时间为10万年(圆圈代表2021年5月22日至2021年11月1日内在研究区发生的震级≥3.0余震序列，其中圆圈颜色代表震源深度，圆圈大小对应不同震级，五角星为主震)； (g) 主震、余震震源位置与江错断层面的空间关系(水平投影).Fig.5 Tectonic stress evolution on the Jiangcuo fault(a) 2,500 years after loading; (b) 7700 years after loading; (c) 13000 years after loading; (d) 21,000 years after loading; (e) 54000 years after loading; (f) 100000 years after loading (circles are earthquake series which greater than or equal to MS3.0, from May 22, 2021 to November 1, 2021, in which colors stand for the depth of earthquake and sizes represent earthquake magnitude, and star means main earthquake); (g) Horizontal projection of epicenters of main shock, aftershocks and Jiangcuo fault.

玛多地震余震重定位的研究结果显示，玛多地震的地震序列呈现出明显的空间分段特征，在发震断层东段上存在着余震活动较少的稀疏带(王未来等，2021；尹欣欣等，2021)，这与图5g中江错断层面上东段的低应力区(K→L)较为相符.由此推测认为：断层面上局部区域应力积累的空间差异导致了余震的空间分段性.结合震级大小和地震数量来看，发震断层面上的地震序列基本可分为4段，其中有3处位于应力集中区，分别对应发震断层面的东、西两段末端及震源区域，有一处位于应力水平较低处(图5f).整体来看，地震序列分布与发震断层面上的应力积累情况具有较好的空间对应关系，积累应力较大位置对应的地震数量丛集、强度大.此外，发震断层面西段与震中之间的地震序列密集段对应的是应力水平较低处；考虑到模型的构建情况，东、西两段以及震源处的应力集中区均位于断层倾角变化的过渡处，在构造运动自西向东推挤作用下，易于在倾角变化处积累应力，从而诱发地震；而低应力处位于两分段的中间位置，并未出现明显的应力积累，这也说明，本研究的模型构建在分段划分上还可进一步细化.

2.3 区域构造应力演化

由图6可知，地表构造应力场演化过程为：整体构造应力场的背景值为0.1～1 MPa；首先在江错断裂西段末端地表破裂分叉处出现了较小的应力集中区(图6a)，大小约为1～2 MPa; 随后断层西段末端应力集中区向南北两侧扩展，且同时在江错断裂东段末端分叉处和西藏大沟—昌麻河断裂也出现了约为1～2 MPa的应力集中区(图6b).随着构造运动持续加载，东、西段末端区域的构造应力进一步积累并向外扩展，其中在西段分叉处构造应力增大至10～100 MPa(图6c)，同时，在鄂陵湖区西南侧也出现了局部的应力集中现象；此外，东、西段末端的应力集中现象，也开始沿着地表断层线向震中方向扩展.当10万年的加载过程结束时，断层东、西段末端的应力集中区都进一步扩展并相互连通(图6d)，整个断层上最大应力集中区为西端分叉处; 震中东、西两侧断层线上的构造应力约1～2 MPa，但震中位置却无显著的构造应力集中现象.

图6 地表构造应力场演化结果(a) 加载时间为3.1万年； (b) 加载时间为5.1万年； (c) 加载时间为8.2万年； (d) 加载时间为10万年.Fig.6 Evolution of tectonic stress on ground surface(a) 31000 years after loading; (b) 51000 years after loading; (c) 82000 years after loading; (d) 100000 years after loading.

作者在玛多地震后的当年6月份开展了为期一周的震区野外科考，沿发震断裂及环鄂陵湖追踪调查了地表破裂与灾情(图7)；结合潘家伟等(2021)报告的地表破裂情况分析，发现本研究揭示的地表构造应力集中区的分布情况与地表破裂带位置基本一致.从整体的分布形态上来看，地表破裂带在东、西两端均形成了大规模分叉的“帚状”结构，但震中位置地表破裂发育并不突出；地表破裂由震中向东、西两段逐渐发育并与东、西两段末端连续；与之对应，图5所示地表应力集中区分布情况也是在东西两段的末端出现分叉，在震中并无明显应力积累，且由震中位置向东、西两段逐渐积累.此外，从局部上来看，西段的地表破裂北支向北西消失于鄂陵湖，南支向西延伸至鄂陵湖南侧后逐渐消失，对应断裂带西段末端积累的应力，也是呈南、北两支消失于鄂陵湖；而断裂带东段的地表破裂北支破裂程度明显大于南支，这与东段末端的应力集中区分布差异也是一致的.

图7 玛多地震地表破裂现场调查结果(a—f分别对应图6d中的1—6处)(a) 沙土液化； (b) 间断式拉挤破裂； (c) 野马滩大桥垮塌； (d) 地表拉张开裂； (e) 黄河乡大桥桥墩破损； (f) 无人机探视.Fig.7 Surface ruptures investigation of Madoi earthquake (a—f correspond to 1—6 in Fig.6d respectively)(a) Liquefaction of sand layer; (b) Extrusion rupture on the ground surface; (c) Bridge on the Wild-horse wetland; (d) Tessive fissure on the ground surface; (e) Damaged bridge at the Yellow River town; (f) Unmanned aerial vehicle (UAV) visitation.

考虑到江错断层面的深度及其倾角分段性特征，计算震源深度处(17 km)水平面上的构造应力场演化，结果如图8所示.震源深度水平面内的构造应力集中现象主要出现在江错断层东段末端以及震源处.从应力演化过程来说，首先在东段末端与分支断层的交汇处出现了约1～2 MPa的应力集中现象(图8a).随后，东段末端的应力集中区域呈南北向扩展态势，其中心区域应力增大至10 MPa左右；同时，震源区域应力集中开始显现(图8b).当东段末端应力集中扩展至西藏大沟—昌麻河断裂和甘德南缘断裂时，震源区域应力集中沿断层线向东、西两侧扩展(图8c).当10万年加载过程结束时，东段末端的应力集中区域已显著影响到北侧的西藏大沟—昌麻河断裂和南侧的甘德南缘断裂，区域中心处最大应力增至100～300 MPa(图8d); 同时，整个发震断层都处于相对的高应力状态，震源位置应力也增至10～100 MPa，可能是玛多地震发震的根本原因.此外，从玛多震源深度平面内江错断层的走向分布来看，震源恰好位于两分段的倾角变化处，因此在构造运动的不断加载下，易于在该位置形成应力闭锁区；同时，整个巴颜喀拉块体在周缘块体的推挤下自西向东运动，这就导致震源处的应力沿断层面向东、向西扩展，由此形成了一个条带状的应力集中区.

图8 震源所在深度水平面内的构造应力场演化结果(虚线为地表断层线)(a) 加载时间为5000年； (b) 加载时间为1.6万年； (c) 加载时间为3.6万年； (d) 加载时间为10万年.Fig.8 Evolution of tectonic stress in the horizontal plane at hypocenter depth (dashed line represents fault line on ground surface)(a) 5000 years after loading; (b) 16000 years after loading; (c) 36000 years after loading; (d) 100000 years after loading.

对比地表和震源所在平面的构造应力时空演化过程(图6和图8)可知：首先，玛多地震主震震源所在水平面内应力集中现象主要体现在东段末端区域以及整个江错断层上，而地表的应力集中主要体现在东、西段末端区域(震中位置未见应力集中现象)，说明断层走向及倾角的变化对整个断层面上构造应力的积累与集中发育过程以及三维空间展布具有显著影响.

其次，在时间上，江错断层在玛多地震震中区域及断层东段末端区域应力演化要比同区域震源深度处的应力演化滞后，即当加载时间为3.6万年时，在震源所在水平面上，震源位置断层面和断层东段均出现了明显的应力集中现象(等效应力>6 MPa, 图8c)，而当加载时间为5.1万年时，地表才在江错断层东、西段的末端出现了一些应力集中现象(等效应力<6 MPa, 图6b).

2.4 讨论

前人针对巴颜喀拉块体及周缘地区的地震模拟中，主要构造断层或发震断层的产状大多数呈单一倾角或直立倾角(杨兴悦等，2013；李平恩等，2019；刘雷等，2021).玛多地震发生在巴颜喀拉块体内部的历史地震空白区，余震丛发且形成了分叉式地表破裂带，采用单一倾角进行模拟难以解释其孕震过程的特征性变化.为此，对图2b所示模型进行简化，进而对比分析断层几何模型的复杂度对计算结果的影响.

图9a所示模型不考虑断层走向及倾角的分段性特征，计算结果仅体现了构造运动自西向东不断推挤所导致的应力积累现象，断层面上积累的最大等效应力为60 Pa.图9b所示模型仅考虑断层走向分段而无倾角变化(均为直立)，计算结果虽然体现了一定的分段性特征，但与图5f所示计算结果相比，在深度方向上应力分布无显著差异，断层面上积累的最大等效应力仅为200 Pa.

图9 简化模型断层面上的应力积累(a) 发震断层直立且无走向分段； (b) 发震断层直立且有与图2a一致的走向分段.Fig.9 Accumulative stress on the fault of simplified model(a) The vertical seismogenic fault without strike sectioning; (b) The vertical seismogenic fault with strike sectioning as in Fig.2a.

显然，在以上两种模拟方式下，发震位置应出现在高等效应力的断层东段侧末端处，且不会在断层两段末端出现分叉式地表破裂带.因此，精细的断层模型方能更好地体现江错断裂东、西两段的应力分布差异以及玛多地震震源处的应力集中现象.综合对比图9与图5f的计算结果，从断层模型的应力积累大小及分布来看，对于玛多地震这种震后依旧持续活跃且具有特征性地表破坏的情况，采用精细化的断层模型进行有限元模拟分析是更为合适和科学合理的.

基于江错断层面上的库仑破裂应力变化(图4)及等效应力演化(图5)、地表及玛多地震主震震源深度水平面内构造应力场演化(图6和图8)、地表破裂发育(图7)、主震及余震序列空间分布(图6)及潘家伟等(2021)所做的地表破裂调查结果，综合分析认为：1)受江错断裂走向及倾角的分段变化影响，玛多地震震源区及断层东、西段末端为构造应力变化活跃区或高应力分布区，与玛多地震主震及余震序列的空间分布相对应；2)地表破裂带位置与地表构造应力集中区具有空间相关性，体现了玛多MS7.4主震会对地表高应力集中区域产生扰动，这也是导致地表分段式破裂的重要原因；3)顾及构造应力、地表破裂、余震活动的空间分布差异，可进一步将沿江错断裂的玛多地震孕震—发震影响划分为5段，即由西向东分别为A-E段、E-G段、G-K段、K-L段、L-N段(分段标识号见图2)：A-E段构造应力主要在10 km以浅局部积累和变化，地表可见形变及破裂(潘家伟等，2021)，且地表高应力区向断层南北两侧的扩展态势明显，余震在断层两侧均有发育，尤其向南侧延伸较远；E-G段构造应力在全深度上近似均匀分布和同步变化，地表形变与破裂发育显著，是余震丛发分布区；G-K段构造应力主要在深部(10～20 km)发育，应力闭锁区靠近主震震源位置，而5 km以浅及地表应力积累很弱，地表形变及破裂发育也不明显；K-L段与A-E段类似，应力积累也发生在10 km以浅，地表形变破裂显著(潘家伟等，2021)，但因其构造应力值小于A-E段，所以余震发育较少；L-N段应力积累作用明显、在深度上近似均匀分布，地表形变破裂及余震活动均显著发育，且在震源所在水平面上应力集中显著并向断层面南北两侧大范围扩展.

需要注意的是：在震源深度水平面内(图8)，震源周围的应力主要沿发震断层面积累和发展，并未向断层面南北两侧扩展，表明玛多地震主震及余震的发生已对震源区断层面上积累的应力有了较充分的释放(He et al., 2021)，体现了震源走滑发震的应力演化特征.而江错断裂东段积累的应力向南、北两侧分别扩展至甘德南缘断裂和西藏大沟—昌麻河断裂，尽管这两处断裂吸收了部分构造应力，但依旧越过这两个断裂继续朝南、北向扩展(图8d)，且其范围远远大于余震序列的空间分布范围(图5g)，表明在发震断层东段末端发生的余震并不能将江错断裂东段积累的应力全部释放.受限于本次模拟设计模型的时空范围，并未进一步获得模型东南侧的应力扩展情况.根据目前的应力扩展趋势以及巴颜喀拉块体内断裂的走向，我们推测：构造应力具有向江错断裂东段的南部地区继续积累的可能，巴颜喀拉块体东南区域可能是未来发震的高风险区.此外，考虑到地表构造应力在江错断层西段末端的显著积累与扩展(图6)，结合该段断层面上的应力积累变化(图4)，推测未来在江错断层西段附近区域可能会发生浅源(深度<10 km)小型余震.本研究虽然未考虑应力演化过程中玛多地震主震及余震对大区域构造应力的扰动，但就未来潜在发震区的空间分析和推测而言，与祝爱玉等(2021)、冯淦等(2021)、He等(2021)分析结果相符.

3 结论

2021年玛多MS7.4地震发生于巴颜喀拉块体内部的历史地震空白区, 震后余震丛发；地震所形成的地表破裂带在震中区并不显著，而在发震断层东、西段两端分叉式发育.本文基于地震序列精定位研究结果建立了断层面的三维精细模型(精细化发震断层的走向与倾角分段性变化)，对发震断层面及巴颜喀拉块体局部区域的构造应力场演化过程进行了数值模拟计算，揭示了模拟区断层及周边应力积累的时空过程，表明震源应力时空变化特征与走滑地震发震机制吻合、地表等效应力空间分布特征与地表破裂发育位置相符.本文研究揭示的现象及分析结论包括：

(1)在板块构造应力作用下，江错断裂面东、西两段末端区域和玛多地震主震震源附近区域的库仑破裂应力变化明显，且等效应力积累作用显著，与主震及余震序列具有较好的空间对应关系，表明江错断层面上的高应力积累是诱发玛多走滑地震主震及余震的主要原因.

(2)从地表等效应力演化的时空过程来看，地表沿江错断层应力积累的空间分布与地表破裂的发育位置基本一致，均为东、西两端较为显著且在末端呈南北向分叉发展的态势，而震中附近不甚发育，表明玛多地震发震过程对地表高应力积累区的扰动是导致地表形变破裂显著的主要原因.

(3)在震源深度处的水平面上，震源周围的应力主要沿发震断层面积累和发展，并未向断层面南北两侧扩展，体现了震源走滑发震的应力演化特征；江错断层东段的构造应力积累与扩展虽然受到了甘德南缘断裂和西藏大沟—昌麻河断裂的阻挡和部分吸收，但构造应力积累区的空间范围远大于玛多地震余震分布范围，表明本次玛多地震并没有将断层东段积累的构造应力全部释放.结合断层东段应力向南、北两侧差异扩展的趋势，推测巴颜喀拉块体东南地区可能是未来发震高风险区.

致谢震区野外科考工作中，得到青海大学史培军校长、刘峰贵教授的指导和帮助；关于余震发育及地表破裂的分段性，与国家自然灾害防治研究院的徐锡伟院长、申旭辉总工，中国地震局地质研究所单新建所长及其团队成员等进行了有益讨论；关于巴颜喀拉块体未来地震风险趋势，与史培军校长、徐锡伟院长进行了有益交流，在此一并表示感谢！同时感谢两位匿名审稿专家提出的有益意见和建议.