0.   引言

地震是在断层面剪切应力驱动下的一种自发破裂过程。从力学角度理解地震破裂，需要了解地震发生前断层的应力状态以及导致断层不稳定的条件。在过去的几十年里，大量学者针对后者做了许多相关研究工作，调查了断层不稳定的因素（Brace and Byerlee，1966），并揭示出各种断层弱化机制，如速率和状态摩擦（Cameron et al.， 1979）、热增压效应（Lachenbruch，1980）。然而，目前探索震间期断层剪切应力积累和分布的研究仍相对较少。地震动力学模拟和大地测量观测研究证明：认识断层剪切应力对于理解地震动力学过程十分重要。近年来基于断层粗糙度和初始应力不均匀性的地震动态破裂数值模拟研究（Ripperger et al.， 2007； Dunham et al., 2011）发现地震破裂的最终滑移分布在很大程度上受断层初始应力分布的影响。俯冲带的大地测量研究揭示出断层锁定程度的不均匀分布，表明沿俯冲带断层滑动亏损的不均匀（Chlieh et al., 2008），并产生了非均匀的断层剪应力分布。Yao and Yang（2022）采用震间闭锁分布估计断层三维应力，并将其作为初始应力开展地震动态破裂模拟，发现数值模拟预测的安宁河断裂破裂分段与历史地震完全一致，进一步证明了断层三维应力对于评估地震灾害的重要性。

基于空间对地观测的断层闭锁程度成像已经成为潜在震源区识别的重要手段。断层闭锁程度系数代表了震间期断裂的耦合程度，从运动学角度刻画了弹性应变积累的空间样式，被广泛运用于推测未来强震的位置和震级。例如，2010年智利Maule M_W 8.8地震同震滑移分布与自1835年M 8.5地震以来的滑动亏损分布具有密切的空间相关性（Moreno et al.， 2012）；2011年3月11日，日本Tohoku-oki大地震同震破裂的滑动分布与地震前俯冲带上的震间闭锁分布具有很好的空间对应关系（Loveless and Meade， 2011）；赵静等（2012）采用Defnode负位错反演程序估算了龙门山断裂在汶川地震前的闭锁程度和滑动亏损分布，其结果较好地解释了汶川地震时从震中向北东向的单向破裂现象。目前从力学角度仍然缺乏有关断层闭锁与应力积累相互关系的了解。力学上的断层闭锁应该由“二元状态”（Herman and Govers，2020）所描述：闭锁或不闭锁；然而大地测量给出的断层闭锁程度系数是一个介于0～1的无量纲数。如何理解二者的差异和联系？当采用断层闭锁程度来识别潜在的地震危险区时，强弱闭锁的界限如何取值？即闭锁程度高于哪个阈值时可以认为是闭锁较为显著？已有研究建议采用闭锁系数0.5作为强弱闭锁的分界（Li et al.， 2023），但这缺乏必要的力学依据。

岩石圈流变特性广泛应用于震后形变解释（Thatcher and Rundle， 1984；Savage et al.，1995；Wang et al.， 2012），同时更多研究发现震间期的壳幔黏弹性效应同样重要，可能会影响对闭锁深度和断层滑动速率的估算（Li et al.，2015； Zhu et al.， 2020 ； Kong et al.， 2022）。近年来相关学者开始关注震间闭锁程度的时空演化特征。例如，Melnick et al.（2017）基于GPS资料发现2010年马乌拉地震发生后，安第斯俯冲带上该地震破裂紧邻段落上的闭锁被强化。Kaneko et al.（2010）基于地震破裂模拟，揭示了凹凸体附近的闭锁程度空间特征，但未考虑壳幔流变作用的影响。考虑岩石圈流变性质时，断层系统会随着时间变化的黏性变形而释放应力，同时也会随着时间的推移而稳定地加载（Sone and Uchide，2016）。Tomoya Muramoto et al. （2023）提出应变和应力在扇断层和俯冲洋壳中具有积累趋势，其积累速率随有效黏度的变化而变化。综上所述，壳幔流变性质如何调节断层闭锁和剪应力积累仍然不清晰。

