慢滑移是一种发生在构造断层的类地震事件，有时会与非火山震颤同时出现。它们同时出现的现象被称为间歇性震颤和滑移（缩写为ETS)。大多数慢滑移会持续几天到几年之久，比典型地震的持续时间要长得多。这类事件往往不会产生可检测的地震信号，因而只能通过GPS等大地测量仪器来观测。部分研究认为，慢滑移会影响浅层孕震区的应力加载，使锁闭区域更趋向破裂（Dragertetal,2001），从而在某些情况下诱发浅层地震。大多数慢滑移发生在俯冲带中力学和水热条件发生改变的深度。Brown等（2005）指出，部分慢滑移事件的产生与断裂带中的流体有关。因此，很多研究认为断裂带的高孔隙压力（导致低有效正应力）可能是慢滑移产生的一个重要机制。除了俯冲带之外，慢滑移也可能发生在转换断层等其它构造环境中。

在过去的几十年里，研究者们从世界各地收集了大量的慢滑移数据。基于这些观测结果，许多地球物理学家从不同角度对慢滑移进行了数值模拟。但是，由于计算能力的限制，大多数模拟中的断层几何结构是基于较为简单的二维平板模型。直到近些年来，基于三维模型的慢滑移数值模拟才得以进行，然而，三维模型对慢滑移的数值模拟会造成何种影响并不十分确定。鉴于此，本文拟使用三角单元对断层的几何形状进行离散，在一个三维的、真实的断裂带几何模型上模拟构造运动所发生的过程，并检验断层几何形状对构造运动所造成的一系列影响。在模拟中，我们使用从实验室中得到的速率-状态依赖型摩擦律和一种准动态算法，选择三种不同的断层几何构造模型进行数值模拟，并对其模拟结果进行对比与讨论，以探究断层模型的几何形态对数值模拟中慢滑移特征数值的影响。

1、阿拉斯加中南部俯冲带构造背景

阿拉斯加俯冲带位于阿拉斯加中南部，北美板块与东太平洋板块交界处（图1）。太平洋板块和亚库塔特地块以55mm/a左右的速度相对于北美板块向西北方向俯冲（Finzeletal,2011），形成了阿拉斯加俯冲带。该区域内所发生的地震较多、释放能量较大，是全球构造运动较为活跃的地区之一。1964年3月28日，阿拉斯加中南部发生了有现代地震记录以来的第二大地震，其矩震级为MW9.2，破裂从震中沿俯冲带向西南方向传播了600—800km后才停止。2002年11月3日，阿拉斯加地区发生MW7.9德纳利大地震，这次大震是近150年来北美板块内发生的最大走滑地震，地表破裂区域长达340km。此次地震的破裂方式以右旋走滑为主，并且沿破裂方向产生了异常强烈的远端效应（例如诱发地震）。

