2001年11 月14 日中国青藏高原北部的东昆仑断裂带上发生了昆仑山8.1 级强烈地震.中国地震台网确定的震中位置为:35.95°N，90.54°E；震源深度10km.根据多位学者研究结果^[1~4]，其震源破裂过程虽略有差异，但均认为这次地震的破裂以左旋走滑为主，断层走向接近东西走向，断层面几乎直立，破裂扩展主体方向为自西向东发展.这个地震是自1951年11月18日西藏当雄8.0级大震以后半个多世纪以来中国大陆发生的最大地震.六年半之后，2008年5月12 日在其东侧位于青藏高原东缘的龙门山断裂带上又发生了汶川8.0级逆冲型强烈破坏性地震.研究这个地震在青藏高原构造块体系统^[5]大的背景中破裂的运动学和动力学过程及对汶川大震孕育发生的影响，对青藏高原构造块体系统稳定性认识，特别是对汶川大震发生影响的认识有着重要意义.

自1976年唐山7.8级大震发生后，我国大陆地区经历了12 年7 级大震的平静，1988 年进入7 级以上大震活跃期，发生了11 次7 级以上大震，这些大震均发生在青藏高原地区各构造块体的边界断层上(见图 1)，表明我国大陆地区这个阶段在青藏高原构造块体系统中块体边界断层的某些部位应变能积累已到相当水平，达到或接近破裂强度.同时也似乎表明，青藏高原地区的大震活动与构造块体系统中块体的运动变形的相互影响有一定关系.特别是从1996年至今，分别在巴颜喀拉构造块体周边边界断层上共发生2次八级大震及4 次七级以上大震.大震活动从未有过的频繁，显示了这些大震活动与该构造块体运动变形有密切关系.在汶川大震发生前，分别在巴颜喀拉块体的西、南、北边界断层上均发生了大震，惟独东边界断层---龙门山断裂带未发生大震.而2001年发生在巴颜喀拉块体北边界断层(东昆仑断裂带)上的昆仑山8.1 级大震，其破裂过程是左旋走滑错动，总体破裂方向自西往东.近东西走向的昆仑山大震发震断层南侧向东的猛烈错动必然引起巴颜喀拉块体进一步东扩，使得龙门山断裂带的应变能进一步加速积累.昆仑山大震的发生引起该构造块体的运动变形有可能对块体东部边界断层---汶川大震发震断层的提前失稳有一定影响，应该进行深入研究.

图 1

Fig. 1

图 1 1988年以来我国大陆地区6.9级以上大震及构造块体分布图 Fig. 1 Distribution map of larger than Ms6.9 earthquakes and tectonic blocks since 1988

非连续变形分析方法是一个分析块体系统中块体相互作用的有效工具^[6]，它可以计算块体系统在动力及静力作用下块体的滑移及不连续变形.对于各个块体允许有滑移，压性或剪切变形，对于整个块体系统，允许块体边界断层有滑动.至今，已有一些学者用这个方法进行多项研究，并取得重要成果^[7~10].汶川大震后我们用三维非连续变形与连续变形相结合的分析方法研究了汶川大震的孕震机理，取得了一些有意义的结果^[11].这个方法适合用来模拟大震的破裂过程和大震活动与构造块体的相互作用的关系.

本文使用三维流变DDA+FEM 方法，考虑青藏高原大的力学与构造背景及通过用GPS 资料和震源机制的约束条件，数值模拟昆仑山大震破裂过程及其对青藏高原及邻近地区构造块体系统稳定性的影响，研究该地区构造块体系统此次大震活动与各构造块体相互作用的关系，特别是昆仑山大震对汶川地震发震断层稳定性的影响.

我们编制了三维DDA+FEM 程序^[10]，关键之处在于将块体之间的相互作用转变为它们之间的接触问题^[12]，同时假设位移是很小的，块体之间的接触点是已知的.

设块体系统包含了许多块体，每个块体又被划分成一些有限单元，有限单元为六面体等参元.每个单元有8个节点，坐标为x_i(i=1，2…，8)，节点位移为u_i.(i=1，2…，8).它们的形函数用N_i(ξ，η，ζ)(i=1，2，…，8)表示，ξ，η，ζ 是局部坐标^[13].

有限单元的位移模式为:

(1)

其中，u，v，w代表单元内任意一点的三个位移分量，I是3×3的单位矩阵，e表示转置.

单元应变ε 与节点位移的关系为:

(2)