数值模型可以较好地考虑空间非均匀的壳幔流变特性（Luo and Liu， 2018； Huang et al.， 2021；Kong et al.， 2022），然而鲜有模型讨论凹凸体附近剪应力的时空变化。以往的数值模型在模拟震间变形时，尽管考虑了岩石圈流变作用，但断层往往被处理为弱化带（Luo and Liu， 2018；Huang et al.， 2021）或摩擦滑动（Kong et al.， 2022），难以考虑凹凸体及其附近的变形和应力积累。采用“负位错”的黏弹性数值模型可以考虑三维闭锁分布（Li et al.， 2015）和非均匀岩石圈流变作用，但其本质是运动学方法。因此，目前亟需发展从力学角度出发的、考虑岩石圈流变作用和凹凸体闭锁的数值模型。

针对上述问题，通过数值模拟方法构建简化的包含单个凹凸体的直立走滑断层数值模型，采用大型商业有限元MSC Marc (www.mscsoftware.com)提供的“黏接”接触算法和Maxwell黏弹性本构关系来模拟黏弹性流变作用下走滑断层凹凸体周缘闭锁程度和应力积累的时空演化。通过数值实验建立断层闭锁和应力积累之间的联系，探究强弱闭锁的阈值并给出解释，确定当闭锁程度低于某个阈值时应力积累效应相对不明显，并进一步讨论岩石圈流变结构对凹凸体附近闭锁程度和应力积累时空演化的影响。与以往基于运动学的Back-slip方法开展的数值研究（Yagi and Kikuchi， 2003； Takeuchi et al.， 2008；Li et al.， 2018a）相比，此次研究基于已有成熟的计算工具构建新的概念模型，从力学角度出发，为理解断层面闭锁和应力积累之间的联系提供新的视角。

1.   研究方法

1.1   有限元模型

1.1.1   控制方程

准静态地壳变形可以由力学平衡方程所描述：

其中${\sigma _{i j}}$为应力张量，ρ为密度，g_i为重力加速度，x为坐标，下标i和j代表张量的维度。研究中采用商用有限元软件MSC MARC 2010 (www.mscsoftware.com)对方程进行求解。

将下地壳和地幔视为各向同性黏弹性材料，在小应变理论下，对于各向同性黏弹性材料，其本构方程通常表示为遗传积分的形式（Zhu and Zhang，2013； Huang et al.，2020）：

其中σ_ij为应力张量，e_kl为偏应变张量，Δ_kl为体积应变张量，t表示当前时刻，τ表示过去时刻，δ_ij为Kronecker符号，${{G_{i jkl}}(t - \tau )}$和${{K_{i jkl}}(t - \tau )}$分别为剪切松弛核函数和体积剪切松弛核函数, d为微分算子。公式（2）中意味着各向同性Maxwell材料的偏量和体积行为是完全独立的。并且对于广义Maxwell体，剪切松弛核函数和体积剪切松弛核函数可以写成如下Prony级数形式：

其中K^∞和G^∞分别为无穷时刻的体积模量和剪切模量，${\lambda }_K^{n}$和${\lambda }_G^{n}$分别为体积和剪切松弛时间，nK和nG分别为用来表示体积和剪切松弛核函数中Prony级数项的数量，其他变量同上。广义Maxwell体是由一个弹簧（刚度用G^∞表示）并联若干条Maxwell体链条所构成，当其退化为各向同性Maxwell体时，Maxwell体链条数为1，此时nK和nG均为1，且G^∞=0。值得注意的是，在大多数地震参考地球模型中，体积模量质量因子通常比剪切模量质量因子大一个数量级。因此，只有剪切模量被视为黏弹性的，而体积模量将保持弹性（Wang et al., 2006），此时nK为0。综上所述，对于地球的Maxwell体黏弹性材料，公式（3）和（4）中的参数设置可以简化为nK=0、nG=1和G^∞=0（Huang et al., 2020）。

1.1.2   模型设置

为探究凹凸体附近闭锁程度以及剪应力积累速率情况，文中构建了概念断层模型，同时为简化计算，在模型中未考虑不同凹凸体之间的相互作用，只针对单个凹凸体开展数值模拟研究。模型由2个尺寸为400 km×400 km×400 km的三维块体构成，其中上地壳深度范围内2个板块的交界处设置为位移不连续面（断层面），上地壳深度以下单元位移连续（图1）。采用六面体单元对模型进行网格划分，在断层附近和上地壳对网格进行了局部加密，模型远场网格分辨率为10 km，断层面分辨率为5 km。模型纵向上分成3层，分别是弹性上地壳、黏弹性下地壳和黏弹性岩石圈地幔。其中，上地壳和下地壳厚度采用数值模型中的典型值20 km（Luo and Liu, 2018），岩石圈地幔厚度为360 km。考虑流变结构的影响，建立了3个模型：模型1为纯弹性模型；模型2为横向均匀、纵向分层的黏弹性模型；模型3为横向非均匀、纵向分层的黏弹性模型。表1和表2分别给出了详细的弹性、黏弹性参数设置。

