据中国地震台网测定，2013年4月20日上午8时02分46秒在四川省雅安市芦山县龙门乡(103.0°E，30.3°N)发生了M_S7.0级地震.此次地震发生在2008年5月12日汶川M_S8.0级地震时没有发生破裂的龙门山断裂带西南段，以逆冲破裂为主.陈运泰等(2013)给出了此次地震的标量地震矩和最佳双力偶解的两个节面，并认为倾向北西的节面是芦山地震的发震断层面.由于芦山地震和汶川地震发生在龙门山断裂带上不同的段落，震中位置很近(87 km，房立华等，2013)，时间间隔仅5年，因此研究这两次地震之间的关系引起广泛关注(陈运泰等，2013；单斌等，2013；Liu et al.，2014，Li et al.，2014).同时，大量研究工作对芦山地震震源区(或龙门山断裂带)复杂的断层构造环境进行了分析探讨，如对震源区附近出露地表的断裂和隐伏断裂带活动特征及其对芦山地震发震破裂的控制作用的研究(徐锡伟等，2013；张岳桥等，2013；李传友等，2013；周荣军等，2013；李渝生等，2013；Chen et al.，2014；《芦山地震科学考察》编委会，2015)；青藏高原物质向东运动速率变化对加速芦山地震孕震过程的影响(赵静等，2013；刘峡等，2014；孟宪纲等，2014)；龙门山断裂带西部中上地壳低速层和滑脱面的存在有利于龙门山断裂带由北西向南东的推覆作用和逆冲作用，及其对芦山地震震源位置的控制作用(朱介寿，2008；Burchfiel et al.，2008；王椿镛等，2008；刘启元等，2009；李勇等，2013；詹艳等，2013； 滕吉文等，2014).上述研究多为在野外调查、地球物理探测和形变观测基础上的定性推测，对震源深浅部构造环境、青藏高原物质向东运动速率变化、低速层和滑脱面对芦山地震孕震过程影响的定量化数值模拟分析，将有助于进一步推进相关研究的分析深度，定量分析各类主要动力学因素所起作用的主次和影响程度.因此本文在上述研究基础上，建立数值模型对上述主要构造条件与芦山地震孕震过程影响进行数值模拟实验研究.

芦山地震的破裂面没有扩展到地表(徐锡伟等，2013；陈运泰等，2013).房立华等(2013)给出的主震和余震定位结果显示，芦山地震的主震震源深度为17.6 km，余震震源深度的优势分布范围为5~20 km，其发震断裂很可能是大川—双石断裂向东到大邑断裂之间存在的多条隐伏断裂.在横跨龙门山断裂带的剖面上芦山地震余震的分布呈“y”型(房立华等，2013；Fang et al.，2015；赵荣涛等，2015).Fang等(2015)和赵荣涛等(2015)依据余震分布图像推测芦山地震的主震破裂较为复杂，包括了两条倾向相对的破裂面，空间分布与余震分布范围基本一致.芦山地震后，中国地震局第二监测中心于2013年6—9月对穿过震区的测线进行了等精密水准复测(郝明等，2014)，并根据2010—2013年穿过震区的水准观测数据进行计算处理，得到了芦山地震同震垂向位移.本文试图以芦山地震同震位移的形变观测资料为约束，通过正演方法判断余震精确定位显示的“y”形分布是否可以确定芦山地震主要破裂面为“y”型双破裂面.

本文采用数值模拟的方法，以多种方法给出的孕震环境地震波速度结构和电磁结构的研究结果确定模型的弹性参数和流变参数；以余震精确定位、震源机制和震源破裂过程研究结果确定地震破裂面的几何形态和位置；以同震形变观测结果为约束检验模型中发震破裂面几何模型的合理性，对芦山地震的孕育条件及发震机理进行了模拟研究.数值模拟实验重点分析在区域高海拔地势引发的重力势能差异影响(祝爱玉等，2015)、上中地壳间低速层和滑脱面、震源区主要断层和汶川地震后青藏高原物质向东运动速率变化等因素对芦山地震震源位置和孕震机理的影响.同时也对芦山地震破裂面具有“y”型复杂破裂特征的可能性进行了验证.具体模型(图 1)的建立依据如下：

