2017, Vol. 60
2. 中国地震局地球物理研究所, 北京 100081
2. Institute of Geophysics, China Earthquake Administration, Beijing 100081, China
库仑应力模型在研究余震触发及主震对后续强震的触发方面发挥着非常重要的作用，国内外很多科学工作者研究了地震触发问题，并取得了引人瞩目的成果 (King et al., 1994; Harris, 1998; Lin and Stein, 2004; Parsons et al., 2008; 石耀霖和曹建玲, 2010; Wan and Shen, 2010; 缪淼和朱守彪, 2012; Zhu and Miao, 2015).例如：King等 (1994)计算了1992年美国Landers地震 (MW7.3) 造成破裂面附近最优方向上的库仑应力变化，并考察了该地震对周边断层造成的应力变化，发现余震广泛分布于库仑应力增加0.01 MPa以上的区域，而库仑应力降低的区域余震活动较少.Deng和Sykes (1997)计算了1812—1995年间发生在南加州的中强震所产生的库仑应力变化总和，发现95%的6.0级以上余震均发生在加载区.尤其是Stein等 (1997)计算了1939—1992年发生在土耳其North Anatolian断裂带10个6.7级以上地震的库仑破裂应力变化，发现其中90%的地震是被先前地震所触发，并且成功地预测了1999年Izmit地区强震的发生.同样，Parsons等 (2008)计算了2008年汶川地震产生的库仑应力变化，发现雅安地区为库仑应力的增加区域，即该区域是地震危险区.果然，在2013年雅安地区发生了芦山MS7.0级地震，这进一步说明了库仑触发模型的有效性和可预测性 (缪淼和朱守彪, 2013).
然而，前人在计算静态库仑应力变化研究地震触发时，多数基于Okada (1985, 1992) 的解析解.我们知道，Okada模型是通过Skempton系数B来体现应力对流体的影响 (Skempton, 1954)，而没有考虑流体对固体骨架的作用.但实际上，流体与固体之间是完全耦合的，流体对于余震的触发起着非常重要的作用 (Nur and Booker, 1972).例如：Bosl和Nur (2002)研究了1992年美国Landers后孔隙压上升造成的库仑应力变化，与余震的时空分布极为吻合.He和Peltzer (2010)使用孔隙弹性介质计算了2008年西藏尼玛MS6.4地震的静态库仑应力变化，同震孔隙压变化使库仑应力上升了0.1 MPa，触发了后续的MS5.9和MS5.4地震.
本文正是基于上述问题，试图利用有限元计算方法，考虑流-固之间的完全耦合效应，针对三种类型的地震 (正断型、逆冲型和走滑型)，分别计算同震静态库仑应力变化；然后与利用传统Okada方法计算得到的结果进行比较，找出它们之间的差别，以期完善库仑应力模型.
2 流-固耦合作用的物理方程地下流体广泛地存在于地壳以及地震震源体中 (Zhao et al., 1996, 2000)，尽管流体与固体之间的相互作用非常复杂，但使用孔隙弹性理论可以很好的计算分析流体与固体骨架之间的力学耦合作用 (Kennedy et al., 1997; Coussy, 2004).
Biot (1941, 1956) 提出使用流体饱和多孔介质来分析流体的瞬态响应.各向同性流体饱和多孔弹性介质的本构方程为 (Rice and Cleary, 1976)：
|
(1) |
|
(2) |
其中σij为饱和孔隙材料的应变分量，m-m0是孔隙材料中单位体积内流体的质量变化，ρ0是孔隙流体的密度，p和σ分别是孔隙压力和平均应力偏量，B是Skempton系数，ν和νu分别是排水和不排水时的泊松比，G为剪切模量，δij为Kronecker符号.
当介质的应力变化极为迅速时，孔隙中的流体仍然处于静止状态，并未发生流动 (Bosl and Nur, 2002)，此时的孔隙弹性介质被认为是不排水的 (undrained state)，p=-Bσ(假定 (2) 式中m-m0=0)，即孔隙压力的变化与地震产生的平均应力相关.Skempton系数B是一个经验常数，用来量化多孔介质骨架传递到孔隙流体上的压应力.这种情况仅适用于断层滑动，流体还没有明显的扩散时.地震之后，孔隙流体将由地震产生的压力较高的压缩区向压力较低的膨胀区域传播.孔隙压力变化与体应变εkk相关 (Roeloffs, 1996)，公式为：
|
(3) |
在各向同性的多孔弹性介质中，如果断层滑动仅有剪应力或剪应变产生，那么孔隙压力并不会有任何变化.假定承压潜水面产生一个体应变变化Δεkk，则水头变化与体应变的比值为：
|
(4) |
若将应变作为一个阶跃函数加入模型中，由方程可知潜水面的水头会出现瞬时反应.然而，如果考虑偏应力的影响，则孔隙压力与骨架应力之间的耦合方程为 (Henkel and Wade, 1966)：
|
(5) |
其中A是由实验得到的系数，对于黏性土壤其范围一般由-0.5至1.5(Skempton, 1954; Wang, 1993)，对于灰岩一般A值为0.22(Hart and Wang, 1995).当A < 1/3时，断层滑动产生的孔隙压力变化小于使用平均应力计算的结果；当A=1/3时，孔隙压力会与偏应力解耦.当介质由不排水变为排水时，可知：
|
(6) |
其中K为排水时岩石的体积模量，Ks为固体骨架的体积模量，Ku为不排水的体积模量.
3 有限元模型及计算结果本研究利用有限元软件ABAQUS/Implicit进行模拟计算，该软件中的接触单元技术可以较好的模拟断层两盘的滑动过程；同时，对于固体骨架与饱和孔隙流体在排水情况下的孔隙弹性耦合问题也有较强的求解能力 (Schoenball et al., 2010).模型中的材料属性参考前人的做法 (Zhang, 2005; Zhou and Burbey, 2014) 来选取，三种类型断层地震模型的物质属性保持一致，假定地下流体为水，具体参数见表 1.
|
表 1 模型的介质参数 Table 1 Material properties used in numerical simulation |
渗透率是岩石介质的固有属性.Ingebritsen和Mannning (1999)认为在构造活跃的大陆地壳中，渗透率随地壳深度逐渐减小，遵从经验公式logk≈-14-3.2log (-z)，其中k为渗透率，单位m2，z为深度，单位为km.因此在地表附近渗透率范围为10-11~10-13m2.Shmonov等 (2002)认为渗透率符合公式logk≈-12.56-3.225(-z)0.223，因此地壳40 km深度范围内的渗透率为10-14~10-20m2.由于本文利用商用软件ABAQUS进行计算，需要将渗透率转换为液导率作为模型的输入参数，两者的关系为：
|
(7) |
其中K为液导率，单位为m·s-1；μ为流体动力学黏度，单位Pa·s；ρ为流体密度.本文假定水的黏度为10-3Pa·s，密度为1000 kg·s-1；g重力加速度为9.8 m·s-2，则渗透率与液导率之间的关系大约为107 m-1.因此，假定液导率在整个模型中是均匀且各向同性的，其值始终为10-5 m·s-1.由于本文的重点在于研究孔隙流体在地震过程中引起的应力变化，因此模型仅使用单一完全耦合孔隙弹性介质，且介质为完全饱和的.分别建立走滑型断层、逆冲断层和正断层三种模型，模拟不同类型地震产生的孔隙压力变化，探讨其对静态库仑应力变化的影响.模型空间几何尺寸统一为250 km×250 km×60 km，仅模拟地壳中的情况.
3.1 同震位错及静态库仑应力变化地震是断层内岩石突然滑动造成的，因此可以通过降低断层内部的静摩擦系数来模拟地震.本文对地震的模拟是通过在形成初始应力场后，突然降低断层中央障碍体的摩擦系数从而产生断层同震错动来实现的.断层两盘为接触单元，大部分区域摩擦系数取0.6；而在两盘的中间有一个闭锁程度较高的障碍体 (Asperity patch)，摩擦系数较高 (1.0).先对模型进行加载，断层及障碍体附近累积弹性应变能.在模型达到均衡后，障碍体的摩擦系数突然降至较低的值 (0.5)，此时障碍体产生错动.模型的边界条件可以根据地震学给出的断层破裂尺度及位错大小作为约束条件进行控制 (Wells and Coppersmith, 1994).这种模拟地震的方式在前人的研究中多有使用 (Mikumo and Miyatake, 1979; 蔡永恩等, 1999; Beeler et al., 2001).
定义静态库仑应力变化ΔCFS为 (Harris, 1998)：
|
(8) |
(8) 式中的Δσn+Δp为有效正应力的变化量Δσeffn.与断层错动方向一致的剪切应力和有效正应力的增加使得ΔCFS增加，促进断层破裂，地震更容易发生.
当使用Okada的解析解计算静态库仑应力变化时，由于无法得知孔隙压力变化，常将 (8) 式简化为：Δσf=Δτ-μ′Δσn，其中视摩擦系数μ′=μ(1-B).在这种情况下，视摩擦系数μ′已经不是岩石的固有介质常数，而是融合了岩石摩擦性质和孔隙压力变化两部分的影响.而本文计算中直接采用公式 (8)，利用岩石固有摩擦性质，独立考虑孔隙压力变化的影响.
3.2 走滑型地震走滑型地震的断层模型和边界条件如图 1所示，有限元网格全部为四面体，单元最小边长为200 m，最大边长为2000 m，共由1926874个单元组成.断层AB中央设置一个摩擦系数较高的障碍体，几何尺寸为13×7 km，这与一次MW6.0走滑型地震的破裂尺度大体相当，其有限元网格如图 2所示.障碍体上方距地表5 km，两个接触面摩擦系数为1.0.
|
图 1 走滑型断层的有限元模型及边界条件 Fig. 1 FEM model of strike-slip fault and boundary conditions |
|
图 2 断层面及障碍体的有限元网格示意图 Fig. 2 Sketch showing FEM mesh of strike-slip fault plane and asperity patch |
首先在模型的四个侧面加滑轮，并对整个模型加载重力，同时在地表施加一个6 MPa的压力和2.2 MPa的孔隙压力，用来模拟一个300 m厚的沉积层.当模型在以上初边值条件的作用下达到稳定后，分别固定模型西盘和东盘的北侧及南侧，由模型的东西两侧施加位移向断层面挤压，用来在断层面上产生正应力.同时，在模型西盘的南侧和东盘的北侧对向施加位移，在断层接触面上将产生不均匀的剪应力，弹性应变能将在模型的断层面处累积.当整个模型达到均衡后，断层上下盘的摩擦系数维持0.6不变，障碍体接触面的摩擦系数由1.0降至0.5，此时模型即发生一次右旋走滑型地震事件.调整东西两盘的南北向位移，使断层两盘的障碍体相对位错平均值为0.4 m，该滑移量符合MW6.0右旋走滑型地震产生的滑移特征，如图 3所示.
|
图 3 走滑型地震障碍体位移分布图 Fig. 3 Slip distribution of strike-slip fault′s asperity patch |
图 4给出了走滑型地震造成的瞬时孔隙压变化，可以看到其图案在空间中呈四象限的正负相间分布.在震源深度上的孔隙压变化幅值是最大的，而距离震源越远，其孔隙压变化越小.另外，孔隙压力在震源位置应该是不变的，而在本文的模拟中，震源位置的孔隙压变化也极小，这与前人的结果保持一致 (Wang, 2000).
|
图 4 走滑型地震造成的孔隙压变化 (a) 10 km深度; (b) 20 km深度. Fig. 4 Pore pressure changes generated by strike-slip earthquake (a) At 10 km depth; (b) At 20 km depth. |
图 5a为使用Coulomb 3.3软件 (Lin and Stein, 2004; Toda et al., 2005) 计算得到的库仑应力变化图案，该软件使用Okada的解析解计算无限半空间完全弹性的介质中给定的位移产生的应力变化，物性参数、几何尺寸及滑移距离与图 1a中有限元模型保持一致.图 5b为考虑流固耦合作用，经过有限元方法模拟得到的结果，即将断层滑动时产生的孔隙压变化Δp带入公式 (2) 中，计算其库仑应力变化.可以看到Okada方法的结果与考虑流固耦合作用时的结果有一定的差异：在考虑流固耦合作用时，静态库仑应力变化在空间中的分布更广泛，颜色更深 (幅值更大) 的区域明显增加.在近场区域图案改变较大，原本红色的应力上升区域变为蓝色的“地震影区”，其作用范围更加广泛.同时，在模型的西北和东南两个部分，蓝色的应力下降区域颜色更深.因此，流固耦合作用对走滑型地震产生的静态库仑应力变化影响较大，近场的静态库仑应力明显下降，降低了其地震危险性.
|
图 5 走滑型地震产生的静态库仑应力变化对比图 (深度10 km) (a) Okada模型得到的结果；(b) 考虑流固耦合作用时得到的结果. Fig. 5 Static Coulomb stress changes of strike-slip earthquake (at 10 km depth) (a) Okada′s method; (b) Fluid-solid coupling model. |
图 6为逆冲型断层的有限元模型及其边界条件.模型网格同样为四面体单元，单元最小边长200 m，最大边长2000 m，共有2 173 221个单元.断层倾角为27°，AB为断层在地表的投影.设置的障碍体中心距地表22 km，几何尺寸为10×7 km，符合MW6.0地震的破裂特征，障碍体接触面的摩擦系数为1.0.模型的边界条件与走滑型地震类似，首先模型在重力和300 m的沉积层作用下达到稳定，然后固定模型的西侧、南侧，并从模型的东侧和北侧施加位移进行挤压.由于断层面与中间的障碍体具有不同的摩擦系数，在东侧的不断挤压过程中将在障碍体上产生应力集中.
|
图 6 逆冲型断层的有限元几何模型及边界条件 Fig. 6 Conceptual model of thrust fault with asperity and boundary conditions |
当整个模型处于均衡状态后，维持断层上下盘的摩擦系数0.6不变，中央障碍体的摩擦系数由1.0突然降至0.5，障碍体产生滑移.不断调整北侧、东侧挤压的位移，使障碍体的平均逆冲位移达到0.52 m，即MW6.0地震产生的滑移特征，如图 7所示.
|
图 7 逆冲型地震障碍体位移分布图 Fig. 7 Slip distribution of thrust fault′s asperity patch |
图 8和9给出了不同剖面和深度上因本次地震造成的瞬时孔隙压变化.图 8a展示的是图 6中的PP′剖面，可以看到断层两边的区域孔隙压变化情况正好相反，以震源为中心，在空间中呈四象限分布，断层下盘整体呈孔隙压下降，仅震源下部较深的区域有一片孔隙压上升的区域.而对于断层的上盘，震源上方有一片孔隙压上升区域，其他区域则为孔隙压下降区域.另外，对比图 8a、b、c可以看出，随着距离断层越来越远，孔隙压变化的幅值逐渐变小.在距离震源20 km的地方，孔隙压变化的最大值仅有0.005 MPa.
|
图 8 逆冲断层滑移造成孔隙压变化的剖面图 (a) 震源中心PP′剖面; (b) PP′剖面往北10 km; (c) PP′剖面往北20 km. Fig. 8 Pore pressure changes on different vertical cross sections (a) Section PP′; (b) 10 km north of PP′ section; (c) 20 km north of PP′ section. |
另外，图 9给出了不同深度上孔隙压的变化情况.可以看到震中附近总是表现为孔隙压上升，上升幅度也较大，呈一椭圆形.在震源深度附近 (20 km)，红色区域最小，而东西两侧的蓝色区域区域最大，且颜色也最深，说明孔隙压在震源深度上主要呈减小趋势.总体来说，震源位置的孔隙压总是增大的，并且主要分布在震源的上下两部分，而断层上下两盘大部分区域为孔隙压减小的区域.
|
图 9 逆冲型断层滑动造成的孔隙压变化 (a) 10 km深度; (b) 20 km深度; (c) 30 km深度; (d) 40 km深度. Fig. 9 Pore pressure changes by thrust fault slipping (a) At 10 km depth; (b) At 20 km depth; (c) At 30 km depth; (d) At 40 km depth. |
流固耦合作用对静态库仑应力变化的影响如图 10所示.由于孔隙压的作用，在震源附近的静态库仑应力增加区域虽然变小了，但是中央的蓝色区域也变小了.另外，通过考察其他深度上考虑流固耦合作用时的静态库仑应力变化后发现，孔隙压变化对逆冲型地震的静态库仑应力变化影响较大，震源附近的“地震影区”在孔隙压力作用下会明显的缩小，而库仑应力增加区域的幅值也会增大.由于近场常常是余震发生的主要区域，且震源深度通常在10~25 km之间，因此，逆冲型地震产生的孔隙压变化对后续余震有一定的触发作用.
|
图 10 逆冲型断层产生的静态库仑应力变化对比图 (a) Okada模型计算结果，深度20 km；(b) 考虑流固耦合作用的结果，深度20 km; (c) Okada模型计算结果，深度30 km；(d) 考虑流固耦合作用的结果，深度30 km. Fig. 10 Static Coulomb stress changes of thrust earthquake (a) Okada′s method at 20 km depth; (b) Fluid-solid coupling model at 20 km depth; (c) Okada′s method at 30 km depth; (d) Fluid-solid coupling model at 30 km depth. |
正断层的有限元模型及边界条件如图 11所示，断层倾角为63°，共由1968127个四面体单元组成.断层中障碍体的几何尺寸为12×8.5 km，其中心距地表10 km.从图 10b中可以看到，该模型的边界条件与逆冲型断层模型的边界条件类似，在重力和沉积层的作用下模型均衡后，在模型的西侧和北侧加载位移进行挤压.不同的是，在该模型西侧施加的位移远远小于逆冲型断层模型中的位移，使得断层附近的最大主应力为y向的挤压.此时，断层的两盘位移呈正断层的特征，即上盘向下滑移，下盘向上滑移.令障碍体的摩擦系数由1.0变为0.5，而断层接触面的摩擦系数保持0.6不变，此时障碍体发生滑动，其平均位移为0.36 m，即MW6.0正断层型地震，如图 12所示.
|
图 11 正断型断层的几何模型及边界条件 Fig. 11 Conceptual model of normal fault model with asperity and boundary conditions |
|
图 12 正断型地震障碍体位移分布图 Fig. 12 Slip distribution of normal fault′s asperity patch |
图 13展示了正断层型地震前后不同剖面上的孔隙压变化，将其对比图 8可知，正断层型地震产生的孔隙压变化图案与逆冲断层产生的图案类似，在空间中呈四象限分布，但正负区域正好相反.图 13a为模型中心PP′剖面的孔隙压分布图，可以看到断层上盘的极大部分区域呈孔隙压上升，仅在震源上方为蓝色的下降区；同时对于断层的下盘，仅在震源下方呈孔隙压下降，而震源深度以上都是红色的孔隙压上升区域.另外由图 13的a、b、c可以看出，孔隙压变化的在空间中的图案样式并不随距离的变化而变化，仅体现为幅值下降.在距离震源20 km的位置，孔隙压变化最大值为0.008 MPa.
|
图 13 正断型断层滑移造成孔隙压变化的剖面图 (a) PP′剖面; (b) PP′剖面往北10 km; (c) PP′剖面往北20 km. Fig. 13 Pore pressure changes of normal faulting slip at different vertical cross sections (a) Section PP′; (b) 10 km north of PP′ section; (c) 20 km north of PP′ section. |
不同深度上孔隙压的变化情况如图 14所示.由图 14a可知在震源上方出现一片蓝色的孔隙压下降，且幅值较大；而在10 km深度上 (震源附近) 基本呈孔隙压上升 (图 14b).在震源深度以上，孔隙压下降区域位于上盘.随着深度的不断下降 (图 14c、d)，蓝色的负值区域逐渐向西移动，并且区域逐渐扩大.总体来说，正断型地震产生的孔隙压变化在震源附近是增大的，而在上下两侧主要是孔隙压减小的区域.
|
图 14 正断型断层滑动造成的孔隙压变化 (a) 5 km深度; (b) 10 km深度; (c) 15 km深度; (d) 20 km深度. Fig. 14 Pore pressure changes of normal fault slipping (a) At 5 km depth; (b) At 10 km depth; (c) At 15 km depth; (d) At 20 km depth. |
考虑流固耦合作用时，正断型地震产生的孔隙压变化对其静态库仑应力变化的影响较小，如图 15所示.在10 km深度上 (图 15a)，考虑流固耦合作用时的静态库仑应力变化与Okada模型的结果极为相似，仅震源部分蓝色区域范围扩大；同时，对于20 km深度，近场的库仑应力变化在流固耦合作用下，断层下盘的蓝色“地震影区”变大，而断层上盘的地震影区也并没有因为孔隙压上升 (图 15d) 而变得更小，仅红色区域的幅值变大.另外，通过考察各深度上的静态库仑应力变化图案可知，大部分情况下流固耦合作用仅会使库仑应力减小的区域增大.因此，正断型地震产生的孔隙压变化并不能使静态库仑应力变化产生较大的变化，对余震的触发效果并不明显.
|
图 15 正断型断层产生的静态库仑应力变化对比图 (a) Okada模型计算结果，深度15 km；(b) 考虑流固耦合作用的结果，深度15 km; (c) Okada模型计算结果，深度20 km；(d) 考虑流固耦合作用的结果，深度20 km. Fig. 15 Static Coulomb stress changes of normal earthquake (a) Okada′s method at 15 km depth; (b) Fluid-solid coupling model at 15 km depth; (c) Okada′s method at 20 km depth; (d) Fluid-solid coupling model at 20 km depth. |
本文基于完全耦合孔隙弹性理论，利用三维有限元方法，分别建立走滑型断层、逆冲断层和正断层三种模型，模拟了地震发生前后的应力状态，考察了Okada模型与流固耦合模型的结果差异.由于完全耦合孔隙弹性理论计算起来较为复杂，前人的相关研究中通常对地震产生的孔隙压变化的作用进行简化，使用Skempton系数B代替孔隙压场的变化，简化为非耦合问题，如Toda等 (2008)，Lin等 (2011)，Luo和Liu (2010)；或对两者分别计算，然后进行加和，简化为假耦合问题，如Lei (2011)，Deng等 (2010)；或只考虑孔隙压变化对岩石骨架的影响，不考虑岩石应力场对孔隙压变化的作用，简化为单项耦合问题，如Gahalaut (2003, 2008)，孙玉军等 (2012)；而周斌等 (2010)和陶玮等 (2014)利用完全耦合孔隙弹性理论对紫坪铺水库蓄水与汶川地震的关系进行了计算，但仅对逆冲型地震进行了研究.
Muir-Wood和King (1993)认为，地下流体受到的影响与地震震级相关性较小，而对于不同类型的地震所产生的断层滑移对流体的影响较大.例如1983年美国Idaho州M6.9的正断层型地震，造成的影响远大于1906年旧金山M7.8的走滑型地震.本文虽然未对不同震级的地震造成的孔隙压变化进行研究，但按照经验公式给出的地震破裂尺度及位错大小 (Wells and Coppersmith, 1994)，建立了三种不同类型的断层模型，模拟同样震级 (MW6.0) 不同类型的地震所产生的断层滑移.表 2给出了三种类型地震造成的孔隙压力变化，可以看到走滑型地震产生的孔隙压力变化远远小于逆冲型及正断型地震，因此，本文的结果与前人保持一致.
|
|
表 2 同震级不同类型的地震造成的孔隙压变化 Table 2 Pore pressure changes of different types of earthquakes with same magnitude |
地震造成的流体变化通常被认为有两种机制 (Rojstaczer and Wolf, 1992; Rojstaczer et al., 1995; Brodsky et al., 2003; Wang et al., 2004, 2013)：(1) 应力场的调整造成了多孔弹性介质中孔隙水的排出；(2) 断裂的张开或闭合导致渗透率的变化.由于本文的研究重点在于孔隙流体在地震过程中引起的应力变化，因此计算中仅考虑了第一种机制：使用各向同性介质，渗透率 (液导率) 在整个模型中取同一值.而实际上，地壳中岩石的渗透率是极为复杂的，与研究区域的岩石类型和主要构造有极为密切的关系.由于局部区域的渗透率极大或极小，造成该区的背景孔隙压极大或极小，当地震发生时，孔隙压的变化将呈更加复杂的图案，对静态库仑应力变化的影响也会不同.另外，地震的发生会使先存裂隙发生变化，使其张开或更加闭锁，渗透率发生变化.并且，在震后的相当长一段时间内，模型由同震时的不排水状态 (Undrained state) 调整为排水状态 (Drained state)，孔隙流体在震后的传播会使得静态库仑应力变化进一步发生改变.如周斌等 (2014)对于龙滩水库附近的断层使用了不同的渗透率，水库加卸载后渗透作用的滞后响应，使得各个地震丛的活动时间及强度都不同.因此，在今后的研究中应对不同的地区使用不同的渗透率，更为细致的研究孔隙压变化对地震的影响.
本文分别建立了三种类型的断层模型计算了同震时的孔隙压变化，并考察了其对静态库仑应力变化的影响.然而实际上，在震后相当长的一段时间内，由于地震的发生改变了区域应力场，使得孔隙流体在地下介质的传播变得更为复杂，而孔隙流体的传播又将反作用于应力场，使其进一步变化.因此，孔隙压场在震后随时间的演化对于我们研究震后静态库仑应力变化，触发后续余震有一定的意义.如陶玮等 (2014)计算了紫坪铺水库蓄水后龙门山断裂带附近的剪切应力、有效正应力和库仑应力的分布及其随时间的演化，三者在空间中不同时间段都有较大的变化.另外，地震之后一段时间内都会出现余滑 (afterslip)，对震源及周围的应力状态都会产生影响.由震后余滑产生的应力扰动引起孔隙压力的变化，有待今后更深入的研究.
4.2 初步结论通过以上数值模拟计算及分析，得出以下初步结论：
本文运用完全耦合孔隙弹性介质的断层滑移模型模拟了地震时流-固耦合作用对静态库仑应力变化的影响，并将结果与传统算法相比较.走滑型断层产生的孔隙压变化图案不随深度变化而变化，流固耦合作用对走滑型地震产生的静态库仑应力变化影响较大，近场的静态库仑应力明显下降，降低了危险性；逆冲型断层滑动产生的孔隙压增加主要分布在震源上下方，会使震源附近的静态库仑应力下降区域变小，更好的触发余震；而正断型地震产生的孔隙压变化则正好相反，增大了震源附近的应力影区范围，这样会降低该区域余震发生的概率.因此，在利用库仑应力变化研究地震触发时，应考虑流-固耦合作用，使得库仑模型的预测结果更符合实际.
致谢感谢匿名审稿专家的宝贵意见.
| Beeler N M, Lockner D L, Hickman S H. 2001. A simple stick-slip and creep-slip model for repeating earthquakes and its implication for microearthquakes at Parkfield. Bull. Seismol. Soc. Am., 91(6): 1797-1804. DOI:10.1785/0120000096 | |
| Biot M A. 1941. General theory of three-dimensional consolidation. Journal of Applied Physics, 12(2): 155-164. DOI:10.1063/1.1712886 | |
| Biot M A. 1956. Theory of propagation of elastic waves in a fluid-saturated porous solid. I. Low-frequency range. The Journal of the acoustical Society of america, 28(2): 168-178. DOI:10.1121/1.1908239 | |
| Bosl W J, Nur A. 2002. Aftershocks and pore fluid diffusion following the 1992 Landers earthquake. J. Geophys. Res., 107(B12): ESE 17-1-ESE 17-12. DOI:10.1029/2001JB000155 | |
| Brodsky E E, Roeloffs E, Woodcock D, et al. 2003. A mechanism for sustained groundwater pressure changes induced by distant earthquakes. J. Geophys. Res., 108(B8): 2390. DOI:10.1029/2002JB002321 | |
| Cai Y, He T, Wang R. 1999. Modelling the source dynamic process of the 1976 Tangshan Earthquake. Acta Seismologica Sinica (in Chinese), 21(5): 469-477. | |
| Coussy O. Poromechanics. Oxford, England: John Wiley & Sons Ltd., 2004. | |
| Deng J S, Sykes L R. 1997. Evolution of the stress field in southern California and triggering of moderate-size earthquakes:A 200-year perspective. J. Geophys. Res., 102(B5): 9859-9886. DOI:10.1029/96JB03897 | |
| Deng K, Zhou S Y, Wang R, et al. 2010. Evidence that the 2008 MW7.9 Wenchuan earthquake could not have been induced by the Zipingpu Reservoir. Bull. Seismol. Soc. Am., 100(5B): 2805-2814. DOI:10.1785/0120090222 | |
| Gahalaut K, Gahalaut V K, Kayal J R. 2008. Poroelastic relaxation and aftershocks of the 2001 Bhuj earthquake, India. Tectonophysics, 460(1-4): 76-82. DOI:10.1016/j.tecto.2008.07.004 | |
| Gahalaut V K, Kalpna, Raju P S. 2003. Rupture mechanism of the 1993 Killari earthquake, India:constraints from aftershocks and static stress change. Tectonophysics, 369(1-2): 71-78. DOI:10.1016/S0040-1951(03)00135-5 | |
| Harris R A. 1998. Introduction to special section:Stress triggers, stress shadows, and implications for seismic hazard. J. Geophys. Res., 103(B10): 24347-24358. DOI:10.1029/98JB01576 | |
| Hart D J, Wang H F. 1995. Laboratory measurements of a complete set of poroelastic moduli for Berea sandstone and Indiana limestone. Journal of Geophysical Research:Solid Earth, 100(B9): 17741-17751. DOI:10.1029/95JB01242 | |
| He J K, Peltzer G. 2010. Poroelastic triggering in the 9-22 January 2008 Nima-Gaize (Tibet) earthquake sequence. Geology, 38(10): 907-910. DOI:10.1130/G31104.1 | |
| Henkel D J, Wade N H. 1966. Plane strain tests on a saturated remolded clay. Journal of Soil Mechanics & Foundations Div, 92: 67-80. | |
| Ingebritsen S E, Manning C E. 1999. Geological implications of a permeability-depth curve for the continental crust. Geology, 27(12): 1107-1110. DOI:10.1130/0091-7613(1999)027<1107:GIOAPD>2.3.CO;2 | |
| Kennedy B M, Kharaka Y K, Evans W C, et al. 1997. Mantle fluids in the San Andreas fault system, California. Science, 278(5341): 1278-1281. DOI:10.1126/science.278.5341.1278 | |
| King G C P, Stein R S, Lin J. 1994. Static stress changes and the triggering of earthquakes. Bull. Seismol. Soc. Am., 84(3): 935-953. | |
| Lei X L. 2011. Possible roles of the Zipingpu Reservoir in triggering the 2008 Wenchuan earthquake. Journal of Asian Earth Sciences, 40(4): 844-854. DOI:10.1016/j.jseaes.2010.05.004 | |
| Lin J, Stein R S. 2004. Stress triggering in thrust and subduction earthquakes and stress interaction between the southern San Andreas and nearby thrust and strike-slip faults. J. Geophys. Res., 109(B2): B02303. | |
| Lin J, Stein R S, Meghraoui M, et al. 2011. Stress transfer among en echelon and opposing thrusts and tear faults:Triggering caused by the 2003 MW=6.9 Zemmouri, Algeria, earthquake. J. Geophys. Res., 116(B3): B03305. | |
| Luo G, Liu M. 2010. Stress evolution and fault interactions before and after the 2008 Great Wenchuan earthquake. Tectonophysics, 491(1-4): 127-140. DOI:10.1016/j.tecto.2009.12.019 | |
| Miao M, Zhu S B. 2012. A study of the impact of static Coulomb stress changes of megathrust earthquakes along subduction zone on the following aftershocks. Chinese J. Geophys. (in Chinese), 55(9): 2982-2993. DOI:10.6038/j.issn.0001-5733.2012.09.017 | |
| Miao M, Zhu S B. 2013. The static Coulomb stress change of the 2013 Lushan MS7.0 earthquake and its impact on the spatial distribution of aftershocks. Acta Seismologica Sinica (in Chinese), 35(5): 619-631. | |
| Mikumo T, Miyatake T. 1979. Earthquake sequences on a frictional fault model with non-uniform strengths and relaxation times. Geophysical Journal International, 59(3): 497-522. DOI:10.1111/gji.1979.59.issue-3 | |
| Muir-Wood R, King G C P. 1993. Hydrological signatures of earthquake strain. J. Geophys. Res., 98(B22): 22035-22068. | |
| Nur A, Booker J R. 1972. Aftershocks caused by pore fluid flow?. Science, 175(4024): 885-887. DOI:10.1126/science.175.4024.885 | |
| Okada Y. 1985. Surface deformation due to shear and tensile faults in a half-space. Bull. Seismol. Soc. Am., 75(4): 1135-1154. | |
| Okada Y. 1992. Internal deformation due to shear and tensile faults in a half-space. Bull. Seismol. Soc. Am., 82(2): 1018-1040. | |
| Parsons T, Ji C, Kirby E. 2008. Stress changes from the 2008 Wenchuan earthquake and increased hazard in the Sichuan basin. Nature, 454(7203): 509-510. DOI:10.1038/nature07177 | |
| Rice J R, Cleary M P. 1976. Some basic stress diffusion solutions for fluid-saturated elastic porous media with compressible constituents. Reviews of Geophysics, 14(2): 227-241. DOI:10.1029/RG014i002p00227 | |
| Roeloffs E. 1996. Poroelastic techniques in the study of earthquake-related hydrologic phenomena. Advances in Geophysics, 37: 135-195. DOI:10.1016/S0065-2687(08)60270-8 | |
| Rojstaczer S, Wolf S. 1992. Permeability changes associated with large earthquakes:An example from Loma Prieta, California. Geology, 20(3): 211-214. DOI:10.1130/0091-7613(1992)020<0211:PCAWLE>2.3.CO;2 | |
| Rojstaczer S, Wolf S, Michel R. 1995. Permeability enhancement in the shallow crust as a cause of earthquake-induced hydrological changes. Nature, 373(6511): 237-239. DOI:10.1038/373237a0 | |
| Schoenball M, Müller T M, Müller B I R, et al. 2010. Fluid-induced microseismicity in pre-stressed rock masses. Geophysical Journal International, 180(2): 813-819. DOI:10.1111/gji.2010.180.issue-2 | |
| Shi Y L, Cao J L. 2010. Some aspects in static stress change calculation——case study on Wenchuan earthquake. Chinese J. Geophys. (in Chinese), 53(1): 102-110. DOI:10.3969/j.issn.0001-5733.2010.01.011 | |
| Shmonov V M, Vitovtova V M, Zharikov A, et al. 2002. Fluid permeability of the continental crust:Estimation from experimental data. Geochemistry International, 40: S3-S13. | |
| Skempton A W. 1954. The pore-pressure coefficients A and B. Geotechnique, 4(4): 143-147. DOI:10.1680/geot.1954.4.4.143 | |
| Stein R S, Barka A A, Dieterich J H. 1997. Progressive failure on the North Anatolian fault since 1939 by earthquake stress triggering. Geophysical Journal International, 128(3): 594-604. DOI:10.1111/gji.1997.128.issue-3 | |
| Sun Y J, Zhang H, Dong S W, et al. 2012. Study on effect of the Zipingpu reservoir on the occurrence of the 2008 Wenchuan earthquake based on a 3D-poroelastic model. Chinese J. Geophys. (in Chinese), 55(7): 2353-2361. DOI:10.6038/j.issn.0001-5733.2012.07.020 | |
| Tao W, Masterlark T, Shen Z K, et al. 2014. Triggering effect of the Zipingpu Reservoir on the 2008 MW7.9 Wenchuan, China, Earthquake due to poroelastic coupling. Chinese J. Geophys. (in Chinese), 57(10): 3318-3331. DOI:10.6038/cjg20141019 | |
| Toda S, Stein R S, Richards-Dinger K, et al. 2005. Forecasting the evolution of seismicity in southern California:Animations built on earthquake stress transfer. J. Geophys. Res., 110(B5): B05S16. | |
| Toda S, Lin J, Meghraoui M, et al. 2008. 12 May 2008 M=7.9 Wenchuan, China, earthquake calculated to increase failure stress and seismicity rate on three major fault systems. Geophys. Res. Lett., 35(17): L17305. DOI:10.1029/2008GL034903 | |
| Wan Y G, Shen Z K. 2010. Static Coulomb stress changes on faults caused by the 2008 MW7.9 Wenchuan, China earthquake. Tectonophysics, 491(1-4): 105-118. DOI:10.1016/j.tecto.2010.03.017 | |
| Wang C Y, Wang C H, Manga M. 2004. Coseismic release of water from mountains:Evidence from the 1999 (MW=7.5) Chi-Chi, Taiwan, earthquake. Geology, 32(9): 769-772. DOI:10.1130/G20753.1 | |
| Wang H F. 1993. Quasi-static poroelastic parameters in rock and their geophysical applications.//Experimental Techniques in Mineral and Rock Physics. Basel:Birkhäuser, 269-286. | |
| Wang H F. Theory of Linear Poroelasticity with Applications to Geomechanics and Hydrogeology. Princeton, NJ: Princeton University Press, 2000. | |
| Wang M M, Jia D, Lin A M, et al. 2013. Late Holocene activity and historical earthquakes of the Qiongxi thrust fault system in the southern Longmen Shan fold-and-thrust belt, eastern Tibetan Plateau. Tectonophysics, 584: 102-113. DOI:10.1016/j.tecto.2012.08.019 | |
| Wells D L, Coppersmith K J. 1994. New empirical relationships among magnitude, rupture length, rupture width, rupture area, and surface displacement. Bull. Seismol. Soc. Am., 84(4): 974-1002. | |
| Zhang L. Engineering Properties of Rocks. New York: Elsevier Health Sciences, 2005. | |
| Zhao D P, Kanamori H, Negishi H, et al. 1996. Tomography of the source area of the 1995 Kobe earthquake:Evidence for fluids at the hypocenter?. Science, 274(5294): 1891-1894. DOI:10.1126/science.274.5294.1891 | |
| Zhao D P, Ochi F, Hasegawa A, et al. 2000. Evidence for the location and cause of large crustal earthquakes in Japan. J. Geophys. Res., 105(B6): 13579-13594. DOI:10.1029/2000JB900026 | |
| Zhou B, Sun F, Yan C H, et al. 2014. 3D-poreelastic finite element numerical simulation of Longtan reservoir-induced seismicity. Chinese J. Geophys. (in Chinese), 57(9): 2846-2868. DOI:10.6038/cjg20140911 | |
| Zhou B, Xue S F, Deng Z H, et al. 2010. Relationship between the evolution of reservoir-induced seismicity in space-time and the process of reservoir water body load-unloading and water infiltration-a case study of Zipingpu reservoir. Chinese J. Geophys. (in Chinese), 53(11): 2651-2670. DOI:10.3969/j.issn.0001-5733.2010.11.013 | |
| Zhou X J, Burbey T J. 2014. Pore-pressure response to sudden fault slip for three typical faulting regimes. Bull. Seismol. Soc. Am., 104(2): 793-808. DOI:10.1785/0120130139 | |
| Zhu S B, Miao M. 2015. How did the 2013 Lushan Earthquake (MS=7.0) trigger its aftershocks? Insights from static coulomb stress change calculations. Pure and Applied Geophysics, 172(10): 2481-2494. DOI:10.1007/s00024-015-1064-3 | |
| 蔡永恩, 何涛, 王仁. 1999. 1976年唐山地震震源动力过程的数值模拟. 地震学报, 21(5): 469–477. | |
| 缪淼, 朱守彪. 2012. 俯冲带上特大地震静态库仑应力变化对后续余震触发效果的研究. 地球物理学报, 55(9): 2982–2993. DOI:10.6038/j.issn.0001-5733.2012.09.017 | |
| 缪淼, 朱守彪. 2013. 2013年芦山MS7.0地震产生的静态库仑应力变化及其对余震空间分布的影响. 地震学报, 35(5): 619–631. | |
| 石耀霖, 曹建玲. 2010. 库仑应力计算及应用过程中若干问题的讨论——以汶川地震为例. 地球物理学报, 53(1): 102–110. DOI:10.3969/j.issn.0001-5733.2010.01.011 | |
| 孙玉军, 张怀, 董树文, 等. 2012. 利用三维孔隙弹性模型探讨紫坪铺水库对汶川地震的影响. 地球物理学报, 55(7): 2353–2361. DOI:10.6038/j.issn.0001-5733.2012.07.020 | |
| 陶玮, MasterlarkT, 沈正康, 等. 2014. 紫坪铺水库造成孔隙弹性耦合变化及其对2008年汶川地震触发作用. 地球物理学报, 57(10): 3318–3331. DOI:10.6038/cjg20141019 | |
| 周斌, 孙峰, 阎春恒, 等. 2014. 龙滩水库诱发地震三维孔隙弹性有限元数值模拟. 地球物理学报, 57(9): 2846–2868. DOI:10.6038/cjg20140911 | |
| 周斌, 薛世峰, 邓志辉, 等. 2010. 水库诱发地震时空演化与库水加卸载及渗透过程的关系——以紫坪铺水库为例. 地球物理学报, 53(11): 2651–2670. DOI:10.3969/j.issn.0001-5733.2010.11.013 | |