图  1  模型网格图

模型坐标为右手系，其中x轴正向代表东，y轴正向代表北，z轴正向代表上

Figure  1.  Mesh of the finite element method (FEM) model

An east-north-up (ENU) system is used in the model.

下载: 全尺寸图片 幻灯片

表  1  数值模拟中使用的弹性参数

Table  1.  Elastic parameters used in the numerical simulation

分层	弹性模量/Pa	泊松比	密度/(kg/m3)
上地壳	8×1010	0.25	2700
下地壳	1×1011	0.25	3000
岩石圈地幔	1.5×1011	0.28	3300

下载: 导出CSV 

| 显示表格

表  2  模型中使用的黏度参数

Table  2.  Viscosity parameters used in the model

	上地壳块体1	上地壳块体2	下地壳块体1/Pa s	下地壳块体2/Pa s	上地幔块体1/Pa s	上地幔块体2/Pa s
模型1	−	−	−	−	−	−
模型2	−	−	1×1020	1×1020	1×1020	1×1020
模型3	−	−	1×1019	1×1020	1×1019	1×1020

下载: 导出CSV 

| 显示表格

模型的边界条件设置如下：模型表面为自由边界，模型底面设置为法向零位移、切向位移自由。在模型的侧面施加方向相反的位移约束，即在断层东盘侧边界施加南向5 mm/a的速度，断层西盘侧边界施加北向5 mm/a的速度。使用MSC Marc中提供的Glued contact方法（该方法通过绑定约束，将单元表面永久地绑定在一起，防止其发生分离或产生相对滑动），在断层面上设置一个长为5 km、宽为4 km的凹凸体。沿断层面任意方向均不能自由滑动。凹凸体中心深度为10 km，凹凸体占据的深度范围为8～12 km。除去凹凸体部分，断层其余的接触部分设置为完全解锁，均可自由滑动。

1.1.3   断层闭锁程度计算

在断层之间发生相对滑动时，如果断层两盘存在凹凸体，凹凸体的存在阻碍了断层局部区域的滑动，导致断层在凹凸体附近的滑动速率低于远场的滑动速率，二者的差异即为滑动亏损。滑动亏损K可用块体相对运动速率S和断层相对滑动速率V计算：

在模型中，断层凹凸体区域由于闭锁，在震间期不发生相对滑动，而断层的其余部分可以自由滑动。在凹凸体附近，虽然可以滑动，但其滑动速率仍然远低于块体相对运动速率。断层面上的闭锁程度ISC（Interseismic Coupling）常常用滑动亏损K与块体相对运动速率S的比值表示：

2.   实验结果

2.1   弹性模型结果

断层完全自由滑动区域不积累剪应力，但由于凹凸体的存在阻止了部分区域的相对滑动，因此凹凸体附近一定距离区域内会积累剪应力。凹凸体附近的闭锁程度分布呈现出围绕凹凸体中心环状衰减的空间样式（图2a）。凹凸体附近的断层闭锁程度和剪应力积累速率对比分析发现（图2），剪应力积累速率非零区的范围明显小于断层闭锁程度非零区的范围。换而言之，剪切应力的衰减相比闭锁程度的衰减更为迅速，这是由于与位移相比，应力是高阶量，随距离的衰减更快。当闭锁程度ISC低于0.5时几乎不积累剪应力，即闭锁程度非零的区域不一定积累剪切应力。

图  2  弹性模型的计算结果

a—闭锁程度；b—剪应力积累速率

Figure  2.  Elastic model results

(a) Interseismic coupling (ISC); (b) Shear stress accumulation rates

下载: 全尺寸图片 幻灯片

2.2   黏弹性模型结果