图 1

Fig. 1

图 1 二维数值模型示意图 模型为芦山地震震中位置(102.983°E，30.291°N)垂直剖面，模型总长为200 km，深度为100 km，并以高原和盆地的分界线为界限，将模型分为两大块，即青藏高原东缘和四川盆地.模型表面自由；底部垂向和水平向都固定；模型的东南部边界水平向固定，垂向自由；模型的西北部受东南向的挤压作用，上地壳加载速度为4 mm·a-1，模型下部分加载速率为2 mm·a-1；整个模型处于重力作用下. Fig. 1 The two-dimensional numerical model The model is the vertical profile of the epicenter of Lushan earthquake(102.983°E,30.291°N). The width of the model is 200 km,and the depth is 100 km. The dividing line between the plateau and the basin divides the model into two parts,i.e,the eastern margin of the Qinghai Tibet Plateau and the Sichuan basin. The surface of the model is free; As for the bottom of model,both of the vertical and horizontal directions are fixed; As for the southeast boundary,the horizontal direction is fixed,and vertical direction is free; As for northwest boundary,the loading rate of the upper crust is 4 mm·a-1 and the rest part of northwest boundary is 2 mm·a-1.

(1)主要断层    根据已有研究结果，本文重点讨论大川—双石断裂和大邑断裂对芦山地震孕震环境和孕震过程的控制作用，在数值模型中采用相对周围材料较低杨氏模量和黏滞系数的软弱夹层方式表示大双—双石断裂和大邑断裂.并设定大川—双石断裂和隐伏的大邑断裂具有上陡下缓倾角特征，大川—双石断裂的倾角在30°～65°范围内(周荣军等，2006)，大邑断裂在近地表的倾角在40°～70°范围内(董绍鹏等，2008)，且断裂带底部连接近水平向的滑脱面；断层宽度为5 km.

(2)发震断层    在讨论余震“y”型分布特征是否为芦山地震的主震破裂面时，以软弱夹层方式表示“y”型地震破裂面，材料参数的选定以陈运泰等(2013)和郝明等(2014)给出的发震断裂和地表同震位错结果为约束，通过试错法确定.“y”型地震破裂面的空间分布形态主要参考Fang等(2015)给出的和本文二维剖面空间位置一致的跨芦山震源位置剖面余震精定位分布图像，破裂面深度为5～20 km.芦山地震发震破裂面的厚度主要依据Fang等(2015)给出的余震分布范围，水平宽度为5 km.

(3)滑脱面和低速层    在川西高原的松洛甘孜褶皱带内，在20 km深处的中地壳之上存在厚约3~5 km的低速高导层，可能代表了一个深部滑脱层(邓起东等，1994).大邑断裂为龙门山前缘扩展变形带的前缘断裂，断层面倾向北西，向下呈铲状，并汇交于滑脱面.该滑脱面很可能就是芦山地震的震源层(李勇等，2013).所以本文在深度20 km处设置了一条接近水平向的滑脱面，该滑脱面和大双—双石断裂，大邑断裂底部相连接.地震波速度结构层析成像(王椿镛等，2003)及大地电磁测深探测(孙洁等，2003)结果发现川西高原的上地壳存在低速层.本文按照层析成像的研究结果(王椿镛等，2003)在滑脱面以上15~20 km建立了一层低速带，并采用低的黏滞系数和低的弹性模量来表示该低速层.在本文开展模拟实验时，将滑脱面和低速层对研究区域应力分布的影响一并进行讨论.

(4)分层结构模型    分层结构模型主要包括青藏高原东部和四川盆地，其深度方向分为上地壳、中地壳、下地壳及上地幔，模型深度为100 km.其中模型表面能反映地形，且高原地区的地形高度约为4 km，四川盆地的地形高度约为400 m.上中地壳底部参考了层析成像(王椿镛等，2003)的研究结果.Moho面深度参照跨龙门山断裂带的爆破地震剖面观测结果(朱介寿，2008).