除地震之外，世界上最大的一次慢滑移也发生在阿拉斯加地区，矩震级约为7.5(Fuetal,2015）。在1995年至2017年期间，研究人员发现，在上库克湾和下库克湾至少发生了两组重复的特征不同的慢滑移事件组，其中：一组发生在上库克湾附近的安克雷奇（Anchorage）城附近，矩震级约为7.2—7.5，在1998—2001年（Ohtaetal,2006）和2009—2013年（Fuetal,2015）各发生一次，平均持续时间为3—4年，间隔约为11年；另一组慢滑移则发生于下库克湾的西南部，是在2010—2012年间发生的一次矩震级约为6.9的慢滑移（Weietal,2012）。基于2000年10月安装的连续GPS台站所记录到的数据，本文作者认为一个类似的慢滑移可能发生在2002—2004年之间，而同组中另一个与2010—2012年事件的大小和空间范围均类似的慢滑移则始于2016年初（Lietal,2018）。由此可见，自二十世纪九十年代以来，下库克湾慢滑移的特征与上库克湾有较大的差异。这些不同的持续时间、强度和空间范围表明，在阿拉斯加中南部俯冲带上至少存在两个区域，其内部的慢滑移特征并不相同。

阿拉斯加中南部俯冲带因其独特的地理位置和复杂的构造环境对该区域下方太平洋俯冲板块的几何形态、俯冲深度以及增生块体在深度方向上的几何形态、构造运动等问题（Fuisetal,1991）的研究具有非常重要的意义。

图1包括两个典型GPS台站（ATW2和AC06）的阿拉斯加中南部地区示意图

2、数值模拟的模型设置

2.1控制方程

在本文的数值模拟中，主要参数由两组控制方程所决定。第一组控制方程描述了断层破裂强度与有效正应力、摩擦系数之间的关系，也即摩擦律，这一组关系决定了数值模拟中每一个网格达到破裂所需剪应力的大小。第二组控制方程描述了每个网格的剪应力与断层滑动之间的关系，这一组关系决定了数值模拟中每一个网格的剪应力演化过程。将两组控制方程结合起来，即可对断层构造运动的演化情况进行数值模拟。

在本文使用的速率-状态依赖型摩擦律中，断层破裂强度τ定义为（Dieterich,1979;Ruina,1983):

式中：V为断层滑动速率；θ代表了滑动的历史，单位与时间单位相同；μ为摩擦系数，是一个取决于V和θ的状态变量；a和b为无量纲的速率-状态依赖型摩擦律参数；Dc为特征滑动距离，表示系统滑移速率受到扰动后从一个稳态速度变化至下一个稳态速度所需要的滑动距离；μ0为参考速度V0下的稳态摩擦系数；为有效正应力，，其中σ为正应力，p为孔隙压力。这里我们使用最常用的速率-状态依赖型演化定律之一的“慢度”定律，即

该定律允许摩擦在固定的粗糙接触面上变化，并得到了实验室实验的广泛支持。