岩石圈的流变特性对于调节地壳变形和应力分布有着重要的作用。近年来研究发现，岩石圈流变性质不仅影响震后变形，同时调节了震间变形和应力积累 （Wang et al.，2012；Kong et al.，2022； Zhu et al.， 2020）。值得注意的是，由于地球深部岩石温度、压力条件以及应变率的差异，岩石圈的流变性质表现出横向非均匀的特点（石耀霖和曹建玲，2008）。因此，文中通过构建黏度横向均匀（模型2）和横向非均匀（模型3）2组模型来讨论在不同流变结构作用下，凹凸体附近闭锁程度和剪应力积累分布的时空演化特征。

2.2.1   横向均匀黏弹性模型

对于Maxwell体，应力松弛的特征时间λ由黏度η 和剪切模量G决定，即λ=η/G。对于模型2，由于壳幔的黏度相同，松弛时间也是相同的（2.5×10⁹ s）。不同时刻（时间分别为0.1λ、0.5λ、1λ和5λ）的凹凸体附近闭锁程度和剪切应力的空间分布如图3所示。壳幔黏弹性松弛效应导致凹凸体附近的闭锁程度和剪应力积累速率随时间进行演化，并且当时间在1倍松弛时间以内时闭锁程度较为稳定，当时间为5倍松弛时间时可以看出闭锁程度为0.2的花瓣范围有较为明显的增大。相较而言，剪应力积累速率的花瓣变化幅度较小，但同样在5倍松弛时间时可以看出花瓣有明显的增大。

图  3  均匀黏弹性模型计算的闭锁程度和剪应力积累速率的时空演化特征

a—闭锁程度；b—剪应力积累速率

Figure  3.  Spatiotemporal evolution of the ISC and shear stress accumulation rates derived from the homogeneous viscoelastic model

(a) Interseismic coupling (ISC); (b) Shear stress accumulation rates

下载: 全尺寸图片 幻灯片

2.2.2   横向非均匀黏度结构

模型3考虑了横向非均匀的黏度结构，即断层西侧的下地壳和上地幔的黏度相比另一侧降低一个数量级，设置为1×10¹⁹ Pa s。此时断层西侧壳幔的松弛时间（λ₁）降为断层东侧壳幔松弛时间（λ₂）的1/10，即λ₁=0.1λ₂=0.1λ（λ=2.5×10⁹ s）。方便起见，在图4中仍然沿用图3的时间（时间分别为0.1 λ、0.5 λ、1 λ和5 λ，其中λ=2.5×10⁹ s）。由于断层两盘松弛时间的差异，断层东盘更快进入松弛状态，断层面上的滑动亏损和剪应力分布主要受到黏度系数较低一侧的控制。

图  4  横向非均匀黏度模型计算的闭锁程度和剪应力积累速率的时空演化特征

a—闭锁程度；b—剪应力积累速率

Figure  4.  Spatiotemporal evolution of the ISC and shear stress accumulation rates derived from the heterogeneous model

(a) Interseismic coupling (ISC); (b) Shear stress accumulation rates

下载: 全尺寸图片 幻灯片

相比横向均匀黏弹性模型，横向非均匀黏弹性模型造成凹凸体附近的闭锁程度的时空演化呈花瓣范围变广的趋势，其中以闭锁程度为0.2和0.4的等值线花瓣尤为明显（图3a，图4a）。剪应力积累速率的等值线花瓣也表现出类似的规律（图3b，图4b）。以往研究也表明对于黏度横向非均匀模型，黏度值较低的一侧往往更容易吸收震间变形并积累应力（Malservisi et al.，2001； Kong et al.，2022）。

2.3   模型综合对比以及流变作用的影响机制

为进一步分析黏弹性效应对凹凸体附近闭锁程度以及剪应力的影响，追踪了特定监测点上的闭锁程度以及剪应力的时间演化特征。与剪应力分布相比，闭锁程度影响的范围较大，因此对其设置的监测点数量略有不同。断层闭锁程度采用5个监测点，而剪应力积累速率采用4个检测点，如图2所示。最终得到不同模型各个监测点处闭锁程度和剪应力积累速率随时间的变化特征，绘于图5。

图  5  监测点闭锁程度及剪应力积累速率随时间的变化

a—闭锁程度；b—剪应力积累速率

Figure  5.  Temporal evolution of the ISC and shear stress accumulation rates at the monitoring points

(a) Interseismic coupling (ISC); (b) Shear stress accumulation rates

下载: 全尺寸图片 幻灯片