模型中密度根据最新的芦山地震科考金川—洪雅剖面(其穿过主震区，北西向测线，长约240 km)，得到模型上地壳密度为2.7 g·cm^-3，中地壳密度为2.85 g·cm^-3，下地壳密度为2.93 g·cm^-3，上地幔密度为3.35 g·cm^-3.王椿镛等(2003)给出了研究区域的P波速度，青藏高原东缘的上、中、下地壳及上地幔P波速度分别为6、6.25、6.8 km·s^-1、7.8 km·s^-1；由于四川盆地地块比较坚硬，其上、中、下地壳及上地幔P波速度分别为6、6.35、6.8 km·s^-1、8.0 km·s^-1；低速层的P波速度为5.8 km·s^-1.根据王椿镛等(2008)的研究结果，青藏高原东缘的上、中、下地壳及上地幔S波速度分别为3.5、3.8、4 km·s^-1、4.2 km·s^-1；四川盆地的上、中、下地壳及上地幔S波速度分别为3.5、3.8、4 km·s^-1、4.3 km·s^-1；低速层S波速度为3.5 km·s^-1.根据P波速度、S波速度及密度，采用公式(1)和(2)计算，分别得到泊松比和杨氏模量(见表 1).

表 1(Table 1)

表 1 模型材料参数表 Table 1 Material parameters of the model 地层 杨氏模量 /GPa 泊松比 黏滞系数 /Pa·s 青藏高原 东部 上地壳 中地壳 下地壳 上地幔 82.165 99.327 115.84 153.15 0.2421 0.2068 0.2354 0.2958 5×1023 1×1023 5×1022 1×1022 四川盆地 上地壳 中地壳 下地壳 上地幔 82.165 100.50 115.84 160.66 0.2421 0.2210 0.2354 0.2969 5×1023 2×1023 1×1023 5×1022 低速层 80.283 0.2137 1×1021 大川—双石断裂和大邑断裂(考虑为活动状态时) 50 0.2421 5×1022 芦山地震发震断层 (地震发生后状态) 80 0.2700 5×1022

表 1 模型材料参数表 Table 1 Material parameters of the model

根据该区域相对偏高的平均大地热流值，可推测川滇菱形地块的下地壳介质强度较低、相对较软，而巴颜喀拉地块的东南边的四川盆地的地壳介质则相对较硬(石耀霖和曹建玲，2008).有限元模型介质的黏滞系数列于表 1.其中确定低速层的具体黏滞系数还缺少足够的探测和实验依据(Clark and Royden，2000).赵国泽等(2008)通过对青藏高原东边缘及其附近地区石棉—乐山剖面大地电磁资料的研究，对青藏高原东边缘中地壳的低阻层黏滞系数做了估计，认为其黏滞系数在10¹⁸~10²¹ Pa·s 之间.根据赵国泽等得到的中地壳低阻层深度判断，其层位对应了王椿镛等给出的上地壳和中地壳之间的地震波速度低速层.有限元模型介质的黏滞系数见表 1.

模型中考虑大川—双石断裂和大邑断裂的活动状态时，其材料参数见表 1；假设断裂带处于不活动状态，材料参数与青藏高原东部上地壳材料参数相同.

青藏高原物质挤出运动被坚硬的四川盆地所阻挡，在龙门山地区形成了一系列的推覆构造，导致芦山地震以逆冲为主的特点.根据江在森等(2009)的结果，龙门山断裂带西侧(青藏高原内部)500 km的跨度上(沿着本文二维模型走向)地壳缩短量为3~5 mm·a^-1.根据Shen等(2005)相对华南地块的GPS形变观测结果显示，青藏高原东缘的平均速率为5 mm·a^-1，并随着与四川盆地距离的减小而逐渐减小.本模型边界距离龙门山断裂带约80 km，因此采用边界位移速率条件时，设定在四川盆地内部的模型南东端部边界水平向固定，在青藏高原内部的北西端施加水平向位移边界条件.上地壳位移速率边界条件选定为4 mm·a^-1，中下地壳和上地幔选定为2 mm·a^-1.为了体现底部地幔物质的拖曳力，同时减少边界效应对计算结果的影响，本文将模型深度设置为100 km，并设定模型底部水平向和垂向为刚性固定.

