1 引言

地球物理观测表明，深度大于70 km的深源地震多发生在俯冲带附近，地震数量随深度继续增加而减少.但在过渡地幔深度(410~670 km)附近先突然增多，在约600 km深度达到峰值，然后在下地幔深度(670 km)附近又突然减少并在约700 km深度几乎消失(Frohlich, 1989, 2006).发生在地壳内的浅源地震一般是脆性破裂，而在深源地震发生区域，由于温度和压力的增加，岩石可能发生塑性变形，其震源机制仍有较多争议.深震震源多显示剪切破裂和双力耦特征(Scholz, 2002).目前，比较流行的深源地震机制主要有三种，即脱水脆裂、反裂隙断层作用和剪切热失稳.脱水破裂机制主要能解释深度小于300 km的深震.反裂隙断层作用在深部受青睐，主要指的是物质在非静水压条件下发生相变时产生体积变化和不稳定的剪切变形，从而形成沿最大应力方向的断裂，但该机制主要跟亚稳态相变有关(Kirby, 1987; Green and Burnley, 1989; Green et al., 1990; Burnley et al., 1991; Kirby et al., 1991, 1996; Green and Zhou, 1996).不过，虽然理论上存在观测到亚稳态相变的可能性(叶国扬等，2008)，但是目前还没有来自地球物理观测数据的支持.因此，亚稳态相能否在该深度下存在仍有很大争议(Koper and Wiens, 2000; Marton et al., 2005; Jiang and Zhao, 2011；王曙光等，2011).尤其是在600 km深度附近，深源地震突然增多，这更是亚稳态相变无法解释的.Karato等(2001)通过分析不同蠕变模型和剪切热失稳解释了最深地震的几个矛盾.剪切热失稳理论认为，在塑性变形失稳的基础上，绝热剪切会引起物质熔融，使物质变形速率加快而产生热量积累，积累的热量进一步促进变形发生，最终又导致物质塑性变形的失稳.剪切热失稳理论很好地解释了深源地震的发生(Hobbs and Ord, 1988; Karato et al., 2001; Wiens, 2001)，但是一直没有实验支持.近年的地震观测结果(Chen et al., 2014; Deseta et al., 2014; Chen and Wen 2015)也越来越倾向于剪切热失稳机制.最近，Liu和Zhang(2015)从宏观应力的角度又提出了一种新的观点，认为应力来源于地壳和上地幔的不同物质在俯冲带条件下体积收缩产生的差异，地震在深部的表现与俯冲板块的俯冲历史有关.该观点能较好地解释东西太平洋和大西洋边缘深源地震的不同区域特征, 这使得深源地震机制的研究越来越多样化.但由于俯冲带环境复杂，深源地震很可能并不是只有一种机制作用的结果(邵同宾和嵇少丞，2015)，而是受多种机制共同作用，并且仍有多个问题需要进一步探讨(干微等，2012).

解释深源地震的几种机制有一个共同特点，即需要一个失稳过程.在这个失稳过程中，应力的变化是一个直接参与量，因此，解释深源地震机制的主要任务是寻找一个在俯冲带的高温高压等复杂环境下能产生破裂的机制.本文通过研究不同深度两相矿物岩石模型在相应条件下的蠕变过程，分析了岩石模型内部剪切应力的发展变化，结果显示林伍德石、石榴石两种矿物的蠕变差异以及林伍德石快速蠕变特点，在600 km深度附近极大地促进了剪切应力的发展，可能是该深度地震突然增多的重要原因.

2 方法和计算 2.1 矿物模型