对于完全弹性模型，断层的闭锁程度及剪应力积累速率不随时间而变化。在黏弹性流变作用下，由于应力的重分布而导致上地壳凹凸体附近监测点的闭锁程度随时间增大而增大，并逐渐趋于稳定。对比黏度横向均匀和黏度横向非均匀模型，黏度横向非均匀模型由于一侧的黏度较低，因此在较短的时间内进入松弛状态，所以表现出断层闭锁程度在前期增加较快。当时间足够长时，黏度横向均匀和黏度横向非均匀模型断层闭锁程度的变化趋于一致且均达到稳定状态。监测点1反映的断层闭锁程度从初始时刻的0.56增加到稳定时刻的0.7左右（图5a）。

在流变作用下，不同监测点的剪应力积累速率变化并不一致。例如，在监测点4剪应力积累速率呈随时间变化而减小的趋势，在监测点2和监测点3剪应力积累速率则表现出随时间增加而变大的特征。不同监测点的剪应力变化和断层闭锁程度变化趋势一致，均随着时间的增加而趋于稳定。当深部物质黏度较低时，流变作用对凹凸体附近的断层闭锁程度和剪应力积累的影响在较短时间内达到稳定，反之则需要较长的调整时间。当不考虑黏弹性时，监测点4的剪应力积累速率可能被高估，因为依据弹性模型估计的剪应力积累速率约为3 kPa/a，但黏弹性模型估计的稳态剪应力积累速率约为2.5 kPa/a（图5b）。

以往研究表明，由于流变效应下地壳或者上地幔的应力松弛将导致应力转移到弹性上地壳（Freed and Lin， 1998）。所以下地壳或者上地幔的黏性流会导致在地震发生之后，弹性上地壳的应力和应变的显著增加，并使其最终成为吸收原始破裂引起应变的主要场所（Freed and Lin，1998; Deng et al.，1999 ）。地球介质的流变效应不止影响震后阶段，在震间期深部物质的应力松弛同样影响壳幔系统的应力重分配，应力被转移到弹性上地壳，从而影响凹凸体附近剪应力积累和闭锁程度的时空演化特征。

2.4   震间变形分析

模型1~3计算得到的穿过凹凸体中心的水平面（深度为8 km）上北向位移速率的分布情况（图6）。距离凹凸体较远区域，位移主要受到远场构造加载的影响。在凹凸体附近，北向位移速率显著小于远场加载速率，即产生了滑动亏损。相较模型1，模型2和3由于流变作用，凹凸体附近的滑动亏损更为明显。对于模型1和模型2，可以看出北向位移速率的花瓣表现出对称性，而模型3产生的花瓣则明显表现出不对称性。在模型3中，断层西侧的下地壳、上地幔黏度较低，凹凸体对断层西侧震间变形的影响范围更大。虽然没有对弹性参数横向非均匀的情形进行讨论，但可以推测该情形下也将产生不对称震间变形。除此之外，在断层非直立而有一定倾角的情况下，同样会产生震间变形的非对称。所以当观测到非对称的震间变形时，可以侧面提供黏度结构或者断层倾角信息，例如，圣安德烈斯断层、Imperial断层和阿尔金断裂带中部均观测到断层两侧震间变形不对称现象，相关学者研究推测可能是由两侧流变系数的差异或断层倾角引起的（Lyons et al.，2002；Houlié and Romanowicz，2011；Ge et al.，2022 ）。

图  6  北向位移速率云图

a—完全弹性模型；b—横向均匀黏弹性模型（时间为5 λ，即1.25×10¹⁰ s）；c—横向非均匀黏弹性模型（时间为5 λ，即1.25×10¹⁰ s）

Figure  6.  Contour map of the northward displacement rate

(a) Elastic model; (b) Homogeneous viscoelastic model; (c) Heterogeneous viscoelastic modelThe observed depth is 8 km for all panels. The snapshots in panels (b) and (c) represent a time of 5 λ (1.25 × 10¹⁰ s).

下载: 全尺寸图片 幻灯片

2.5   断层显著闭锁阈值的探索

