引言

龙门山断裂带位于青藏高原东缘和四川盆地交汇处。在不到5年时间内，龙门山断裂带的中段和南段相继发生了2008年5月12日汶川M_S8.0地震和2013年4月20日芦山M_S7.0地震，给人民的生命财产造成了巨大损失。地震发生后，国内外学者从不同角度对这2次地震进行了深入研究。其中，库仑破裂应力作为一个重要的指标，在分析地震之间的相互关系、探讨震后断层应力水平的变化以及估算地震的发震周期方面具有独特的意义 (Harris，1998；Parsons等，2008；陈连旺等，2008；王辉等，2008；邵志刚等，2010；徐晶等，2013；李玉江等，2013a；2013b；2014；Wang等，2014；Jia等，2014)。因此，许多学者从库仑破裂应力角度研究地震的影响。Parsons等 (2008)采用地震波反演得到的均匀弹性有限断层模型，通过计算后发现汶川地震引起的同震应力变化使龙门山断裂带南段库仑破裂应力增加了0.1MPa，该断裂带正是芦山地震发生的位置。Toda等 (2008)和万永革等 (2009)采用弹性位错模型分别计算了汶川地震后周围断层上库仑破裂应力的变化，认为龙门山断裂北部和最南端、鲜水河断裂带南段、岷江断裂的地震危险性增加。在此基础上，单斌等 (2009)进一步考虑了地壳自重及分层结构。董培育等 (2013)计算了芦山地震后周边断层的库仑应力变化。单斌等 (2013)也采用同样的方法计算了汶川地震和芦山地震引起的同震和震后粘弹性松弛应力场变化，结果表明，汶川地震导致芦山地震震源处库仑破裂应力增加0.37-1.14×10^-2MPa，汶川地震可能促进了芦山地震的发生。Wang等 (2014)采用三维粘弹性分层位错模型进行分析也得到了相同的结论。缪淼等 (2013)采用有限断层模型计算了汶川地震的静态库仑应力变化对芦山地震的影响以及芦山地震对余震的影响，指出芦山地震是由汶川地震触发造成的。

上述工作使人们对汶川和芦山这2次地震相互关系的认识取得了重要进展。然而在计算中，他们采用的弹性位错模型或有限断层模型都属于简单的物理模型，不能反映不同区域的地形、地质构造、介质性质以及在横向和纵向上的非均匀性等。对此问题，数值模拟方法则具有天然的优势。首先，它能根据已有的观测资料和研究结果建立更符合实际情况的物理模型；其次，能针对所研究的问题，通过调整参数进行数值实验，从而分析不同参数的影响。因此，许多学者采用数值模拟方法，通过建立合适的有限元模型计算同震库仑破裂应力，分析地震对周边地区的影响 (Luo等，2010；Liu等，2014)。李玉江等 (2013a)充分考虑地形差异和粘弹性松弛等因素，并用接触摩擦方法处理活动断裂，研究了汶川地震对周围断层的影响。随后，他又采用相同的方法进一步分析了芦山地震对川滇地区主要断裂带的同震加载和卸载效应 (李玉江等, 2013b, 2014)，但他没有讨论这2次地震之间的关系。

因此，本文采用数值模拟方法，根据地质构造、深部反演结果以及GPS观测资料等，建立能反映地表起伏和岩石圈分层结构的龙门山地区三维粘弹性有限元模型。采用Maxwell粘弹性本构关系模拟岩石圈在漫长的地质演化过程中的流变效应以及地震发生瞬间的弹性效应。以GPS观测值为约束条件重建研究区现今构造应力场。采用准静态方法依次模拟了汶川地震和芦山地震的发生过程。通过分析同震库仑破裂应力与余震分布的关系，探讨了这2次地震主震对余震的触发作用以及汶川地震对芦山地震的影响。最后计算了这2次地震的共同作用对周边断层应力变化的影响。

1.   三维粘弹性有限元模型