单元应力σ 由应力应变关系得到:

(3)

则最小位能原理的泛函总位能的表达式为:

(4)

对于离散模型，系统的位能是各单元位能的和，即

(5)

假设块体A表面上的点P与块体B表面的点Q接触.x_P，u_P 和x_Q，u_Q 代表接触点的坐标和位移.x_Q，u_Q 与点Q所属于的有限单元8个节点的坐标和位移有如下关系:

(6)

(7)

点P和点Q在接触面上局部坐标的位移分别由u^A和u^B表示.当接触状态是粘结时，接触力F^A由下式表示

(8)

其中λ 为虚拟弹簧系数，通常取得很大，以保证接触面不会嵌入.根据面力化为等效结点力的公式，可将接触面上n个接触点对有限单元结点的贡献结合到有限元计算格式中的刚度矩阵中去.

当满足莫尔-库仑准则，即处在摩擦滑动接触状态时，接触力F^A表示为:

(9)

(10)

其中f为摩擦系数，位移下标x，y表示在局部坐标系下平行接触面的位移，n表示局部坐标下垂直接触面的法向方向.

对于接触问题，我们将利用罚函数法将附加约束条件引入泛函以进行求解.此时，泛函表示为

(11)

其中Π_CP是用罚函数引入定解条件的附加泛函.对于黏结接触状态，

(12)

对于有摩擦的滑动状态

(13)

罚函数法的有限元求解方程如下:

(14)

其中^tK_L 是通常有限元的刚度矩阵，K_ca是接触面对刚度矩阵的贡献，^t+ΔtQ_c 是接触面滑移引起的附加载荷力，^tF是通常的外载荷等效结点力.

如果考虑岩石圈运动变形是一个长期过程，且假定应变率在一定时间段为常速率，则流变过程可采用类似解弹性力学问题的方法.流变本构关系取为下列形式:

(15)

这里ε_c 是流变应变张量，σ 是应力张量，而

(16)

其中η 是黏性系数，μ 是泊松比.

对于计算地震后应力的蠕变松弛过程，采用初应变方法.当时间步长较小时，可进行一次迭代计算.

本文取包括青藏高原及其东侧的四川盆地和鄂尔多斯块体作为研究区，见图 2a.地表面构造块体划分参考张培震等^[5]对我国大陆地区构造块体划分结果，将研究区划分为:塔里木、柴达木、阿拉善、祁连、巴颜喀拉、羌塘、拉萨、川滇、滇西南、鄂尔多斯、四川盆地及华南等十二个块体.块体边界断层包括阿尔金断裂带，祁连山断裂带、西秦岭断裂带，东昆仑断裂带、西昆仑断裂带、六盘山断裂带、贺兰山断裂带、玛尼-玉树带、喀拉昆仑断裂带，喜玛拉雅断裂带、龙门山断裂带、鲜水河断裂带、三江断裂带、红河断裂带、安宁河-小江断裂带、澜沧江断裂带、秦岭大别山断裂带、汾渭断裂带和阴山断裂带等近二十条断裂带.

图 2

Fig. 2

图 2 (a)构造块体系统模型及昆仑山和汶川大震发震断层分布；(b)龙门山断裂带地区地壳分层垂直剖面分布示意图 Fig. 2 (a) Tectonic block system model and distribution of ^eismogenic faults of Kunlun^han earthquake and Wenchuan earthquake； (b) Schematic diagram of distribution of vertical section of the crust on Longmenshan fault zone

对上述区域垂直深度120km 的空间划分为5层，在青藏高原东缘以西地区划分上地壳层20km, 中地壳层20km, 下地壳层30km 及上地幔.在青藏高原东缘及以东地区，考虑了龙门山断裂带东西两侧地势及地壳厚度和分层的明显变化^[14~17].地势从四川盆地海拔500 m 向西到龙门山地区陡升至4000m.地壳厚度在四川盆地为45km 向西逐渐增厚，在龙门山断裂带东侧地区为60km, 并逐渐向西增厚，其西侧的青藏高原地区地壳厚度为70km.在龙门山断裂带地区上地壳层及中地壳层各为15km, 下地壳层为30km.由于资料限制，为简化起见，除龙门山断裂带以外，其他边界断层均考虑为直立.龙门山断裂带的倾角在中地壳层为40°，在上地壳层为60°，见图 2b.为了提高块体内应力应变计算的精度，在每个块体内划分一些有限元网格.有限单元为六面体.整个块体系统共划分有4236个结点，3205个单元.底面垂直方向约束，水平方向可移动.在块体内，计算速度场时采用流变本构关系，假定应变率在一定时间段为常速率，则流变过程可采用类似解弹性力学问题的方法.由此计算的应力场由边界速度确定，而重力作用按弹性半空间模型计算的静岩应力简化处理，当作初始应力加入模型计算的应力中去.