根据模型应变和应力变化趋势探讨各种动力学因素对芦山地震孕育过程的影响，初始条件在模拟计算过程中不能被忽视.在模拟初始应力时，模型考虑了重力作用，以表示高海拔的地形对区域构造应力的影响，以及由于长时间的岩石圈的流变作用形成的 “静岩压力”.形成“静岩压力”状态的计算总时间为40万年.同时本文采用最新的初始位移消除算法(ADINA R & D，Inc.，2010)将长时间的重力作用所造成的网格严重变形恢复到模型的初始网格状态，同时保证地应力几乎不变(Zhu et al.，2016).

众多研究结果均将青藏高原高海拔地势蓄积的重力位能作为影响高原物质向东水平运动的重要动力学因素之一(张健和石耀霖，2002；熊熊和滕吉文，2002；Zhu et al.，2016)，为体现重力作用造成的不同地块重力势能差别和重力作用在形成现今构造应力场过程中的重要作用，本模型加入重力作用，并且将东南边界和西北边界设置为水平向固定，垂直向自由.本模型的计算采用ADINA有限元软件，经过40万年的计算，模型几乎趋于静岩压力状态(尹祥础，1985)，然后采用位移消除算法将网格复原到初始网格并保持应力几乎不变，以此应力作为模型初始应力场(Zhu et al.，2016).然后施加上述2.3节中的位移边界条件，采用时间步长为500年，计算了100个时间步，得到5万年加载后的研究区稳定地壳构造应变场和应力场.本文对该区域的垂向沉降速率及水平向形变等计算结果的可靠性进行检验(见图 2，模拟结果考虑了低速带和滑脱面、大川—双石断裂、大邑断裂带的存在).

图 2

Fig. 2

图 2 研究区域的主要形变约束 (a)水平向形变速率；(b)垂向形变速率. Fig. 2 The main constrains of the model (a)The rate of the horizontal deformation;(b)The rate of the vertical deformation.

模拟计算结果如图 2a显示，模型西北部边界到四川盆地，其水平向速率逐渐减小.水平向速率的计算结果和GPS形变观测结果几乎一致；同时本文也考虑了地表隆升速率测量结果对模型地表隆升速率模拟结果的约束(图 2b).根据区域精密水准测量得到青藏高原东缘平均隆升速率达6 mm·a^-1，而四川盆地平均下降速率为2 mm·a^-1(孟宪纲等，2014).由于控制地表隆升速率的因素非常复杂(李祖宁等，2002)，考虑本模型在地壳分层、地幔对流对地壳影响、剥蚀作用(Kirby et al.，2002)、重力均衡状态等方面(李勇等，2005)均做了一定程度的简化，因此本文以模拟结果在量级上和自西向东逐渐减小等指标上定性满足测量条件为约束.

本文在建立模型部分介绍了龙门山断裂带西侧青藏高原可能存在低黏滞度的低阻层，及其对龙门山断裂带推覆构造形成影响的研究结果.本文试图通过比较存在低速层和滑脱面与否的两个模型最大剪切应力计算结果，初步讨论低速层和滑脱面对龙门山断裂带逆冲推覆构造运动和芦山地震震源区断裂带发生地震剪切破裂滑动力学条件的影响作用(见图 3).由于龙门山断裂带发育有大量断裂，并具有不同倾向，在这部分的模拟实验中，本文没有采用针对特定断层面投影的库伦破裂应力分布图像.分析图 3a和图 3b中最大剪切应力分布特征的差异可以看出，青藏高原物质相对四川盆地向东挤压运动时，无论是否存在低速层和滑脱面，均在龙门山中央断裂带及其西侧的上地壳形成了高剪切应力区，有利于形成龙门山逆冲推覆构造带；同时也在芦山地震震源位置形成了高剪切应力环境，有利于断裂带的地震剪切破裂滑动.

图 3

Fig. 3

图 3 模型左侧边界上地壳水平位移速率4 mm·a-1时，考虑低速带和滑脱面存在(a)和不考虑低速带和滑脱面存在(b)的最大剪切应力模拟结果(为突出本文探讨的芦山地震震源部分，本图仅给出模型中间部位地表至50 km深、100 km宽范围的应力状态) Fig. 3 The simulated maximum shear stress of considering(a)and without considering(b)the low velocity zone and detachment surface,when the loading rate of the upper crust of the northwest boundary is 4 mm·a-1(In order to highlight the source of Lushan earthquake,this paper presents the stress state of the surface to 50 km depth and a wide range of 100 km in the middle part of the model)