1.1   模型的建立

根据龙门山地区的活动构造特征，本文选取101.5°-106.8°E，29°-33.7°N的矩形范围为研究区域，如图 1所示。三维模型包括了龙门山断裂带 (后山、中央和前山断裂)、岷江断裂、虎牙断裂、龙日坝断裂、鲜水河断裂带南段等主要活动断裂，如图 2所示。根据深部反演结果，青藏高原东缘和四川盆地在地壳厚度以及物性参数上存在巨大的差异，因此在横向上将模型分为东侧的青藏高原东北缘和西侧的四川盆地2个地质部分。三维有限元模型沿纵向 (Z轴) 从地表到地下200km共分为5层，分别是地表、上地壳、下地壳、岩石圈上地幔和软流圈上地幔，包含活动地块和主要断裂带2种地质单元。断层宽度大约取5km，从地表至地下20km，为简化起见，所有断层均为直立。由于研究区域不大，采用笛卡尔坐标系代替球坐标系进行计算，以东向为X轴正向，北向为Y轴正向，垂直向上为Z轴正向，采用Guess-Kruger投影算法将研究区的经纬度坐标投影为XYZ直角坐标。引入ETOP高程数据¹得到模型地表垂直方向坐标，根据Moho面深度值 (来自MapSIS软件提供的陈学波关于中国及外围地区莫霍面深度分布图) 得到地壳与地幔的分界面。为考虑地形附加的重力影响，在地表单元，通过设定密度值引入体力作用。使用六节点三棱柱单元对三维模型进行剖分，最终得到69003个节点，120210个单元。采用通用有限元计算程序ADINA作为计算平台。

图 1  研究区域 (根据邓起东等，2011)

Figure 1.  Research area of this study (after Deng et al., 2011)

下载: 全尺寸图片 幻灯片

图 2  龙门山地区三维有限元模型

Figure 2.  3-D finite element model of the Longmenshan area

下载: 全尺寸图片 幻灯片

1.2   本构关系和材料参数

在漫长的地质演化历史中，岩石圈介质表现为流变效应，可以用粘弹性介质来描述。本文采用Maxwell粘弹性本构模型，其长期表现为流变性质，瞬时表现为弹性性质，可以用来模拟岩石圈万年尺度的应力应变演化过程以及地震发生时的瞬时弹性效应，有本构关系为：

其中，K(t) 和G(t) 分别为Maxwell模型的体积模量和剪切模量，η为粘滞系数，K为弹性体积模量，G为弹性剪切模量。后2个参数可由线弹性本构模型的杨氏模量E和Possion比υ换算得到 (周光泉等，1996)。根据研究区深部反演获得的三维P波、S波速度结构 (吴建平等，2006；王椿镛等，2008) 可以计算得到杨氏模量E和泊松比υ。根据前人在该地区的研究成果 (石耀霖等，2008；杨强等，2010；曹建玲等，2009) 确定粘滞系数。考虑到断层带介质一般相对周边活动地块介质较软，将断层的杨氏模量取为其两侧介质杨氏模量平均值的1/3，取泊松比为比两侧介质的平均值高0.02(杨兴悦等，2013；李玉江等，2013a)。地表单元的密度取为2.7×103kg/m³。最后得到的研究区介质分层材料参数，如表 1所示。

表 1  研究区介质分层材料参数

Table 1.  Material parameters of the layered medium in the research area

分层	深度/km	四川盆地			青藏高原东缘
		E/1010·Pa	υ	η/1022·Pa·s	E/1010·Pa	υ	η/1022·Pa·s
地表	0-0.5	3.75	0.21	0.8	3.75	0.21	0.8
上地壳	0.5-20	7.2	0.246	9.96	7.39	0.25	10.3
下地壳	20-Moho面	12.5	0.253	0.925	11.9	0.254	2.02
岩石圈	Moho面-100	17.5	0.265	0.05	17.50	0.265	0.05
软流圈	100-200	17.5	0.265	0.5	17.50	0.265	0.5