我们采取俯冲带不同深度的主要矿物组合建立两相矿物模型，计算矿物模型在一定温度和差应力下剪切应力随时间的变化，分析矿物模型对俯冲带不同深度剪切应力发展的影响.首先，参考地幔不同深度主要矿物的含量(Ringwood, 1975; Irifune and Ringwood, 1987; Mitra, 2004)，建立二维两相矿物分布模型，在200~410 km深度为模型1，组成为：橄榄石-单斜辉石；410~520 km深度为模型2，组成为：瓦兹利石-单斜辉石；520~670 km深度为模型3，组成为：林伍德石-石榴石.三个模型中矿物体积比都为60%~40%，如图 1所示，模型中共含20×20块正方形，每个正方形代表一个矿物单元，其中浅灰色方块与深灰色方块个数比约为6:4.对模型施加温度和差应力条件，利用流变方程，计算模型内部剪切应力随时间的变化.分析剪切应力在模型内部的发展变化.

模型的应力应变本构方程为弹-黏塑结构，总应变ε_ij分为弹性应变ε_ij^e和非弹性(黏塑性)应变ε_ij^p两部分，其中弹性部分的应力应变关系满足胡克定律(σ_ij为应力张量，c_ijkl为应力应变关系的本构系数)：

(1)

由于俯冲带主要矿物的变形机制(尤其在低温区域)从理论和实验上都以幂指数形式的蠕变为主(Shimojuku et al., 2009; Xu et al., 2003; Kawazoe et al., 2013)，因此非弹性应变部分采用幂指数流变方程：

(2)

其中，为应变率, A为常数，σ为差应力, n为应力指数, H是蠕变活化焓, R是气体常数，T是绝对温度.橄榄石、单斜辉石、瓦兹利石、林伍德石和石榴石在相应稳定压力和干的环境下(本研究中不考虑水的影响，采用矿物不含水条件下的流变方程)流变方程见表 1.密度及弹性模量见表 2.

2.2 初始和边界条件

考虑到俯冲带深源地震的应力降在几十到上百个兆帕(Wiens, 2001)，而俯冲带环境应力可能达到几百兆帕甚至更高(Liu and Zhang, 2015)，因此我们选取100 MPa作为模型外加差应力.图 1中两相矿物模型应力和位移边界条件如下：下边界固定，上边界施加恒定应力200 MPa，两个侧边应力均为100 MPa.矿物相块边界相互耦合，位移连续，没有空隙和叠加.矿物相的局部力学行为假设为各向同性.蠕变过程中模型保持恒温，不考虑摩擦热和剪切热的影响，也不考虑剪切过程中的力不稳定性.因为首先，模型中的矿物蠕变过程非常缓慢，近似为准静态过程；其次，我们目的是分析模型内部应力在岩石破裂或颗粒滑动等不稳定力学行为发生之前的阶段性发展.为确定合适的网格划分密度，将图 1中每块矿物相中的有限单元个数分别划分为：5×5，10×10，20×20和50×50进行试验.结果显示, 当每块正方形内网格密度达到10×10(即100)个单元后，应力的计算结果相差不大.因此，本文中的计算结果均为网格密度采用每块矿物相中10×10(即100)个单元得到的，模型中有限单元总数为40000个.

2.3 求解方法和过程简介