另一个广泛使用的描述状态量变化的速率-状态依赖型演化定律是“滑动”定律。在状态变量的演化方程中，为了使方程有意义，滑动不能为零。在使用相同代码对卡斯卡底的慢滑移进行模拟的一项研究中，Li和Liu(2017）既使用了“滑动”定律又使用了“慢度”定律，其研究结果表明，“慢度”定律在参数选择方面具有更大的灵活性。因此，本文选取“慢度”定律。

除速率-状态依赖型摩擦律之外，本文研究的第二组非常重要的演化关系为不同网格间的剪切应力与断层滑动之间的关系。由于断层被设定嵌入在一个弹性体中，因此任意一个网格的运动都将导致其它网格剪应力的变化。在本文的数值模拟中，关于网格位错与其它网格剪应力的变化关系，我们采用Rice(1993）的准动态模型，其离散形式为

式中：τi(t）和δi(t）分别为第i个单元的剪切应力和滑移量；η为辐射阻尼因子，η=G/2cS，式中G为弹性剪切模量，cS为剪切波速度；刚度矩阵Kij表示第i个单元上剪切应力由于第j个单元网格的位错所导致的变化，这个格林函数Kij是在弹性半无限空间中计算的（Okada,1992），由Stuart等（1997）首次在三角离散元中实现，本文中以Meade(2007）开发的Matlab代码为基准进行计算.辐射阻尼因子的引入是为了防止地震时滑移速度会无限制的增大（Rice,1993）。Liu和Rice(2007,2009）认为这种准动态方法适用于完全锁定的逆冲地震带的慢滑移过程模拟。

2.2模拟参数

在本文的数值模拟中，若参数值与深度无关，例如特征滑动距离Dc和网格的大小，那么为了避免这些参数对结果的影响，它们均将被设置成相同的数值。所有与深度相关的参数，如有效正应力、本构方程无量纲参数a和b、沿倾向速度弱化区的宽度W，与深度的关系演变在每一次模拟中均相同。除Dc和，与Dc成比例）外，大多数参数值的选取均与实验室所得的实验值或者全球地震活动性研究、理论所推断的数值相一致。

无量纲参数a-b在速率-状态依赖型摩擦律的框架下决定了断层的滑动稳定性：当a-b>0（速度强化）时，断层会连续地、抗震地滑动，这种情况称为蠕变；当a-b<0（速度弱化）时，滑动可能不稳定或是有条件的稳定（Ruina,1983）。而在此情况下，断层具体的行为则会取决于另一个无量纲参数W/h*，其中W为沿倾向的慢滑移区域（同时满足速度弱化、低有效正应力的条件）的宽度，而h*是不稳定断层滑移的最小破裂尺度，也被称为特征成核尺寸。已有研究表明，长度比W/h*是决定俯冲带和转换断层发生抗震还是类地震滑动事件的一个关键参数。对于较大的有效正应力或较宽的速度弱化区，当W/h*>1.0时，断层的滑动方式是类地震的；当0.6

在室内岩石滑动实验中，摩擦参数a-b往往取决于温度、岩石类型和加载条件（Blanpiedetal,1998）。通常采用断层的热模型将实验室所测得的a-b与温度的关系转换为与断层深度有关的a-b关系。a-b=0的深度定义了类地震滑移到非震滑移的过渡线，对于不同类型的岩石，往往会有不同的过渡等温线，例如湿花岗岩的过渡等温线是350℃（Blanpiedetal,1998），在超临界水条件下湿辉长岩泥的过渡等温线是500℃（Heetal,2006）。在本文的数值模拟中，我们采用Abers等（2006）通过接收函数与温度和应变率依赖型橄榄岩流变律所获取的阿拉斯加俯冲带热模型。为了与GPS的观测结果保持一致，我们设置一个由（0℃，0.0035),(100℃，-0.0035),(344℃，-0.0035),(448℃，0.001）这几个点所连成温度与a-b的关系。我们所采用的参数a-b的值并不能严格地反映实验室中测量到的某种岩石类型。相反，我们允许a-b有适当的变化以符合一些观测数据和理论研究的限制（Lietal,2018），且使a-b保持在与实验室所测量数据相同的数量级（约0.001）。

特征滑移距离Dc在实验室中测得的数值往往处于1—100μm的范围内，但Beeler(2004）指出，在实际断层的大尺度破裂中，Dc是随着断层滑动距离而增长的。因此，在本文的数值模拟中，取Dc=19.25mm（此时h*=7km）以确保W/h*处于一个合理的范围，同时也与其它研究中的Dc值非常接近。

已有研究显示，慢滑移事件的发生机制通常与断裂带的高孔隙压力（导致低有效正应力）有关。在我们的模拟中，只有设置了低有效正应力才能产生在特征（例如事件间隔、持续时间、等效震级等）上较为符合GPS观测结果的慢滑移（Liu,Rice,2005,2007）。这个现象与该地区远震接收函数的结果相一致（Kimetal,2014）。在该深度上，孔隙压力的增加可能与靠近莫霍面大陆的矿物脱水有关。为保证慢滑移的发生区域与实际GPS观测到的事件发生区域大致一致，我们在模拟中对39—60km的深度范围设置低有效正应力为20MPa，以便在研究区产生与观测到的慢滑移发生区域相一致的慢滑移，这也与之前的一些研究中低有效正应力的量级相一致。为了避免在GPS观测到的慢滑移区域（深度为39—60km）以外产生慢滑，0—39km深度区域的应力随深度从0MPa开始迅速增加至600MPa，以此锁定住比39km更浅的断层；而比60km更深区域的应力也设置为600MPa，以防止事件沿断层向更深的区域进一步发展。同理，我们将浅于39km、深于60km的断层的Dc值也设为2000mm。

虽然沿阿拉斯加—阿留申弧向西的俯冲速度逐渐增加，但在本文的研究区域内，俯冲速度约为55mm/a，变化并不大，故设置剪切模量G=30GPa(Suito,Freymueller,2009），泊松比ν=0.25；将慢滑移的阈值定义为断层上的最大速度达到板块俯冲速度的三倍，这也与该区域的GPS检测阈值较为接近。

在之前的一些数值模拟研究中，一部分研究将模型设置为弹性，而另一部分则在模型设置中考虑到了黏性，这些研究大都成功地再现了一些观测到的慢滑移特征。本文的目标是探究断层几何构造对慢滑移数值模拟的影响，因此，我们无需纠结黏性和弹性的选择，将模型设置为完全弹性。数值模拟关键参数的具体数值详见表1。

表1阿拉斯加慢滑移数值模拟的关键参数

2.3断层几何构造模型

在过去的几十年中，许多研究人员试图通过地震的位置来描绘阿拉斯加俯冲断层的几何形状，从而创建了很多模型。本文选取了其中应用广泛的两个模型—Slab1.0模型（Hayesetal,2012）和Li等（2013）的模型，后者简称为Li模型。Slab1.0是主要俯冲带三维几何模型的汇总，其中每个俯冲带的模型都是基于对历史地震目录、CMT解决方案、活动地震剖面、全球板块边界、测深和沉积物厚度信息等独立数据集进行加权平均而得。Li模型的确定则基于数千次微震的重定位。2006年至2009年期间，多学科俯冲带观测项目在基奈半岛部署了一个密集的临时地震台网，通过震源和速度结构的联合反演会自动检测出发生的地震。随后他们使用双差定位法hypoDD(Waldhauser,Ellsworth,2000）重新定位这些地震从而确定俯冲带的几何形状。此外，本文还使用了一个倾角为13°的二维平面俯冲带几何模型作为控制模型来进行对比。图2给出了上述三个俯冲带模型的研究区域。

在以往的慢滑移数值模拟中，囿于复杂度与计算成本，断层的几何模型大多是基于离散矩形网格的二维平面模型。本文采用三角单元对断层的几何形状进行离散，这使得我们能够在一个三维的、更接近实际情况的断裂带几何模型上模拟所发生的构造运动过程，并检验断层几何形状对构造运动所造成的一系列影响。为了简化计算，以（151°W,60°N）为原点将研究区域顺时针旋转50°，使x轴方向与断层走向大致一致，y轴方向与俯冲板块倾角方向在地表的投影大致一致，z轴为深度。考虑到这一地区慢滑移的具体发生位置，研究区域设定为（143°W—157°W,55°N—64°N），该区域内包含了绝大多数有记录的慢滑移。本文选取了x轴从-150km到300km,y轴从-180km到200km,z轴从20km到100km作为最终的研究区域。首先将研究区域内的俯冲断层离散为3万多个边长为2km的正方形单元格，随后将每个单元格划分为两个全等的三角形来更好地拟合俯冲带的真实形状。虽然利用有限厚度的破裂带对慢滑移的一系列模拟取得了很好的效果，但在本文的模拟中，为了聚焦于俯冲带的几何形态对慢滑移特征的影响，我们并未将上覆和俯冲板块的影响引入模拟中，这一处理方式与很多研究都相同。

图2三个俯冲带几何模型的研究区域（彩色）.粉色线条为等深线

3、模拟结果

为了解断层几何形状对慢滑移时空演化的影响，我们采用两种不同的三维俯冲断层模型和一个倾角为13°的二维平面俯冲断层模型来进行慢滑移模拟。图3给出了三个模型在俯冲断层上速度弱化、高孔隙压力（低有效正应力）区域（慢滑移区域）上的平均速度，黄色区域表示慢滑移所发生的时空位置（x轴坐标）。由图3可以看出：平面断层在研究区域左右两侧的对称位置和中间位置均产生了非常相似的重复的慢滑移（图3c）；而在两种非平面模型中，慢滑移行为则较为复杂。Li模型中出现的慢滑移往往集中在中间位置，且大多不相似（图3b）；而在Slab1.0模型中，研究区域的各处均会出现慢滑移（图3a）。尽管在Slab1.0模型中，其东部发生的慢滑移远远大于西部发生的慢滑移，但与另外两个模型相比，Slab1.0模型中发生在东部和西部的慢滑移之间的相互影响较小，事件内部的相似性和周期性更好。

除此之外，我们还计算了上、下库克湾中两个典型GPS站点（AC06和ATW2）的合成地表位移，并计算提取了慢滑移事件地表位移的平均间隔、持续时间、振幅和最大滑移速率等特征值以分析三种模型在地表的位移差异。结果如图4和表2所示，可见：三个模型在AC06台站的四组特征较ATW2台站更加接近；在平板模型中，两个台站的慢滑移地表特征值非常接近，而且如前文所述，平板模型的结果非常对称，造成这种台站间微小差距的唯一原因是它们与对称轴的距离不同。在Li模型和Slab1.0模型中，尽管区域平缓程度以及因此导致的事件大小程度有所不同，但普遍而言，更快、更大的事件发生在更平缓的区域（例如Li模型中的AC06台站和Slab1.0模型中的ATW2台站）。与Slab1.0模型中发生的大而快的慢滑移相比，在二维平面模型和Li模型中发生的慢滑移事件相对更平缓。与Slab1.0模型中发生的大而快的事件相比，在平板模型和Li模型中发生的事件相对比较平静。总的来说，对于所有的结果，更快的速度和更大的振幅对应于更长的平均间隔和更短的持续时间。

图3基于Slab1.0模型（a),Li等（2013）模型（b）和二维平板模型（c）的慢滑移区域（高孔隙压、速度弱化）内20—80年平均滑动速率的模拟结果

图4基于三个模型两台站y方向的地表合成位移（线性趋势已经去除）和速度

(a)AC06台站；（b)ATW2台站；（c)ATW2台站移除了Slab1.0模型的曲线

表2AC06和ATW2台站的慢滑移地表特征

4、讨论与结论

我们的结果表明，慢滑移区域沿倾角的宽度W沿走向的变化对三维俯冲断层的数值模拟具有很大的影响。由图3可知，无论基于哪种模型，在俯冲带相对平坦的区域（更平坦的区域将导致更宽的慢滑移区域）更有可能发生尺度和震级更大、滑动速率更快的慢滑移。一个非常明显的例子是在Slab1.0模型中东部平缓断层所发生的非常大的慢滑移。不同深度的等深线在地表投影之间的距离反映了断层几何形态的平缓情况。在Slab1.0模型中，39—60km等深线之间的距离在俯冲带东部比西部要大得多，这表明在Slab1.0模型中，断层的东部比西部要更平坦。Slab1.0模型中东部的ATW2站的综合地表位移较大（图4，表2所示）同样证明了这一点。我们通过数值模拟所得出的这一结论与其它研究一致。此外，还有一些研究表明，长度比W/h*对于发生在俯冲带和转换断层的事件到底是震颤滑动还是非震滑动起到了关键的作用。对于一个细长的断层，W/h*也会决定慢滑移的分段情况（Liu,2014）。在大多数此类数值模拟中，h*往往是固定的，因此W为决定慢滑移时空演化行为的唯一断层几何参数。换言之，若想在断层中重现慢滑移现象，低有效正应力区域的宽度需要足够大，以使不稳定滑动能够成核，但也不能太大，以免使事件进一步发展成为地震。

本文的工作还表明，数值模型中断层几何模型的选择将对模拟的结果产生非常重要的影响。在数值模拟中，俯冲断层的几何形状既影响地震成核也影响慢滑移。在卡斯卡底北部，倾斜角度较小的平坦断层上，慢滑移倾向于累积较大的滑移量。本文的模拟结果为后续模拟提供了很好的指导，即采用更真实的断层模型是非常必要和重要的。

李昊天,周仕勇.断层几何形态对阿拉斯加中南部俯冲带慢滑移特征的影响[J].地震学报,2019,41(06):681-694+808.

基金:国家自然科学基金(41674047)资助.