下载: 导出CSV 

| 显示表格

1.3   模型加载条件和边界条件

为了进行计算，需要对模型施加加载条件和边界条件。根据地表GPS观测值 (张培震等，2008) 在模型4个侧边界通过插值得到位移加载条件，如图 3所示。由于GPS值在深部与地表的差异尚无定论，因此假设加载条件从地表到地下200km深度保持一致 (李玉江等，2013a)。考虑到GPS观测值代表地壳水平位移的年变化值，因此设置加载时间函数为每年增加1个载荷单位，这样就保证了时间间隔为1年的相邻两个状态，加载条件的增加量为地壳水平位移的年变化值，从而使加载条件在物理意义上是合理的。模型底边界条件为沿垂向 (Z) 固定，沿水平 (X，Y) 向自由；上表面为自由边界。

图 3  模型加载条件

Figure 3.  Loading condition of the model

下载: 全尺寸图片 幻灯片

1.4   模型检验

在以上介质参数、加载条件和边界条件下，对模型加载40万年，并使用ADINA自带的位移消除算法解决因施加重力导致的模型塌陷问题，最后得到稳定的应力状态，将其作为研究区现今构造背景应力场，在此基础上模拟地震的发生。为检验模型的可靠性，将模拟得到的地壳水平运动速率与GPS观测值进行比较，如图 4所示。对比模拟值与观测值可以看出，两者的方向和大小在大部分地区都非常接近，仅在鲜水河断裂带附近存在一些差异，这可能与鲜水河断裂带活动性较强有关，也不排除GPS观测值本身存在误差。定量计算得到模拟值与观测值方向夹角的平均值A为3.934°，长度相对差异的平均值L为0.142。因此，我们认为，模拟结果基本反映了研究区域的现今构造变形状态，将此状态下的应力场作为研究区现今构造背景应力场，即汶川和芦山地震发生的背景应力场状态。

图 4  地壳水平运动速率模拟值与观测值的对比

Figure 4.  Comparison of the simulated and the observed values of the crustal horizontal movement velocity

下载: 全尺寸图片 幻灯片

2.   库仑破裂应力累积速率

2.1   库仑破裂应力计算方法

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

其中，μ为摩擦系数，τ为主震产生的应力张量在断层面上的剪切应力大小，以沿断层滑动方向为正，σ_n为断层面上的正应力，以张应力为正，p为孔隙压力。通常更关心的是库仑破裂应力的变化 (Harris，1998)，当μ不随时间变化时，由式 (2) 可得库仑破裂应力的变化ΔF_CFS为：

孔隙压力变化控制着断层面上的有效正应力，对于各向同性均匀介质，孔隙压力对摩擦系数的影响可用等效摩擦系数μ'=μ(1-B) 表示，其中B为Skempton系数，范围为0-1(万永革等，2000；缪淼等，2013)，则式 (3) 变为：

等效摩擦系数μ'取值通常为0-0.75，平均值约0.4(Freed，2006)。数值实验表明，取不同μ'值计算得到的库仑破裂应力的空间分布相差不大，但应力的大小有一定的差别 (万永革等，2000)。本文计算中取μ'=0.4。一般认为，库仑破裂应力变化大于0.01MPa是触发后继地震的阈值 (Harris，1998)。

2.2   断裂带库仑破裂应力累积速率