Li et al. （2023）提出用0.5作为断层显著闭锁阈值，认为超过该值可视为中强闭锁，低于该值为弱闭锁状态。但这个阈值的提出缺乏力学依据，因此笔者试图通过力学模型来找到一个临界的闭锁程度值，当超过该值时断层面剪应力积累明显。由于闭锁程度是介于0～1的无量纲量，为了更好地对比分析闭锁程度和剪应力积累的空间分布样式，在用剪应力积累速率（dτ/dt）除以剪应力积累速率的最大值(dτ/dt)_max后得到一个归一化的剪应力积累速率因子（Φ_τ），由闭锁程度（ISC）和归一化剪应力积累速率因子（Φ_τ）的等值线对比可知（图7），凹凸体附近的闭锁程度非零区域比剪应力积累非零区域的范围更大。因此，存在一部分区域的闭锁程度不为0，但剪应力积累值很低，几乎可以忽略不计。

图  7  由弹性模型得出的凹凸体附近闭锁程度和剪应力归一化因子等值线

红线表示ISC等值线；黑线表示Φ_τ等值线，等值线间隔为0.2，由凹凸体向外扩展的4根黑线分别为Φ_τ=0.8、0.6、0.4、0.2

Figure  7.  Contour map of the ISC and shear stress normalization factor near an asperity derived from the elastic model

Red lines represent the ISC contour lines at an interval of 0.1. Black lines represent the Φτ contour lines at an interval of 0.2. The four black lines extending outward from the asperity were 0.8, 0.6, 0.4 and 0.2.

下载: 全尺寸图片 幻灯片

根据完全弹性模型计算结果，当断层闭锁程度小于0.3的情况下，剪应力的积累速率小于最大值的0.2，即Φ_τ<0.2，此时剪应力的积累效应相对不明显。当闭锁程度ISC=0.5，对应的Φ_τ约为0.4。根据横向均匀黏弹性模型，在一倍松弛时间内，黏弹性模型和弹性模型计算得到的闭锁程度ISC和归一化剪应力积累速率因子Φ_τ的等值线样式接近。当时间足够长，产生了明显差异：闭锁程度为0.4对应的区域上剪应力积累速率就会降低到最大值的0.2以内。通过对比图7和图8发现，对于模型3断层一侧的黏滞系数降低一个数量级，因此该模型更快达到应力松弛状态。根据图8，当时间足够长（5 λ），在闭锁程度为0.5对应的区域上，在凹凸体上方部分区域Φ_τ<0.2。由此可见，闭锁程度的阈值受到流变结构影响，并且和当前时刻在地震循环中所处的阶段有关。例如，在震间早期，黏弹性效应尚不明显，作为一阶近似，闭锁程度0.5对应的归一化剪应力积累速率因子Φ_τ约为0.4，将其作为中强闭锁程度的阈值不失为一种简化的方法（Li et al.， 2023）。

图  8  由黏弹性模型给出的凹凸体附近闭锁程度和剪应力归一化因子的时空演化特征

a—均匀黏弹性模型；b—非均匀黏弹性模型红线表示ISC等值线；黑线表示Φ_τ等值线，由凹凸体向外扩展的4根黑线分别为Φ_τ=0.8、0.6、0.4、0.2

Figure  8.  Spatiotemporal evolution of the ISC and shear stress normalization factor near an asperity derived from the viscoelastic model

(a) Homogeneous viscoelastic model; (b) Heterogeneous viscoelastic model Red lines represent the ISC contour lines at an interval of 0.1 Black lines represent the Φ_τ contour lines at an interval of 0.2. The four black lines extending outward from the asperity were 0.8, 0.6, 0.4 and 0.2.

下载: 全尺寸图片 幻灯片

3.   讨论

3.1   鲜水河断裂带模型验证

为验证“黏接接触”考虑断层锁固的震间变形计算方法，选取鲜水河断裂带炉霍至康定段作为应用实例，计算在断层闭锁调节下的速度场。鲜水河断裂是一条分割松潘−甘孜地块和川滇地块的大型左旋走滑断裂，自西向东呈现出滑动速率的过渡。西北段滑动速率为10～20 mm/a，东南段的滑动速率为5～9 mm/a。鲜水河断裂带不但穿越人口稠密区，还与川藏铁路等国家重大基础设施交汇，该断裂的地震危险性长期受到关注。历史上鲜水河断裂带地震频发（Wen et al.，2008 ）。白永健等（2019）统计了青藏高原东南缘的地震发生情况以及由此引发的地质灾害，该地区受灾范围甚广，受灾人数众多，因此该区域的地震危险性值得关注。2022年泸定M_W 6.7地震再次引起学者对鲜水河断裂的关注，讨论该断裂的运动学特征和应力积累对于地震灾害评估具有重要意义。