模型分析和求解是利用安装在国家超算天津中心上的Ansys软件来完成，我们采用plane182单元，该单元具有适用于二维实体的平面结构模型，为一阶单元，具有四个节点，节点在X和Y方向上为线性位移，具有弹性、塑性、黏塑性等功能.求解过程分两个载荷步进行：第一载荷步，加载时间为10^-6s，得到弹性的应力应变，作为第二个载荷步施加的初始应力应变；第二个载荷步，保持模型外加应力和位移不变且在第一个载荷步得到结果为初始条件下，分多个时间子步(共100步)求蠕变行为造成的应力应变变化，得到应力应变随时间的变化和在不同蠕变率矿物中的重新分配结果.我们选择隐式蠕变积分，没有屈服的判定，认为只要有应力存在，在经历足够长的时间，变形就会发生.同时，为了避免发生剪切锁定，对于plane182单元我们利用“选择性高斯积分(B-Bar方法)”，该积分对体积项积分时，在单元内部积分点为1个，对偏差项积分时单元内部积分点为4个.高斯积分点出现于刚度矩阵和单元内部，对于Ansys的有限元求解，高斯积分点上的应变是最准确的，所以程序默认采用高斯点上的应力应变解，然后外推到节点上求应变和应力.因此，Ansys计算蠕变的应力应变关系的流程如下：使用Von Mises方程计算蠕变应变率方程中使用的标量等效应力和等效应变，通过节点载荷向量[F^(m)]、总体刚度矩阵[K^(m)]、节点位移向量[u^(m)]的关系式(3)，在引入边界条件后，首先通过解方程组求节点的位移向量.对蠕变方程积分时，程序会按照蠕变方程求解等效应变增量，通过调整合适的时间步，使得每一步的应变保持在10^-5量级附近，通过牛顿-拉普森迭代线性逼近，并满足设置的力和位移收敛准则，计算积分点的蠕变应变，然后使用经过修改的总应变，得到等效总应变，从而计算出等效应力.等效应变增量基于程序设定的时间步和蠕变方程求出，然后从总应变中减去蠕变应变来得到弹性应变，接着利用Prandtl-Reuss方程，计算蠕变应变的分量(只有偏差分量).然后，根据物理方程(胡克定律，方程(1))计算出高斯积分点上的应力，最后通过形函数将高斯积分点上的应力外推到单元的节点，对于不同矿物相之间边界节点(由于相边界节点为公共节点，来自不同矿物相积分点求得的应力值不同)，则将该节点关联的所有单元计算得到的应力值进行平均.其中总体刚度矩阵方程为：

(3)

本研究主要从应力发展的角度分析蠕变差异对模型内部应力重新分配的影响，因此，我们只从计算结果中提取了节点的应力数值进行分析讨论.

3 结果 3.1 剪切应力在不同矿物模型中的发展

2.1节矿物模型中所述的三种模型：模型1(橄榄石-单斜辉石)，模型2(瓦兹利石-单斜辉石)，模型3(林伍德石-石榴石)，由于蠕变差异较大(图 2)，结果显示蠕变差异造成应力(即剪切应力，下文中的应力均指剪切应力)在模型内部重新分布(图 3)，应力主要集中在蠕变较慢的矿物上，而蠕变较快的矿物上应力明显较小(图 4)，该结果与前人结果(Madi et al., 2005; Bystricky et al., 2006)类似.作为对比，如果模型中两种矿物参数一致，则应力和剪切应力分布相对均匀，变化很小(如图 3d).模型2瓦兹利石-单斜辉石在相同时间内得到剪切应力比模型1橄榄石-单斜辉石的剪切应力大.如果应变率太小(图 3a)，接近于弹性，应力聚集效应不明显，或者需要很长的蠕变时间.如模型3林伍德石-石榴石由于两种矿物都蠕变较快，并且差异很大，虽然温度(600 ℃)比模型1、模型2的温度(800 ℃)低，但是剪切应力发展却快很多.我们的多次试验结果表明，只要保持图 1中两相的体积比例(60%~40%)不变，不同的矿物相随机分布，得到剪切应力的量级相似.

3.2 矿物蠕变差异造成的剪切应力对俯冲带应力的影响

虽然相对于俯冲带的复杂环境，本文的模型还比较简单，但是揭示了矿物蠕变差异在俯冲带剪切应力发展过程中的重要作用.我们模型的边界条件需要持续的应力来源，俯冲带产生的弹性应变率约为10^-15s^-1(许俊闪，2014).当蠕变应变与弹性应变可比时，才能对应力产生较大影响.因此，如果模型中两相矿物的应变率都远小于10^-15 s^-1(如本文中800 ℃下的橄榄石-单斜辉石模型)，则模型主要表现为弹性加载，蠕变效应可以忽略不计；如果模型中两相矿物的应变率都远大于10^-15 s^-1，接近于流变，则模型中的弹性应变能被消耗.只有当模型中一个矿物相的蠕变率大于10^-15 s^-1，而另一个小于10^-15 s^-1时(如本文中600 ℃下的林伍德石-石榴石模型)，即两相矿物模型的一种矿物表现为黏性，另一种表现为弹性，模型内部应力重新分配才能使得应力集中在蠕变较慢的矿物中，最终在硬的矿物相(表现为弹性的矿物相)中产生滑动或破裂.