将模拟得到的节点应力张量的6个分量投影到所选断裂带的断层面上，再根据式 (4) 计算断层上库仑应力的年变化，可得到该断裂带上库仑破裂应力的年累积速率 (陈连旺等，2001；李玉江等，2013a)。对研究区内主要断裂带均采用以上方法进行计算，最后得到龙门山及周边地区主要断裂带库仑破裂应力年累积速率分布，如图 5所示。由图可以看出，库仑破裂应力年积累速率最大的是鲜水河断裂带，为1.2-1.6×10^-3MPa·a^-1，表明鲜水河断裂带上的应力积累速率很快，发生较强地震的周期较短，这与鲜水河断裂带历史强震频发的规律基本一致。对于龙门山断裂带，中段和北东段累积速率很小，约为0.2-0.3×10^-3MPa·a^-1；南段的累积速率比北东段大，约为0.4-0.6×10^-3MPa·a^-1，因此本文认为龙门山断裂带南段的活动性要强于北东段，这与史翔等 (2009)和陈立春等 (2013)通过历史地震研究得到龙门山断裂带北东段在第四纪活动不明显、第四纪以来无强震记录、南段至少有3次6级以上强震记录的结论是相符的。

图 5  龙门山及周边地区主要断裂带上的库仑破裂应力年累积速率

Figure 5.  Cumulative Coulomb failure stress of year in the main fault zone in Longmenshan and its surrounding area

下载: 全尺寸图片 幻灯片

3.   计算结果

3.1   汶川地震和芦山地震的数值模拟

地震在一定背景应力场条件下发生，与漫长的地质演化历史相比，地震的发生是一个瞬态过程，在力学上可视为一个弹性过程。如果忽略地震中发震断层的动态破裂过程，而仅仅研究地震前后的静应力状态变化，可以采用准静态方法模拟地震发生。本文使用ADINA的重启动算法引入研究区现今构造背景应力场，采用降低发震断层介质参数的方法模拟实际地震的发生。根据弹性力学理论 (王敏中等，2011)，对于弹性材料，泊松比υ与体积模量K和剪切模量G的关系为：

一般认为介质的泊松υ比是常量，因此我们在模拟地震时，按相同比例同时降低K和G。为了使模拟地震的断层错动尺度与实际情况一致，选取发生地震时断层的水平滑动位移作为约束条件。关于地震模拟中发震断层的设置，对于汶川地震，以实际地表破裂为发震断层 (徐锡伟等，2008；2013a；邓起东等，2011)，即沿龙门山中央断裂的北川-映秀破裂 (长约240km)，沿前山断裂的汉旺-白鹿地表破裂 (长约72km)，沿小鱼洞的地表破裂 (长约6km)，共设置3条发震断层。而对于芦山地震，由于震后地震地质调查没有发现明显的地表破裂痕迹 (徐锡伟等，2013b)，且对于发震断层的确定尚未形成共识 (许冲等，2014；陈立春等，2013)，因此我们选择震中位置所在的双石-大川断裂作为芦山地震的发震断层，根据地震类型、破裂长度与震级的统计关系 (Wells等，1994) 计算得到断层破裂长约46.7km，并且在设置发震断层时，使破裂未到达地表。表 2给出了模拟2次地震所使用的计算参数。

表 2  汶川地震和芦山地震参数

Table 2.  Earthquake parameters of the Wenchuan earthquake and Lushan earthquake

事件	经度	纬度	走向	倾角	滑动角	最大水平滑动量/m	破裂长度/km	数据来源
汶川	103.4°E	31.0°N	231°	35°	138°	4.9	> 300	中国地震台网中心；Global CMT；徐锡伟等，2008
芦山	102.89°E	30.31°N	212°	44°	92°	1.3	46.7	USGS；曾祥方等，2013；徐锡伟，2013；张勇等，2013

下载: 导出CSV 

| 显示表格

3.2   主震对余震的影响