针对鲜水河断裂应用实例建立了三维有限元模型，模型设置如表3所示，将模型分为上地壳、下地壳、上地幔三层，上地壳设置为完全弹性，下地壳和地幔设置为黏弹性，断层闭锁设置在断层的上地壳处。基于大地测量的鲜水河断裂闭锁程度反演揭示了该断裂广泛分布的浅部蠕滑（Qiao and Zhou，2021）和潜在的凹凸体（Jiang et al.，2015 ），也为数值模型的凹凸体设置提供了依据。如图9a所示，将闭锁程度高于0.5的区域设置为凹凸体(图9a中黑线框圈中部分)，即“黏接接触”。鲜水河断裂朱倭到炉霍段和康定附近的高闭锁区主要集中在浅部，而炉霍以南的段落高闭锁区主要集中在深部。

表  3  鲜水河断裂带数值模拟中使用的参数

Table  3.  Parameters used in the numerical simulation of the Xianshuihe fault zone

分层	深度/km	弹性模量/Pa	泊松比	密度/(kg/m3)	黏度/（Pa s）
上地壳	18	8× 1010	0.25	2700
下地壳	18	1× 1011	0.25	3000	1×1020
岩石圈地幔	64	1.5× 1011	0.28	3300	1×1020

下载: 导出CSV 

| 显示表格

图  9  鲜水河断裂带邻区的震间变形

a—鲜水河断裂带闭锁程度分布图（Jiang et al., 2015；黑色划线区域为此次模型中设置的凹凸体的位置）；b—地表速度模拟值和GPS速度的比较（带误差椭圆的红色箭头表示模拟的地表地震间速度，带误差椭圆的蓝色箭头表示模型内部的GPS速度观测值，带误差椭圆的白色箭头表示模型外的GPS速度观测值）

Figure  9.  Interseismic deformation near the Xianshuihe fault

(a) Map of the locking degree of the Xianshuihe fault (Jiang et al., 2015; The area outlined by the black line represents the asperities used in this model.) (b) Comparison of simulated surface and GPS velocities (The red arrow with an error ellipse represents the simulated interseismic velocity, blue arrow with an error ellipse represents the GPS velocity, and white arrow with an error ellipse represents the GPS velocity outside the model domain.)

下载: 全尺寸图片 幻灯片

模型底部无法向位移，侧面节点的速度边界条件根据GPS速度（Wang and Shen，2020）插值确定，如图6b所示。模型中下地壳和地幔黏度均按照1×10¹⁹ Pa s施加，弹性参数参考已往数值模型（Luo and Liu，2018），下地壳的松弛时间约为 8 年。由于模型排除了2022年泸定M_W 6.7地震和1973年四川炉霍M 7.6地震，最近的一次地震为1955年康定折多塘7.5级地震，距今时间约为80 年。所以计算时间选为80年，在数倍松弛时间作用下，最近一次地震引起的应力调整和变形扰动可以近似忽略，作为简化，模型未考虑历史地震的影响。

模拟的地表运动速率与GPS观测速率基本一致，如图9b所示。同时，为了定量评价观测误差，文章计算了模拟结果和观测结果之间的误差均方根为1.8 mm/a。根据图9b显示的模拟吻合度以及较低的误差均方根说明了断层闭锁和参数设置的合理性。但值得注意的是，距离断裂较远的一个观测点数值模拟的地表运动矢量和观测速度矢量之间有较大的夹角，这可能是因为模型只考虑了一条主要断裂，一些次级断裂未被考虑。而距离主断裂较远的观测点可能主要受到次级断裂的影响。同时，数值模型的构建势必要进行合理地简化，这些简化也不可避免造成模拟结果和观测结果之间的差异。

3.2   二元闭锁状态和闭锁程度之间的联系