图 4给出了低速层和滑脱面存在与否时，模型地壳的垂直运动速率模拟结果，比较图 4a和图 4b可以看出，存在低速层和滑脱面时，龙门山地区的地表抬升速率模拟结果较没有低速层和滑脱面的模拟结果更接近现今地形特征，也就是说，低速层和滑脱面的存在更有利于形成龙门山现今地形.

图 4

Fig. 4

图 4 模型左侧边界上地壳水平位移速率4 mm·a-1时，考虑低速带和滑脱面存在(a)和不考虑低速带和滑脱面存在(b)的模型垂直位移速率模拟结果 Fig. 4 The simulated rate of vertical deformation considering(a)and without considering(b)the low velocity zone and detachment surface,when the loading rate of the upper crust of the northwest boundary is 4 mm·a-1

2013年4月20日芦山M_S7.0级地震发生在龙门山断裂带西南段，主震震源和余震活动范围基本限于大川—双石断裂到大邑断裂之间(Fang et al.，2015)，这两条断裂带与芦山地震的关系引起广泛关注和讨论.本文在考虑长期稳定边界位移的条件下，通过分别模拟计算两条断裂带处于无震滑动状态和不活动状态的长期稳定构造条件下芦山地震震源区的库伦破裂应力状态年变化率，讨论了这两条断裂带对芦山地震孕震过程的长期控制作用(图 5).在模拟各种断裂带活动方式时均考虑了滑脱面和低速层的存在.库伦破裂应力计算所需要的断层面几何参数以陈运泰等(2013)给出的芦山地震主震破裂面几何参数为依据.计算时断层面摩擦系数为0.6(Byerlee，1978).

图 5

Fig. 5

图 5 模型左侧边界上地壳水平位移速率4 mm·a-1时，考虑大川—双石断裂和大邑断裂活动与否的库伦破裂应力年变化率分布图(各情况均考虑了滑脱面和低速层的存在) (a)大川—双石断裂和大邑断裂都活动；(b)大川—双石断裂和大邑断裂都不活动；(c)大邑断裂活动而大川—双石断裂不活动；(d)大川—双石断裂活动而大邑断裂不活动. Fig. 5 The annual change rate of the Coulomb failure stress,when the loading rate of the upper crust of the northwest boundary is 4 mm·a-1 (a)Both of Dachuan-Shuangshi fault and Dayi fault are active;(b)Neither of Dachuan-Shuangshi fault and Dayi fault is active;(c)Dayi fault is active and Dachuan-Shuangshi fault isn′t active;(d)Dachuan-Shuangshi fault is active and Dayi fault isn′t active.

比较依据芦山地震主震破裂面倾向和倾角计算得到的库伦破裂应力变化图 5a和图 5b可以看出，在考虑青藏高原物质东流形成长期稳定挤压作用的边界条件控制下，受青藏高原高海拔地形蓄积的重力位能影响、上地壳底部滑脱面和低速层的存在影响，无论大川—双石断裂和大邑断裂是否处于活动状态，芦山地震震源区均为模型中库伦破裂应力增加幅度最大的部位.同时也注意到，滑脱面端部与大邑断裂交汇部位是另一个库伦破裂应力增加幅度高值区.两个断裂带的活动对断裂带临近区域的库伦破裂应力变化有一定的吸收和隔离作用，但对整个龙门山断裂带库伦破裂应力变化的影响并不明显.比较图 5b和图 5c可以看出，芦山地震震源东南侧的大邑断裂活动与否，对龙门山断裂带东部和四川盆地山前地带上地壳倾向北西的低角度逆冲断裂带库伦破裂应力变化图像格局有一定影响；比较图 5b和图 5d可以看出，芦山地震震源西北侧的大川—双石断裂活动与否，对芦山地震震源区上部的库伦破裂应力变化图像格局有一定影响.