为研究主震对余震的影响，本文收集了汶川地震主震发生后至2008年7月8日共2553次2.5级以上余震，以及芦山地震后至2013年4月26日共445次2.5及以上余震的定位结果。分别计算了汶川和芦山2次主震的同震库仑破裂应力变化，并分析其与余震空间分布的关系。图 6是汶川地震的同震库仑破裂应力变化及余震的空间分布。结果表明，汶川地震有88.2%的余震发生在主震库仑破裂应力大于0.01MPa的加载区，说明主震触发了大部分余震。计算得到汶川主震震中位置的同震库仑破裂应力为0.442MPa。根据前文计算的龙门山断裂带中段的库仑破裂应力年积累速率为0.2-0.3×10^-3MPa·a^-1，如果按照该累积速率，并忽略其它影响龙门山断裂带库仑破裂应力累积的因素，则汶川地震的发震周期约为1473-2210年，这与张培震等 (2008)通过地震地质方法估算的复发周期为2233-5952年，以及Ran等 (2013)通过古地震研究得出的复发周期约为2865-2960年从量级上是一致的。图 7是芦山地震的同震库仑破裂应力变化及余震的空间分布。计算表明，芦山地震有80.3%的余震发生在主震库仑破裂应力的加载区，其中有75.42%的余震位于库仑破裂应力大于0.01MPa的区域，说明芦山地震的大部分余震也由主震触发。相比汶川地震中主震对余震的影响，芦山地震的余震被主震触发的比例较小，这可能是由于芦山地震的部分余震发生在与发震断层相交成Y字型的反冲断层上，该断层上的地震活动由主震触发 (房立华等，2013)，故这部分余震相当于是被主震间接而非直接触发。另外也可能与芦山地震的震级相对汶川地震较低、发震构造特殊以及破裂没有贯穿到地表有关。总体而言，从库仑破裂应力角度，汶川地震和芦山地震的大部分余震都是被其主震所触发的。

图 6  汶川地震的同震库仑破裂应力变化与余震的空间分布

Figure 6.  Coseismic Coulomb stress change of the Wenchuan earthquake and the spatial distribution of aftershock

红色实心圆圈表示汶川地震主震，白色圆圈表示余震

下载: 全尺寸图片 幻灯片

图 7  芦山地震的同震库仑破裂应力变化与余震的空间分布

Figure 7.  Coseismic Coulomb stress change of the Lushan earthquake and the spatial distribution of aftershock (Solid red star represents the main shock of Lushan earthquake, the white circle represents the aftershock)

红色实心五角星表示芦山地震主震，白色圆圈表示余震

下载: 全尺寸图片 幻灯片

3.3   汶川地震和芦山地震的关系

芦山地震发生在汶川地震之后不足5年的时间内，2次地震之间存在怎样的关系，学术界还有争议 (陈运泰等，2013；刘杰等，2013；徐锡伟等，2013b；陈立春等，2013；Jia等，2014)。为此，本文将从库仑破裂应力和余震分布角度进行一些探讨。图 8显示了汶川地震产生的库仑破裂应力变化在芦山地震破裂面上的投影。可以看出，芦山地震位于汶川地震同震库仑破裂应力增加的区域，增加量约为0.016MPa，这与其他学者采用有限断层模型或位错模型的计算结果一致 (Parsons等，2008；Toda等，2008；万永革等，2009；缪淼等，2013；Wang等，2014；Jia等，2014)。因此，从库仑破裂应力角度，汶川地震对芦山地震有触发作用。如果按照龙门山断裂带南段库仑破裂应力年累积速率为0.4-0.6×10^-3MPa·a^-1计算，汶川地震使芦山地震提前发生约27-40年，这与前人的研究结果是一致的 (缪淼等，2013；Wang等，2014)。

图 8  汶川地震的同震库仑破裂应力在芦山地震破裂面上的投影

Figure 8.  The projection of the coseismic Coulomb stress of the Wenchuan earthquake on the fracture plane of the Lushan earthquake

红色实心圆圈表示汶川地震，红色实心五角星表示芦山地震

下载: 全尺寸图片 幻灯片