考虑到青藏高原及其邻近地区P、S波速度结构的不均匀，参考大量青藏高原地区地震波三维速度结构的研究结果^[18~23]，本文对不同块体各层采用不同的物性参数，见表 1.其中拉萨块体上地壳层杨氏模量相对其他块体较小，而塔里木块体上地壳的杨氏模量相对较大，考虑到下地壳层处于温度较高的状态，速度低且高导层发育，介质呈较强的黏塑性，青藏高原整个下地壳层均考虑为软弱层，特别是羌塘及川滇两块体的下地壳层更软.对青藏高原东缘及以东地区，考虑了四川盆地及鄂尔多斯块体的地壳介质强度大，而龙门山断裂带西侧的松潘、甘孜一带地壳层的介质强度相对弱，特别是中下地壳层、在鄂尔多斯块体南侧与四川盆地北部之间存在一个介质强度较弱的地带；考虑到前人的波速反演弹性性质结果^[24]，本文各块体所取的泊松比，除了上地幔部分，取值均比文献[11]参数表中的值小了0.05，以便更符合实际情况.计算表明，两者会有一定差别，但位移场和应力场的基本特征基本一致.

表 1(Table 1)

表 1 计算模型中材料力学参数 Table 1 The mechanical parameters of materials in the model

表 1 计算模型中材料力学参数 Table 1 The mechanical parameters of materials in the model

为了研究青藏高原及邻区构造块体边界断层上的危险度分布及模拟昆仑山大震的破裂过程，我们要首先满足模拟计算的上地壳层速度场和应力场与青藏高原及邻近地区GPS资料获得的速度场^[25]及用地震资料获得的该区构造应力场^[26]一致.本文采用位移边界条件，地表速度值参考张培震GPS资料得到的相应地区地表速度场(参考框架为稳定的欧亚板块)结果^[25]，通过尝试，我们在一组边界条件下得到了如图 3a所示的地表速度场和部分点模拟计算的震源机制分布图(图 3b).模拟震源机制通过计算获得的主应力投影到乌尔夫网来表示.计算的地表速度场的大小从喜马拉雅断层带向北到阿尔金断裂带及祁连断裂带逐渐由大变小.速度方向从南边的喜马拉雅断裂带西南部的北北东方向随着向北推移逐渐向东偏转，在羌塘块体，巴彦喀拉块体及柴达木块体的东部位移方向接近向东的方向.并从龙门山断裂的西部向南到川滇块体位移方向逐渐再向南偏转，由南东方向逐渐转为南南东向.整个图像与由GPS资料获得的速度场的分布图像相似.计算的模拟震源机制图表明青藏高原周边为逆断层型震源机制，高原内部特别是喀拉昆仑断裂带主要以拉张性质的正断层震源机制为主，但在东昆仑断裂及玛尼-玉树带存在部分左旋走滑型的震源机制，与由地震波资料得到的震源机制解图及近年来发生大震的机制结果^[26]很类似.

图 3

Fig. 3

图 3 (a)计算得到的地表速度场分布；(b)计算得到的上地壳层模拟震源机制分布 Fig. 3 (a) Distribution of surficial velocity calculated in study zone；(b) Distribution of focal mechanism calculated from stresses tensor in upper crust

本文在青藏高原及邻近地区三维构造块体相互制约的环境中用三维DDA+FEM 方法数值模拟昆仑山大震的破裂过程.

参考许力生、陈运泰等^[1]对这个地震破裂过程的研究成果，本文考虑昆仑山大震的发震断层模型分布在上地壳层中接近东西走向的东昆仑断裂带，由三个相连接子事件组成，长度由西往东分别为110、170和180km(见图 2a).本文在上述计算初始速度场及应力场与GPS 资料得到的相应地区速度场及用大量震源机制结果获得的构造应力场结果大致一致的基础上，在初始应力场上附加一个扰动，即降低发震断层的摩擦系数，使得在上地壳层中的昆仑山大震发震断层处满足莫尔-库仑准则发生失稳，产生滑动，从而发生昆仑山大震破裂过程.