汶川地震后，由于龙门山断裂带中北段由闭锁状态转为震后蠕滑状态，龙门山断裂带西侧青藏高原地块相对华南地块垂直于龙门山断裂带走向的水平运动速率由汶川地震前的4 mm·a^-1增加至8 mm·a^-1(赵静等，2013).这种横跨龙门山断裂带两侧的地壳缩短速率增加，势必造成龙门山断裂带南段压缩速率变大.为分析这种横跨龙门山断裂带的压缩速率突然变化与芦山地震孕震环境的关系，本文分别考虑了大川—双石断裂和大邑断裂处于活动状态和不活动状态的情况，模拟计算了芦山地震主破裂断层面(陈运泰等，2013)投影的库伦破裂应力1年尺度变化率图像(图 6).模拟计算考虑了低速层和滑脱面的存在，模型左侧边界上地壳边界的水平位移速率由4 mm·a^-1增加至8 mm·a^-1.计算时断层面摩擦系数为0.6.

图 6

Fig. 6

图 6 考虑模型左侧边界上地壳水平位移速率由4 mm·a-1增至8 mm·a-1的芦山地震主破裂面库伦破裂应力年尺度变化率模拟结果 (a)大川—双石断裂和大邑断裂处于活动状态；(b)大川—双石断裂和大邑断裂处于不活动状态(图中五角星为房立华等(2013))给出的芦山地震主震震源位置. Fig. 6 The annual change rate of the Coulomb failure stress in the main rupture surface of the Lushan earthquake,when the loading rate of the upper crust of the northwest boundary is changed from 4 mm·a-1 to 8 mm·a-1 (a)Both of Dachuan-Shuangshi fault and Dayi fault are active; (b)Neither of Dachuan-Shuangshi fault and Dayi fault is active.

从图 6可以看出，作为龙门山断裂带两侧地块相对挤压速率增加的响应，模拟结果图 6a和6b的芦山地震主破裂面投影库伦破裂应力变化率最高值的部位均处于模型中滑脱面与大邑断裂交汇位置.Fang等(2015)给出的芦山地震震源也落在这个库伦破裂应力变化率最大区域范围内.而Fang等(2015)给出的芦山地震余震则主要分布在这个高值区上部上地壳的库伦破裂应力增加率较大区域.

比较图 6a和6b可以看出，芦山地震震源东南侧的大邑断裂活动与否，对龙门山断裂带东部和四川盆地山前地带上地壳倾向北西的低角度逆冲断裂带库伦破裂应力变化图像格局有一定影响；大川—双石断裂活动与否，对芦山地震震源区上部的库伦破裂应力变化图像格局有一定影响.这和上一节中长期尺度库伦破裂应力变化的趋势是一致的.

以芦山地震后几个月内基于水准数据得到的芦山7.0级地震震间和同震位移场特征为地表同震位移约束(郝明等，2014)，以陈运泰等(2013)给出的芦山地震标量地震矩和震源破裂过程反演结果作为地震主震破裂面平均滑动量的约束，模拟了不同形态主震破裂面条件下的模型垂直位移分布特征(图 7)，对余震精确定位结果给出的余震分布“y”型分布特征是否可以作为推测主震具有“y” 型复杂破裂特征的依据进行了初步探讨.

图 7

Fig. 7

图 7 假定芦山地震破裂形态的不同空间分布条件下，模型及地表垂向位移分布特征的模拟结果 (a)“y”型破裂且破裂未达到地表条件下垂直位移模拟结果；(b)同(a)相同条件下地表垂直位移的模拟结果和观测结果比较；(c)“/”型破裂且破裂未达到地表条件下垂直位移模拟结果；(d)同(c)相同条件下地表垂直位移的模拟结果和观测结果比较；(e)“\”型破裂且破裂未达到地表条件下垂直位移模拟结果；(f)同(e)相同条件下地表垂直位移的模拟结果和观测结果比较；(g)“y”型破裂且破裂达到地表条件下垂直位移模拟结果；(h)同(g)相同条件下地表垂直位移的模拟结果和观测结果比较. Fig. 7 The vertical displacement and surface displacement at different spatial distributions of Lushan earthquake rupture (a)The simulated vertical displacement of the “y” type rupture and the rupture didn′t reach at the surface;(b)The surface vertical displacement of the same condition of(a);(c)The simulated vertical displacement of the“/”type rupture and the rupture didn′t reach at the surface;(d)The surface vertical displacement of the same condition of(c);(e)The simulated vertical displacement of the “\”type rupture and the rupture didn′t reach at the surface;(f)The surface vertical displacement of the same condition of(e);(g)The simulated vertical displacement of the “y” type rupture and the rupture reached at the surface;(h)The surface vertical displacement of the same condition of(g).