为了进一步分析汶川地震主震对芦山地震余震的影响，将汶川地震和芦山地震的同震库仑破裂应力叠加再投影到芦山地震破裂面上，其结果与芦山地震余震的空间分布如图 6所示。计算表明，芦山地震的余震有85.8%发生在2次地震总的库仑破裂应力的加载区，其中有82.2%的余震位于库仑破裂应力大于0.01MPa的区域，比单独考虑芦山地震同震库仑破裂应力的情况分别增加了5.5%和6.78%，说明汶川地震对芦山地震余震也有触发作用，但影响非常有限，芦山地震的余震大部分还是被其自身主震触发的。

图 9  汶川地震和芦山地震共同产生的同震库仑破裂应力变化与余震的空间分布

Figure 9.  Coseismic Coulomb stress change of the Wenchuan earthquake and Lushan earthquake as well as the spatial distribution of aftershock

红色实心五角星表示芦山地震，白色圆圈表示余震

下载: 全尺寸图片 幻灯片

3.4   汶川地震和芦山地震对周边断层的影响

本文计算了汶川地震和芦山地震发生后产生的总库仑破裂应力对周边断层的影响，并以此来评估周边断层未来的地震危险性。表 3给出了研究区内主要活动断层的参数 (邓起东等，1994；李传友等，2004；Toda等，2008；万永革等，2009；单斌等，2009；缪淼等，2013) 以及最后计算得到的在该断层上的库仑破裂应力变化。图 10为汶川地震和芦山地震共同导致的龙门山地区主要活动断层的库仑破裂应力变化。可以看出，汶川和芦山地震后，龙门山断裂带库仑破裂应力变化较大，在北东段，大部分地区应力减小，仅在中央断裂部分应力增加，变化范围为-4.2-3.2MPa。在龙门山断裂带南段，后山断裂和中央断裂应力增加，前山断裂部分应力减小，总的库仑破裂应力变化范围为-2-1.8MPa。龙日坝断裂和岷江断裂整体呈应力减小趋势。虎牙断裂应力增加，最高值为0.08MPa。鲜水河断裂南段带应力增加约0.005-0.15MPa。在这些地区未来发生地震的可能性增大。

表 3  研究区内主要活动断层参数及汶川地震和芦山地震共同引起的库仑破裂应力变化

Table 3.  Parameters of the main active faults in research area and the Coulomb stress change caused by the Wenchuan earthquake and Lushan earthquake

断裂名称	走向/°	倾角/°	滑动角/°	库仑破裂应力变化/MPa
龙门山断裂带东北段	225	60	180	-4.2-3.2
龙门山断裂带南段	220	60	90	-2-1.8
鲜水河断裂带南段	142-159.5	90	0-45	0.005-0.15
龙日坝断裂	205-229	60	135	-0.067-0
岷江断裂	180	45-60	45-70	-2.96-0
虎牙断裂	150	75	45	0~0.08

下载: 导出CSV 

| 显示表格

图 10  汶川地震和芦山地震共同引起的龙门山地区主要活动断层的库仑破裂应力变化

Figure 10.  Coulomb stress change of the main active faults in Longmenshan and its surrounding area caused by the Wenchuan earthquake and Lushan earthquake

下载: 全尺寸图片 幻灯片

4.   结论与讨论

本文基于地质构造、深部反演结果以及GPS观测资料等，建立能反映地表起伏和岩石圈分层结构的龙门山地区三维粘弹性有限元模型。考虑地形附加重力影响重建了研究区现今构造背景应力场。在此基础上，采用数值模拟方法计算了汶川地震和芦山地震的同震库仑破裂应力。通过分析其与余震分布的关系，探讨了汶川和芦山2次主震对余震的触发作用以及汶川地震对芦山地震的影响。研究表明，汶川地震和芦山地震的余震大部分是由各自的主震触发的，汶川地震对芦山地震的余震有微弱的触发作用，比例约为6.78%。从库仑破裂应力角度，汶川地震对芦山地震有触发作用，并使后者提前约27-40年发生。汶川和芦山地震的发生使周边的鲜水河断裂带南段、虎牙断裂、龙门山断裂带部分区域库仑破裂应力增加，加大了这些断层发生地震的可能性。