莫尔-库仑准则为:当τ_s ≥f·σ_n 时断层面发生滑动，其中τ_s和σ_n为断层面上的剪切应力和法向应力，f为摩擦系数.在模拟昆仑山地震发生过程时，发震断层的摩擦系数从0.4 降为0.1，后者使得发震断层满足莫尔-库仑准则，从而该处断层产生错动.

图 4a是模拟计算得到的昆仑山大震发生前后邻近地区主压应力变化分布与2004 年许力生和陈运泰给出的该震震源机制对比图^[1].图中明显显示在北西-南东方位上主压应力增加，与震源机制解中P波初动向上的区域一致；在北东-南西方位上主压应力减少与震源机制解中P 波初动向下的区域一致.由于发震断层方位为近东西走向，表明计算得到的昆仑山发震断层发生了左旋走滑错动，地震释放的主压应力轴方向为北东向.与图中震源机制解的结果基本一致.图 4b是大震发生前后最大剪应力变化等值线与徐锡伟等2002 年给出的昆仑山地震地表破裂带同震位移分布对比结果^[27].由图可见等值线图由三个大小不等并排的等值线图像组成，总体分布范围近东西向，与徐锡伟给出的近东西走向地震破裂带上三个峰值同震位移分布基本对应.图中还显示:发震断层南北两侧的等值线分布具有不对称特点，南侧最大剪应力变化梯度明显大于北侧.图 4c是发震断层垂直面上位错等值线图，图中显示大震破裂主要发生在30多公里深度范围内，与2004年周云好等用地震波资料反演得到大震发生在约30km 深度范围内的结果基本一致^[3].图 4d是计算获得的大震引起垂直位移变化三维分布结果，图中明显显示垂直位移变化在发震断层西段大体表现为南盘上升，而东段表现为北盘上升，与万永革等2004年用GPS 资料反演的结果基本一致^[28].本文计算获得的发震断层最大水平错动为4.5m, 与周云好等反演得到的3.6m 比较接近.计算得到的最大应力降为18MPa, 与周云好等反演得到的18MPa一致.

图 4

Fig. 4

图 4 昆仑山大震破裂过程模拟结果对比 (a)主压应力变化分布与震震源机制[1](a1)对比图；(b)最大剪应力场变化等值线与大震地表破裂带同震位移分布[27](b1)对比图，色标为主应力值；(c)发震断层垂直面上位错等值线(数字单位m)图；(d)大震引起垂直位移(色标数字单位m)变化三维分布图. Fig. 4 Contrast diagram of emulation results of the rupture process of the Kunlunshan earthquake (a) Contrast diagram of the principal stress change distribution and the focal mechanism [1](a1)；(b) Contrast diagram of maximum stress field change contour and coseismic displacement distribution of surfical rupture zone [11](b1) ； (c) Contour map of dislocation on the vertical plane of seismogenic fault； (d) 3D distribution of vertical displacement caused by the earthquake.

为了探讨昆仑山大震发生引起的各构造块体边界断层稳定性的影响，本文进一步计算昆仑山地震发生对各构造块体边界断层上库仑破裂应力的变化，为了直观和便于分析，我们用不同颜色表示不同变化范围.红色，黄色及深蓝色分别表示不同范围的库仑破裂应力增加，其中红色为明显增加.浅蓝色及绿色表示库仑应力减小的不同范围.这里所说的库仑破裂应力是地震学界通常用来描述地震破坏状态的一个量.它的定义如下:

(17)

其中τ 是断层面上的剪应力，σ_n断层面上的法向应力，μ 为摩擦系数，本文取为0.4.正的Δτ 和Δσ_n 使得断层接近破裂，负的Δτ 和Δσ_n使得断层远离破裂.因此，库仑破裂应力变化可以定量表示断层接近破裂的程度.

图 5(a, b, c)分别是计算得到的昆仑山大震发生引起青藏高原及邻近地区构造块体边界断层上、中和下地壳层边界断层上库伦破裂应力变化.图中显示昆仑山大震的发生使得其发震断层两端靠近发震断层某些地段库仑应力变化为明显增加区(红色区)，表明发震断层两个端部应力进一步集中.图中还显示，北侧的阿尔金断裂带中段，南侧的玛尼-玉树断裂带中段，喀拉昆仑断裂带中偏东部分地段及喜马拉雅断裂带偏东部分地段均有较明显库仑应力增加.南侧各边界断层上出现的红色区随着该带与东昆仑断裂带距离的增加出现在该带的位置也越往东移动，长度也越小.同时图中还显示在大震东侧龙门山断裂带的汶川大震发震断层上，库仑破裂应力增加0.016 MPa(上地壳层).我们以前的工作^[11]表明，龙门山断裂带属于地震危险地区.虽然汶川大震库仑破裂应力增加幅度相对红色区略小，但仍然有一定影响.看来，昆仑山大震的发生对汶川大震发震断层的孕震与失稳起到促进的作用.