考虑主震破裂为复杂“y”型双破裂面模型的垂直位移模拟结果显示(图 7a和7b)，在破裂面平均滑动量与陈运泰等(2013)给出的芦山地震主震破裂面平均滑动量基本一致条件下，地表最大位移约为0.2 m，地表以下位移上升区域主要集中在“y”型区域的上岔口，并延续到了大川—双石断裂，模拟得到的地表垂直位移图像在数值和形态上和同震地表垂直位移的形变观测结果(郝明等，2014)较为一致；图 7c和图 7d给出了主震破裂为倾向北西的简单“/”型破裂面模型垂直位移模拟结果：与“y”型破裂面模拟结果比较，地表最大位移模拟结果约为0.13 m，且影响范围明显减小；图 7e和7f给出了主震破裂为倾向南东的简单“\”型破裂面模型垂直位移模拟结果：地表垂直位移最大值约为0.08 m，最大位移集中在大川—双石断裂的中部，垂直位移的幅度和形态相对“y”型和“\”型破裂面模型的模拟结果显著偏离同震形变观测结果.图 7g和7h给出了假定破裂面出露地表时的“y”型双破裂面模型垂直位移模拟结果：地表最大位移约为0.3 m，相对破裂面未出露地表的模型(图 7a和7b)，其地表最大位移加大了50％，并且影响范围明显增加，偏离同震形变观测结果.

(1)主要断裂活动与否对芦山地震孕震过程控制作用

对中上地壳内已经存在的断裂带，可以采用Byerlee摩擦定律作为判断其进入地震破裂状态的初步近似准则(Byerlee，1978).在模拟实验中，选择了依据芦山地震主震破裂面倾向和倾角计算得到的库伦破裂应力变化图像作为讨论的参考依据(图 5).本文在讨论相邻断层对芦山地震孕震条件的长期演化过程影响时，所建立的模型仅考虑了大川—双石断裂和大邑断裂，实际上该区域构造活动非常复杂，存在多条断裂带(李勇等，2013)，这些断裂带的存在均会影响芦山地震震源区的应力演化，更全面考虑这些断裂带影响还需要进一步研究.

(2)青藏高原相对四川盆地位移速率变化与芦山地震的关系

本文在分析青藏高原相对四川盆地位移速率变化与芦山地震关系时，位移速率变化参考了赵静等(2013)的研究结果，设定边界位移速率由4 mm·a^-1增加至8 mm·a^-1.讨论主要依据模型边界位移速率变化引起的库伦破裂应力变化是否有利于芦山地震发震断层的破裂，计算时考虑的受影响破裂面依据陈运泰等(2013)给出的倾向北西的芦山地震主震破裂面和滑动方向.图 6a和6b显示大川—双石断裂和大邑断裂之间的上地壳和中地壳过渡深度(15~25 km范围)是投影在芦山地震主破裂面的库伦破裂应力增加值最大部位，这实际上也是芦山地震的发震位置(Fang et al.，2015).模拟结果显示，汶川地震后龙门山断裂带中北段由闭锁状态转为震后蠕滑状态，导致青藏高原物质相对四川盆地挤压速率的增加，确实形成了有利于芦山地震发生的构造应力变化条件.从库伦破裂应力年变化速率高值区的分布看，低速层和滑脱面的存在是控制库伦破裂应力显著增加区正好处于芦山地震震源区的主要控制因素，比较图 6a和6b可以看出，大川—双石断裂和大邑断裂的活动与否，并没有影响芦山地震震源区的库伦破裂应力增加.