图 5给出了几种矿物的蠕变率在俯冲带随深度的变化曲线(选取的是俯冲带不同深度的最低温度).随着深度的增加，俯冲板块中心温度上升较慢, 同时由于压力升高较快，蠕变活化焓增大，部分地抵消了温度升高而增大的蠕变率，所以总的来看，蠕变率随深度的变化都不太大.相对于模型2和模型3的矿物，模型1的橄榄石、单斜辉石两种矿物蠕变活化焓都较大，蠕变率都较慢.所以橄榄石-单斜辉石模型在相当长的时间内，对应力的重新分布影响较小.只在俯冲板块的上下边界处，由于温度较高，应力可能发展较快.而瓦兹利石和林伍德石在俯冲带温度应力条件下蠕变相对较快，使得应力更易集中在蠕变较慢的单斜辉石或石榴石中.

由于应变是个相对量，相同温度和应变率下由于模型尺度不同也会导致结果不同.图 1二维矿物模型中，矿物相块大小尺度为毫米(mm).如果尺度为米(m)，即不同矿物相分布的不均匀尺度是米，那么在此尺度上要达到图 3c(600 ℃)的应力结果，则需要约千年的时间.同样在米的尺度下，而如果温度增大为800 ℃，则时间又缩为1年左右.

4 讨论 4.1 俯冲带可能的应力不均匀区域

图 5显示，无论是橄榄石-单斜辉石，瓦兹利石-单斜辉石还是林伍德石-石榴石的矿物组合模型在各自稳定的压力和温度条件下，两相矿物蠕变率都相差数个量级.但在俯冲板块内部的橄榄石-单斜辉石相区，由于温度较低，最低处可能只有200~400 ℃(许俊闪，2014)，橄榄石和单斜辉石的蠕变差异对应力发展变化的影响很小，当橄榄石相变为瓦兹利石后，蠕变速率在一定温度范围内(600~900 ℃)较橄榄石快，应力应该较快地向较硬的单斜辉石集中.但是，橄榄石、单斜辉石的应变率都远小于板块俯冲可能产生的弹性应变率(约10^-15s^-1)，这说明蠕变的影响还远小于弹性应变的积累，板块变形还是以俯冲引起的弹性变形为主.不过，瓦兹利石是一种含水能力很高的矿物，含水的瓦兹利石的蠕变率远远大于无水瓦兹利石(Kohlstedt, 2006；Kawazoe et al., 2013)，在实际过渡地幔中，水的存在会加速应力重新分布的过程.而当瓦兹利石相变到林伍德石后，单斜辉石也变为石榴石，蠕变率都会很快.这时林伍德石的蠕变率远大于10^-15s^-1，而石榴石小于10^-15s^-1，使得应力很快地向石榴石积累.当温度大于1000 ℃后，即使是干的石榴石蠕变率也超过了10^-15s^-1，石榴石和林伍德石的较快蠕变率使得应力不但难以积累，还能消耗掉板块俯冲产生的弹性应变，这时变形以蠕变为主，不会再积累应力而发生地震.图 6给出了俯冲带可能发生剪切应力的温度区域分布.其中俯冲带温度模型基于Negredo等(2004)，其他参数见许俊闪(2014).图 6给出了受矿物蠕变差异影响较大的可能发生应力重新分布的温度范围和深度区域，对于520 km以下的林伍德石和石榴石区域约为600~1000 ℃，是俯冲带的核心区(红色区域，图 6).而对于瓦兹利石-单斜辉石相区约为800~1200 ℃，在俯冲带的边缘(绿色区域，图 6).如果考虑到瓦兹利石的含水可能性，将扩展到俯冲板块的核心区域，即图 6中的蓝色区域.对于橄榄石-单斜辉石相区，橙色为其受蠕变差异影响的应力发展区域.而板块中心的灰色区域是弹性控制的应变区域.因此，在上地幔深度，俯冲带的核心区由于温度较低相对稳定，仅在板块上下边缘的高温区应力发展较快，考虑到板块的厚度，高温区已接近边界，俯冲产生的应力可能难以维持.而在过渡地幔深度矿物发生相变后，由于瓦兹利石和林伍德石的应变率较快，在温度较低的俯冲板块核心区域，仍可发生蠕变差异造成的应力重新分布，使得应力急剧发展，对地震的孕育提供了应力积累的基础.尤其是林伍德-石榴石相区，应力发展要快几个量级，会对地震的孕育起到极大的促进作用，可能是地震数量在600 km深度附近增多的原因.而当温度太高时(>1000 ℃)，应变率太快，反而消耗了应力的积累，地震无法发生.另外，从深源双震带的复杂情况看，其可能与矿物相变关系密切(张克亮和魏东平，2008；Guest et al., 2004)，这与瓦兹利石-单斜辉石相区的应力发展区域(图 6，绿色区域)分布类似.