图 5

Fig. 5

图 5 昆仑山大震发生对各构造块体边界断层库仑破裂应力影响 (a)上地壳层；(b)中地壳层. Fig. 5 Distribution of Coulomb failure stress change on the boundary faults of Qingzang and Chuandian block system after Kunlunshan earthquake (a) Upper crust； (b) Middle crust.

我们知道，岩石在应力加载后要经历三个阶段:首先应变速率递减；随后应变率保持常量；最后应变率加速增加.计算的地震后引起库仑破裂应力随时间的变化，反映蠕变的初始和稳态两阶段的过程，我们计算中所采用的黏度值实际上是这两过程的平均值，因而要小于计算稳态位移场时所取的黏度值.如果取黏度为1.58×10¹⁵Pa· s, 所计算得到的昆仑山大震发生导致的汶川大震发震断层中4 个点(见图 2)上地壳和中地壳的库仑破裂应力变化随时间变化如图 6和图 7所示.由图可看到:由于本文计算模型的复杂性，昆仑山大震发生引起地壳的黏弹性松弛及流变导致汶川大震发震断层库仑破裂应力的变化也比较复杂.表现在最初2~3年中变化呈现起伏，各点变化总趋势也略有不同.图中也明显显示，昆仑山大震发生后，从第2 或第3 年开始，1，3 和4点库仑破裂应力(上地壳层)有较明显小幅度增加，分别至第3年或第5 年达最高值，其后随时间变化不大.其中第3点(汶川大震初始破裂点)在第5年时库仑破裂应力为最高点，其后随时间变化不大.考虑到汶川大震是在昆仑山大震之后6 年半发生的，看来是处于变化最高点后1 年半的状态下发生的.在中地壳层，1，2和4点在经过初期起伏后，在第2年至第3年时段有明显小幅增加，然后趋于稳定.我们还注意到参数选择对结果有直接的影响.黏度量值选取与变化最长时间密切相关.

图 6

Fig. 6

图 6 昆仑山大震引起汶川大震断层处库仑破裂应力随时间的变化(上地壳) (a, b, c)和(d)分别对应1，2，3，4号点(见图 2). Fig. 6 Coulomb failure stress change with time at the fault of Wenchuan earthquake induced by Kunlunshan earthquake (upper crust) (a, b, c) and (d) corresponding to number 1，2，3，and 4 site (shown as Fig.2).

图 7

Fig. 7

图 7 昆仑山大震引起汶川大震断层处库仑破裂应力随时间的变化(中地壳) (a, b, c)和(d)分别对应1，2，3，4号点(见图 2). Fig. 7 Coulomb failure stress change with time at the fault of Wenchuan earthquake induced by Kunlunshan earthquake (middle crust) (a, b, c) and (d) corresponding to number 1，2，3，and 4 site (shown as Fig.2).

本文尝试在几大板块特别是印度板块的强烈俯冲及青藏高原和邻近地区构造块体的相互制约的构造大环境中，用三维流变非连续变形与连续变形相结合的方法，数值模拟计算昆仑山大震的破裂过程及对青藏高原构造块体系统稳定性的影响，取得了一些有意义的结果.

(1) 计算得到青藏及其邻近地区表层速度场分布图像与张培震等^[25]用GPS资料获得的该区速度场总体分布图像基本相似.计算得到这个地区表层震源机制的分布图像与许忠淮用大量地震震源机制结果得到的分布图像及近年来发生大震的震源机制也基本一致.尽管计算结果与GPS资料仍未达到完全一致，但图像的趋势基本相近.我们的结果仍然为昆仑山地震模拟结果接近实际情况打下基础.使得计算得到的昆仑山大震引起青藏高原及邻近地区构造块体边界断层上库仑破裂应力的变化结果具有一定的可靠性.