本文分析模拟结果时也注意到，边界条件突然变化条件下计算得到的库伦破裂应力年变化量最大值达到2 MPa的量级.这个数值已经接近地震过程中应力降的量级，远远大于已有多数研究结果给出的判断应力触发的地震库伦破裂应力计算值(平均为0.01 MPa，万永革等，2002)，也远远大于单斌等(2013)计算得到的汶川地震破裂导致的芦山地震发震断层面上的库伦应力变化值.分析造成模拟结果偏高的原因，可能是：本文模拟实验直接在西侧边界采用GPS观测给出的青藏高原物质相对四川盆地运动速率由4 mm·a^-1增加至8 mm·a^-1的结果作为边界速率变化条件，建立模型时根据保证计算易于收敛和计算量的需求忽略了大川—双石断裂以西的龙门山断裂带其他多条断裂带和褶皱带对压缩应变的吸收作用，这可能造成芦山地震震源区域压缩应变偏大，进而导致库伦破裂应力变化值偏大的计算结果.

(3)同震形变和余震定位结果与主震破裂过程的关系

本文以震后数个月内开展的水准观测给出的地表同震垂直位移观测结果为约束，探讨了“y”型余震分布是否反映出芦山地震的主震破裂面为两条倾向相对的破裂面.从观测资料的可靠性分析，本文依据的水准观测资料是地震后数个月后获取的，郝明等(2015)在分析了不同影响因素对同震形变的影响量级后认为，虽然无法严格地将震后滑移和同震位移区分开来，但是利用震后数个月内观测资料进行同震变形计算研究，可暂不考虑震后变形的影响.

关于“y”型构造是否反映了主震破裂空间分布，Fang等(2015)给出了芦山地震后180 h内不同时间段余震的跨震源区北西-南东走向剖面空间分布图像，可以看出震后3~12 h、12~24 h、24~48 h、48~96 h、96~180 h的余震空间分布图像均显示出“y”型的余震分布特征.这种短时间内给出的余震分布特征，应该与主震的破裂面形态分布非常接近.Scholz(2002)也认为目前通常的研究还是利用余震分布确定主震破裂范围.Fang等(2015)利用芦山地震后一年的余震观测资料开展的余震精确定位研究，和赵荣涛等(2015)等利用芦山地震后一年多的微、小余震开展的精确定位研究，更加清晰地显示出余震的“y”型分布特征.赵荣涛等(2015)认为，余震的“y”型分布说明芦山地震可能是两条余震带所对应的两条断裂同时活动所产生.实际上，正如房立华等(2013)所分析的，在龙门山山前地区，断裂带的逆冲滑动受阻而反向逆冲形成反冲断层，共同构成“y”型构造组合的例子并不局限于此次芦山地震的震源破裂.

综合上述的模拟实验结果及讨论，本文得到结论如下.

(1)在龙门山断裂带西侧的青藏高原存在的低速带和滑脱面，在长期地质尺度对龙门山断裂带地区的地表隆升有明显的控制作用，同时也对在芦山地震震源区的较大范围形成高剪切应力区有长期尺度的控制作用.

(2)从长期尺度的应力环境演化过程看，大邑断裂处于无震蠕滑状态或不活动，对龙门山断裂带东部和四川盆地山前地带上地壳的最大剪切应力分布格局有一定影响；大川—双石断裂和大邑断裂处于无震蠕滑状态或不活动状态，模拟结果均显示芦山地震震源区处于高剪切应力环境，两个断裂带活动与否对芦山地震震源位置长期的应力积累过程影响不大.

(3)汶川地震后短期内龙门山断裂带中北段由闭锁状态转为震后蠕滑状态，青藏高原东部物质相对四川盆地向东的运动速率增加，同时由于龙门山断裂带西侧上中地壳低速层和滑脱面的存在，导致龙门山断裂带的芦山地震震中附近区域在汶川地震后承受了更大的水平方向应变压缩作用.在低速层和滑脱面的控制下，这种水平向挤压应力的增强，促进了芦山地震的孕育和发生过程.而大川—双石断裂和芦山断裂的活动与否，对芦山地震的孕育和发生影响并不大.

(4)在本文开展的数值模拟实验中，根据余震“y”型分布给出两个倾向相对的破裂面建立模型所得到的地表垂直位移计算结果与形变观测结果符合较好.支持了房立华等(2013)提出的芦山地震主震可能有两个破裂面的推测.至于两个破裂面是同时发生破裂，或具有一定的先后顺序(赵荣涛等，2015)，本模型利用目前的观测研究结果所进行的模拟实验验证还不能给出确定的结论.