基于空间大地测量观测数据，采用Back-slip方法的运动学反演给出的闭锁程度往往是一个连续分布的值，范围在0～1之间。而在力学上，断层闭锁状态只有闭锁和非闭锁两种情况：在闭锁区，断层无滑动；非闭锁区域，断层可以滑动。通过数值计算，二元闭锁状态下产生的闭锁程度在空间上同样是连续的，在0～1之间分布，并非凹凸体上为1，离开凹凸体立刻阶跃为0的样式。这是由于块体介质在内部是连续的，因此在每个块体内部滑动速率也应该是连续的，不会出现阶跃。导致凹凸体附近的滑动速率是连续变化的，中间存在一个缓冲的地带，这样就形成了滑动亏损过渡区，也就是闭锁程度的过渡区。同时，根据计算结果，闭锁程度高值区（例如，闭锁程度大于0.8）可以较好地勾勒出凹凸体的形状。而闭锁程度低值区的分布范围则取决于所处的时间和黏度结构。

计算结果揭示了在走滑断层中不能简单地把运动学的闭锁程度和力学闭锁混为一谈。对于俯冲带，情况则变得更为复杂。例如，实验室摩擦实验和热结构研究 （Byrne et al., 1988；Vrolijk, 1990；Oleskevich et al., 1999）认为在俯冲带浅部，因其摩擦性质表现为速度强化（velocity strengthening），而处于无震蠕滑状态。在这种条件下，俯冲带浅部可以被认为是摩擦解锁断层，但其运动学闭锁程度一定很低吗？答案是否定的。在俯冲带浅部，虽然没有摩擦闭锁，但是却有可能处于附近闭锁区域产生的应力阴影区，从而可能导致高的运动学闭锁（Almeida et al.，2018； Lindsey et al.，2021 ）。

3.3   模型限制

弹簧−滑块模型（即单自由度弹性动力系统）可以较好解释岩石摩擦的实验结果，并持续发展为现今应用广泛的速度−状态依赖的摩擦本构定律（Dieterich，1979；Carter et al.，1981；Ruina，1983）。在此研究中，未考虑摩擦系数的演化，而是采用“黏接接触”算法强制地让凹凸体难以自由滑动。在模型中，由于主要关注震间期凹凸体闭锁程度和剪应力积累速率的时空演化，并不关注地震的孕育成核过程，作为一阶近似暂不考虑速度−状态依赖的摩擦本构定律。真实断层往往存在多个凹凸体，例如，基于黏弹性分层模型的断层闭锁程度反演（Jiang et al.，2015）表明鲜水河−安宁河−则木河断裂系统上存在6个潜在凹凸体。最新的青藏高原东南缘主要断层滑动速率和闭锁系数反演研究，识别出数十个闭锁系数>0.5的强耦合断层段（Li et al.，2023）。目前，此研究仅考虑单凹凸体，凹凸体之间如何相互作用仍未得到很好的理解，将在下一步工作中进行重点探究。同时，闭锁程度以及地应力的时空变化，并不局限于文中讨论的流变作用。例如，2008年汶川地震、2010年智利Maule地震的发生，改变了地震破裂面附近断层面的闭锁程度以及应力分布（Melnick et al.，2017; Li et al.，2018b; 丰成君等，2018）。对于多凹凸体系统，当其中一个凹凸体破裂，如何影响断层面上的断层闭锁和应力分布？下一步工作可以通过考虑“黏接接触”的失效来近似模拟这一过程。篇幅限制，鲜水河断裂的应用仅给出了地表运动速度场，未来的工作中还应进一步分析断层面的应力积累分布，同时流变参数和断层闭锁的不确定性还需要被进一步讨论，以保证结果的鲁棒性。

4.   结论

通过建立包含接触单元的三维黏弹性有限元模型，讨论流变作用下直立走滑断层上单一凹凸体附近的应力积累速率和闭锁程度的时空演化特征，得到以下主要结论。

（1）由于介质的连续性，凹凸体邻近区域虽然不闭锁，但其滑动速率仍低于块体运动速度，从而导致闭锁程度的空间样式表现为由凹凸体向外逐渐衰减。

（2）黏弹性效应调节凹凸体附近的变形，并进一步影响凹凸体附近闭锁程度的时空演化，断层闭锁程度等值线空间范围随时间增大，且在同一位置闭锁程度也随时间加载而不断增大。

（3）凹凸体附近的断层闭锁程度与剪应力积累速率的变化趋势相似，闭锁程度可以近似作为衡量剪应力积累速率的指标，当不考虑黏弹性效应时，可以取0.5作为中强和弱闭锁的分界值，低于该值剪应力积累不明显。

（4）考虑鲜水河断裂带炉霍至康定段的空间非均匀闭锁分布，模拟得到的地表速度与GPS观测吻合较好，验证了方法的可靠性。