(2) 由大震释放的主压应力场图像最大剪应力变化等值线图、大震前后位移变化矢量图、大震发震断层垂直面上位错等值线图及大震引起垂直位移变化三维图分别与大震的震源机制，地表破裂带同震位移分布，GPS同震位移图及地震波反演和GPS反演的结果总体上均比较相近，计算获得的最大水平错距和最大应力降与用地震波资料反演的结果也基本接近.

(3) 昆仑山大震的发生引起青藏高原及邻近地区构造块体边界断层上库仑破裂应力的变化，除了使其发震断层两端靠近发震断层某些地段库仑破裂应力明显增加，应力在端部进一步集中外，北侧的阿尔金断裂带中段，南侧的玛尼-玉树断裂带中段，喀拉昆仑断裂带中偏东部分地段及喜马拉雅断裂带东部部分地段均有较明显库仑应力增加.大震南侧各边界断层上出现的红色区随着该带与东昆仑断裂带距离的增加出现在该带的位置也越往东移动，长度也越小，显示了该震左旋走滑错动的影响.同时图中还显示在大震东侧龙门山断裂带中的汶川大震发震断层库仑破裂应力增加0.016 MPa(上地壳层).虽然增加幅度相对前者略小，但这个地段属于同时具有失稳危险度较高^[11]及引起正的库仑破裂应力变化两个特点的地区.看来，对于失稳危险度较高，应力积累已接近破裂临界状态的边界断层，任何小的库仑破裂应力增加都有可能导致大震的发生.但是如果边界断层未接近临界状态，前面地震产生的应力变化就不可能触发后续地震的发生.昆仑山大震也引起其他个别地段库仑破裂应力增加较大，但至今没发生大震，可能其中有这个原因.类似于沈正康等^[29]的估算方法，如果假设断层每年错动2mm, 假设在100km 的空间范围内产生形变，地壳剪切模量取为4.4 ×10¹¹Pa, 则每年会产生10^-8量级的应变，相当于8800Pa的应力积累，这样0.016 MPa的库仑破裂应力变化相当于约2年的应力积累.

(4) 由昆仑山大震发生引起地壳的黏弹性松弛及流变导致汶川大震发震断层库仑破裂应力变化的计算结果表明，由于计算模型的复杂性，各点变化趋势也比较复杂.但图中仍可看到汶川大震是在昆仑山大震引起同震库仑破裂应力变化基础上，再经历小的起伏，并达到最高值后，趋于稳定的阶段发生的.值得指出的是，黏性参数的选择仍然是一个需要进一步探讨的问题.

(5) 昆仑山大震破裂过程是在近东西走向的发震断层上的左旋走滑错动，最大水平错距达4.5m, 以及大震的总体破裂方向自西向东发展，均引起发震断层南侧巴彦喀拉块体的向东运动和一定规模的变形.图 4b显示，发震断层南侧最大剪应力变化梯度明显大于北侧，间接地显示南侧巴彦喀拉块体剪应变变化增强明显大于北侧.巴彦喀拉块体的向东运动及变形的增强使得其东边界断层中低倾角的汶川发震断层库仑破裂应力增大，应变能积累增强，对处于危险度较高状态的汶川大震发震断层发生逆冲型失稳起了促进作用.说明某个构造块体边界断层上的失稳发生大震，不仅在该发震断层上产生应力降，释放应变能，同时由于各块体的相互制约，有可能向该块体其他边界断层及其周围的块体传递应力和引起运动变形，使其他边界断层上应力状态及应变能积累改变，加速或延缓边界断层上失稳.

(6) 2010年4月14日发生的玉树M_S=7.3大震发生在巴彦喀拉块体的南边界断层上.玉树大震前我们计算得到的昆仑山大震引起青藏高原及邻近地区构造块体边界断层上库仑破裂应力的变化图中，在距玉树大震地区相近的地段，也是库仑破裂应力增加比较大的地段，计算结果似乎有所反映，当然这方面情况还需要后续更多研究来证明.

(7) 本文采用的三维流变非连续变形与连续变形相结合的力学模型，相对以往的数值模拟研究模型更接近青藏高原地区的孕震环境.要说明的是，本文采用的岩石圈结构及物性选取比较简单，特别是龙门山断裂带以西地区.块体边界断层除汶川大震发震断层外，均为直立.参数的选取也只能参考现有的地震波三维速度结构的结果.这些简化处理，对模拟结果会有影响.相信随着今后资料的进一步积累，这项数值模拟研究将会更加完善.

对沈正康研究员的帮助和有益的讨论，以及审稿人的修改建议表示衷心感谢.