4.2 对反裂隙断层作用和剪切热失稳作用的补充

反裂隙断层作用由于有实验支持，受到了较大的认可，但一直未有来自地震观测的证据支持(Green and Houston, 1995; Koper and Wiens, 2000; Frohlich, 2006).我们的模型显示，即使没有亚稳态相存在，不同矿物之间的蠕变差异使得应力重新分配，也可引起剪切应力的聚集，甚至随着应力的继续发展，可能发生不稳定的剪切，形成微破裂或小尺度的滑动.这在没有亚稳态相的情况下，从理论上提供了剪切不稳定的来源.剪切热失稳机制虽然没有实验数据支持，但是最新的地震学观测证据表明，深震的震源机制与剪切热失稳机制相符(Chen and Wen, 2015).剪切热失稳的发生需具备两个条件：(1)变形速率快；(2)热反馈大，并对失稳起促进作用.而实验中林伍德石原位变形实验的结果Xu等(2003)确实显示了其较快的变形速率，远远大于Weetman理论模型和佩尔斯机制(许俊闪，2014)的计算结果.同时，由我们矿物蠕变差异对剪切应力发展起到的促进作用看，正是由于林伍德石较快的蠕变率加大了其与石榴石的蠕变差异，为600 km深度附近的剪切失稳提供了更大的可能性，加速了失稳的进程.

Liu和Zhang(2015)提出的壳幔层在高温高压下的体积收缩差异机制，从宏观尺度上解释了应力的来源和产生机理.该机制从相变产生的应力随深度变化的包络线角度类比解释俯冲带地震数量随深度的分布.由于俯冲带环境非常复杂，地震发生的数量还是应该以应力发展快慢来衡量，当然，我们的模型与其并不矛盾，他们的结果解决了俯冲带深部的应力来源问题，也为我们的模型提供了外部应力条件，使我们的模型在较小尺度上成为可能.由于受到较大的整体外部应力作用，矿物相之间的蠕变差异将会更明显，更有利于应力的重新分布，即应力向蠕变较慢的矿物聚集，从应力的角度为600 km附近地震数量的增多提供了更大的可能.

5 结论

本文从俯冲带主要矿物的蠕变差异角度分析了俯冲带地震可能产生的机制.提出了矿物蠕变差异会对俯冲带内部应力的发展起到重要作用的观点.为进一步理解深源地震机制提供一个新参考.无论从温度、压力等背景环境的角度看，还是从弹性变形、塑性变形等力学角度看，或者从俯冲带历史和矿物相变角度看，俯冲带环境都非常复杂，单一的机制很难从根本上解释深源地震的形成问题.因此，本文的矿物蠕变差异加速剪切应力发展的观点能为更好地理解俯冲带深源地震机制提供一种新参考.

致谢

感谢国家超级计算天津中心和天河一号对程序运行的支持.特别感谢审稿专家和美国洛斯阿拉莫斯国家实验室陈彧博士对本文提出的意见和建议.