中国科学院大学 中国科学院计算地球动力学重点实验室, 北京 100049
2016年03月07日 收稿; 2016年04月28日 收修改稿
基金项目: 国家自然科学基金重大项目（41086965）与中国科学院国际合作创新团队项目（KZZD-EW-TZ-19）资助
通信作者: 石耀霖, E-mail:shiyl@ucas.ac.cn
摘要: 2004年12月26日苏门答腊Mw9.3地震产生的应力扰动，强烈影响了苏门答腊及其周边地区的地震活动性，甚至对中国西南部地区的地震活动性也有明显的影响。基于断层滑动模型，采用等效体力处理断层错动有限单元方法，分别用分层球对称地球模型与三维横向不均匀性地球模型计算在全球地壳和地幔中产生的同震效应。据此，对在苏门答腊和中国西南等地区产生的同震效应进行重点分析。结果表明：在同震破裂面附近，使用三维不均匀性地球模型与使用分层球对称模型得到的结果有明显差别，前者更加合理。在同震破裂面附近有必要使用三维不均匀性地球模型。使用球对称模型计算，有67.6%的余震发生在库仑应力增加区；而使用横向不均匀性地球模型时，计算的这一比率可高达72.3%，这说明余震与断层库仑应力增加区域有很好的对应关系。在西安达曼断层北段余震的发生与库仑应力为正值的主要断层有很好的对应关系，余震受到其附近的主要断层的控制作用。
关键词: 2004年苏门答腊地震有限元方法三维横向不均匀性库仑应力地震触发
Coulomb stresses induced by the 2004 Mw9.3 Sumatra earthquake
QU Wulin, ZHANG Bei, HUANG Luyuan, CHENG Huihong, ZHANG Huai, SHI Yaolin
Key Laboratory of Computational Geodynamics of Chinese Academy of Sciences, University of Chinese Academy of Sciences, Beijing 100049, China
Abstract: The stress perturbation induced by 2004 Sumatra earthquake has greatly changed seismicity around Sumatra and its vicinity, including Southwest China. In this study, we calculated the co-seismic stress changes induced by this earthquake using a layered spherically symmetric earth model and a 3D laterally heterogeneous earth model. Both the calculations are based on the fault slip model, which fits well with the GPS observation and has been widely quoted by other researchers, and the equivalent body force finite element method is applied to deal with the co-seismic fault slip effect. The entire volume of a spherical shell from earth surface to core-mantle boundary is included in our investigation, and boundary conditions can be easily assigned to the two boundaries. In the spherical symmetric model, there are 67.6% of aftershocks occurring in regions of positive Coulomb stresses. In the lateral heterogeneous model, 72.3% of aftershocks occurring in regions of positive Coulomb stresss. The aftershocks in the northern part of west Andaman fault have good correlation with the main faults which have positive Coulomb stress. So the aftershocks were apparently controlled by adjacent main faults.
Key words: 2004 Sumatra earthquakefinite element method3D lateral heterogeneityCoulomb stressearthquake triggering
2004年12月26日在苏门答腊地区发生的矩震级为9.1到9.3级的地震[1]是现代全球宽频带、高动态范围地震仪与全球定位系统(GPS)铺设后, 记录到的第一个高质量地球物理数据的大地震。地震从锡默卢岛(Simeulue)西北部海域开始破裂, 沿着巽他海沟(Sunda trench)向北传播到安达曼群岛(Andaman Islands)(图 1)。破裂面长度为1 200~1 500 km，破裂速度为2.0~2.8 km/s, 持续时间超过500 s, 破裂面宽度为150~200 km[1-7]。其断层破裂长度与破裂持续时间都是有记录以来之最[6, 8]。该地震引起的灾害程度很大, 影响范围也很广。
Fig. 1
 | Download: JPG larger image |
图 1 苏门答腊地区构造背景与余震分布 Fig. 1 Tectonic context of Sumatra-Andaman and aftershocks of the 2004 earthquake | |
图 1为苏门答腊地区构造背景与余震分布。红色大箭头与红色数字表示印度板块相对巽他板块的运动速率; 品红色线表示板块边界; 红色带三角形线是巽他海沟俯冲断层; 红色细线是苏门答腊走滑断层; 圆点代表震级大于4的从2004年12月26日到12月31日的余震。数据来源:余震数据NEIC; 板块边界数据来自Bird[9]。断层数据改自Subarya等[10]; 印度板块相对巽他板块的运动速率引自Subarya等[10], 由Bock等[11]的区域动力学模型将计算得来; 洋壳年龄数据来自Muller等[12]。
苏门答腊地震对周围地区产生了明显的影响, 在川滇地区观测到明显的同震位移[13]。该大地震触发了周围区域的地震活动, 震后当天云南与缅甸就分别触发5.0级与5.6级地震[14]。在云南发生了持续约15天的高频地震活动[15]; 柴达木盆地的地震活动也明显增加[16]。
Cattin等[17]计算苏门答腊地震对周边断层的库仑应力的影响, 并探讨对余震的触发效应。缪淼和朱守彪[18]计算得到的余震处于库仑应力增加区域只有49.8%, 远小于板内地震对余震的触发比例。其他****[14, 19]也对苏门答腊地震对周边及更远地区的应力改变与地震触发做了详细的探讨。但是他们都忽略了地球曲率、地球分层介质与横向不均匀性的影响。
苏门答腊地震发生于俯冲带, 岩石物性参数从洋壳到陆壳变化十分剧烈, 有必要在计算中考虑岩石的横向不均匀性。目前广泛用于计算同震效应的方法以Okada[20]的半无限解析解和Wang等[21]的分层模型为代表, 但在这两种方法中都不能考虑真实地球的横向不均匀性。Fu和Sun[22]基于微扰动理论, 研究三维不均匀性地球模型的形变, 但实际地球模型的不均匀性已超出微扰动模型所能接受的范围。张贝等[23]开发的等效体力并行有限元方法可考虑横向不均匀性等复杂情况。
本文利用等效体力并行有限元方法[23], 计算苏门答腊地震引起的全球同震位移与应力变化, 然后计算主要断层和余震处的库仑应力, 据此探讨苏门答腊地震对周边主要断裂的应力变化和对余震的触发作用。
1 计算方法和参数采用张贝等[23]提出的等效体力有限元方法, 计算全球的同震位移与应力。在该方法中, 将断层上的位错量等效为体力, 添加到平衡方程的右端项。对核幔边界以上至地表的球壳部分采用自适应有限元网格剖分。在二维与三维问题的计算中, 在保证计算精度的前提下, 自适应网格分别节省100倍与1 000倍的计算量[24], 能明显提高计算效率。离发震断层越远, 同震位移的梯度越小, 单个网格分辨率就越低, 网格越粗糙。如图 2所示, 有限元网格在发震断层附近较密, 其大小约4 km×4 km×4 km; 网格尺寸随距发震断层距离的增大而逐渐增大, 在中国西南地区的网格大小约32 km×32 km×32 km。整个地球的地壳和地幔剖分的单元数量共约480万。
Fig. 2
 | Download: JPG larger image |
图 2 全球网格剖分 Fig. 2 Global mesh | |
半无限空间均匀模型与分层模型都忽略了地球曲率的影响, 但研究显示地球曲率对远场有显著的影响, 对走滑地震的影响可达20%, 而对近场的影响可以忽略不计[25]。Sun和Okubo[26]对比半无限空间介质与球形地球模型的同震位移后, 发现曲率与层状构造的影响都不容忽视。Dong等[27]发现地球曲率对同震位移的影响为5%, 而分层模型的影响达到25%, 对重力的影响为11%。总之, 在近场可以使用半无限空间均匀模型与分层模型, 但对于远场的同震变化, 地球曲率的影响不能忽略。
另外, 巽他海沟西侧是大洋地壳, 而向东跨过海沟后则为大陆地壳(图 1), 岩石的物性参数变化很大, 因此有必要在计算中考虑地球模型的横向不均匀性。Fu和Sun[22]基于微扰动理论, 研究三维不均匀性地球模型的形变。但实际的横向不均匀性已超出微扰动模型所能接受的范围。Qu等[28]探讨三维横向不均匀性地球模型对2004年苏门答腊地震的同震位移的影响, 发现三维横向不均匀性地壳结构对同震水平、垂直位移的影响分别为23%和40%, 因此三维横向不均匀性地球模型的影响不可忽略。
本文在考虑地球曲率的情况下, 分别计算分层球对称地球模型与全球横向不均匀性模型下的同震效应。分层球对称模型采用的是PREM模型[29]。全球三维横向不均匀性地球模型由两部分组成, 在地壳部分使用Crust 1.0[30]模型, 在地幔部分使用PyGSuM[31]横向不均匀性模型。在核幔边界采用弹簧边界条件, 由核幔密度差值决定; 地表为自由边界条件, 即应力为零。
我们使用Chlieh等[6]给出的苏门答腊地震断层滑动模型。瞿武林等[32]对比研究不同的断层滑动模型, 得出Chlieh模型(图 3)不仅在近场与观测数据符合比较好, 并且在全球的同震位移也甚为合理。并且该模型也被广泛应用于一些研究中[15, 17-18, 33]。因此我们在计算中使用Chlieh等[6]给出的断层滑动模型。该断层滑动模型是由近场GPS数据、近场珊瑚礁测高数据、近场遥感影像数据以及远场GPS数据等作为观测, 使用3个平面来拟合发震断层的几何面(图 3)反演得到的。如图 3(a)所示, 断层滑动位移在中段与南段较大, 在北段较小; 中段6.5°N附近取得最大滑动位移为17 m。断层的中段与南段滑动矢量基本垂直于海沟(图 3(b)), 而北段带有部分右旋走滑分量, 这与印度洋板块相对于缅甸微板块的斜向俯冲有关。从尼克巴群岛(Nicobar Islands)到苏门答腊岛西南海域处的右旋走滑分量基本上都发生在苏门答腊大断层(Great Sumatra Fault)上, 而俯冲带上基本只有逆冲分量[34]。
Fig. 3
 | Download: JPG larger image |
图 3 断层滑动模型(Chlieh等[6]) Fig. 3 Fault slip model (Chlieh et al.[6]) | |
2 同震应力变化虽然我们计算的是苏门答腊地震对全球的同震影响, 但由于距发震断层越远, 同震影响越小, 故这里只显示苏门答腊地震引起周边地区在15 km深度的同震应力变化(图 4)。图 4中, 南北向和东西向正应力以拉张为正, 压缩为负; 图中数字为应力的量级, 红色表示应力为正值, 蓝色表示应力为负值, 图中等值线数字表示应力的量级, 负号表示应力为负。与同震位移一样, 同震应力变化随着离开断层距离的增加而迅速减小。在断层破裂面及其附近, 同震应力变化超过1 MPa。在中国西南云南地区的同震应力变化有数千帕, 向北迅速减小, 只有数百帕。同震应力变化超过1 kPa的范围相当大, 包括断层破裂面及其附近区域、中国云南地区与中南半岛大部分区域。南北向水平正应力σNN(图 4(a))与水平剪应力σNE(图 4(c))成正负相间, 东西向水平正应力σEE(图 4(b))在断层破裂面以及断层破裂面南北方向上为负值, 在破裂面东西方向为正值。
Fig. 4
 | Download: JPG larger image |
| (a)与(d)南北向正应力σNN; (b)与(e)东西方向正应力σEE; (c)与(f)水平剪应力σNE。(a)-(c)使用分层球对称介质; (d)-(f)使用三维横向不均匀性介质。 图 4 苏门答腊地震引起的15 km深度同震应力变化 Fig. 4 Coseismic stress changes induced by Sumatra earthquake at a depth of 15 km | |
为更清晰地看到使用分层球对称地球模型与三维横向不均匀性地球模型所造成的近场同震应力变化差异(由分层球对称模型得到的应力分量减去由三维横向不均匀性地球模型得到的相应应力分量), 我们将近场部分(图 4(a)中方框区域部分)放大, 如图 5所示。图 5中黑色方框为2004年苏门答腊Mw9.3地震的破裂范围[6], 绿色方框为2005年3月28日尼亚斯(Nias)地震的破裂范围[35]。两个模型造成的应力分量在正负分布形状上区别不大(图 5(a)与(d), (b)与(e), (c)与(f)), 但二者的差值(图 5(g)-(i))却较大, 甚至在断层破裂面范围内超过1 MPa。大的差异集中在断层破裂面中部尼克巴群岛附近与南部锡默卢岛西北区域, 这些地区附近的断层滑动量也相应较大(图 3(a))。地壳和地幔的横向不均匀非常强, 这对结果也有明显影响。类似于前人对同震位移差异的统计[22, 26], 以分层球对称地球模型计算得到同震应力变化中绝对值最大为标准, 计算同震应力变化各个分量影响的百分比。图 5中所示为15 km深度的水平切面以内, 南北方向同震正应力σNN的最大差异为137%, 东西方向同震正应力σEE的最大差异为201%, 水平剪切应σNE的差异为65%。而三维全空间中, 同震应力变化的各个分量的最大差异更大, 对应的分别为147%, 243%和110%。
Fig. 5
 | Download: JPG larger image |
| 显示图 4(a)中方框区域。(a)-(f)同图 4; (g)-(i)分层球形模型与三维横向不均匀性模型同震正应力差异。(g)南北向同震正应力σNN差异; (h)东西向同震正应力σEE差异; (i)水平同震剪应力σNE差异。 图 5 近场15 km深度同震应力变化 Fig. 5 Coseismic stress changes in near field at the depth of 15 km | |
3 库仑应力3.1 计算公式和参数地震引起的永久变形可以改变周围弹性介质的应力分布, 静态库仑应力通常被用来解释余震的分布[36-37]与探讨地震的触发机制[38]。静态库仑应力公式为
ΔCFS=Δτ+μ′Δσ,
其中, Δτ为地震引起的剪应力变化, 与断层滑动方向一致为正; μ′为等效摩擦系数, 一般取值为0.2~0.8[37, 39]; 对走滑型断层与未知方向的断层, μ′通常取0.4左右[38]。Δσ为地震引起的正应力变化, 拉张为正。当库仑应力ΔCFS为正时, 则表示本次地震的发生有利于断层的相对错动, 地震发生的可能性增加; 相反, 库仑应力为负值时, 则表示抑制断层的相对错动, 地震发生的可能性减小。总之, 库仑应力表示的是一次地震活动对断层发生错动的贡献, 这个贡献是正应力变化与剪应力变化这两个贡献的加权和。当正应力变化是拉张的时候, 有利于断层错动; 当剪应力变化为正时, 也有利于断层错动。
3.2 苏门答腊地区断层的库仑应力Cattin等[17]给出苏门答腊地区及其周边的断层的走向、倾角与滑动角数据。该地区大多数走滑断层都近直立, 如实皆断层(Sagaing Fault)、西安达曼断层(West Andaman Fault)的北段与苏门答腊大断层等; 逆冲断层向西倾斜, 如西安达曼断层南段。由于印度-澳大利亚板块与欧亚板块在此处的斜向聚敛, 该地区断层没有左旋分量, 走滑断层基本上只有右旋分量, 逆冲分量一般也都很小, 只有靠近海沟的西安达曼断层的南段的逆冲分量较大。
使用Cattin等[17]给出的断层数据, 分别在分层球对称地球模型与三维横向不均匀性地球模型下计算15 km深度断层上的正应力、剪应力和在等效摩擦系数为0.4时的库仑应力(图 6)。等效摩擦系数0.4被广泛使用[36, 38, 40], 这里先使用该值作为等效摩擦系数。严格来说, 计算断层上的库仑应力需要考虑背景应力场的影响[41], 但真实背景应力场的资料很难获取, 这里假设苏门答腊地震并没有对这些断层的属性发生很大的改变。
Fig. 6
 | Download: JPG larger image |
| (a)-(c)分层球对称地球模型; (d)-(f)三维横向不均匀性地球模型。(a) (d)正应力, 以拉张为正; (b) (e)剪应力, 沿断层滑动方向为正; (c) (f)库仑应力, 以促进破裂为正。NAR (North Andaman Rift), 北安达曼洋中脊; ATZ (Andaman Transform Zone), 安达曼转换断层带; CAR (Central Andaman Rift), 中安达曼洋中脊; WAF (West Andaman Fault), 西安达曼断层; SEU (Seuliman Fault), 苏里曼断层; SGS (Sumatra Fault System), 苏门答腊断层系; 图中黑色和绿色方框同图 5。 图 6 苏门答腊地区及其周边断层15 km深度应力变化 Fig. 6 Coulomb stress changes on the faults in Sumatra area and its vicinity at the depth of 15 km | |
使用三维横向不均匀性地球模型得到的正应力, 除在西安达曼断层的南段的西支为负值之外, 其他部分几乎都为正值(图 6(d)); 而剪应力在西安达曼断层的南段均为负值(图 6(e)), 其他部分与分层球对称地球模型基本相同。由三维横向不均匀性地球模型得到库仑应力在西安达曼断裂的南段为负值, 其他部分与分层球对称地球模型基本相同。
对比分层球对称地球模型与三维横向不均匀性地球模型库仑应力, 在西安达曼断裂的南段(锡默卢岛以西北)差异较大, 前者为较大的正值(部分超过1 MPa), 而后者为负值(部分绝对值超过1 MPa)。这个差异应该是三维横向不均均性地球模型和分层球对称地球模型差异导致了该区域的计算应力差异(图 5(g)-(i)); 在其他断层相差不大。整个区域只有西安达曼断裂的南段的差异比较明显。图 6(b)和6(e)分别是分层球对称地球模型和三维不均匀性地球模型下的剪应力, 二者在西安达曼断裂的南段差异比较明显, 而其他部分的差异很小。这显示了该区域逆冲断层对介质模型的敏感性较大。图 6(a)和6(d)是断层面上的正应力, 也是在西安达曼断层的南段有较大的差异。只与断层的走向和倾角有关, 与断层的滑动方向无关。而剪应力的计算, 不仅与断层的走向和倾角有关, 还与断层的滑动方向相关。而走滑型断层和正断层的正应力和剪应力的变化很小(图 6(a)和(d), 图 6(b)和(e))。
由于库仑应力由剪应力和正应力两部分组成, 考虑等效摩擦系数(0.2, 0.4, 0.6)的不同对库仑应力的影响(图 7)。随等效摩擦系数的增加, 正应力部分的影响增加, 库仑应力会在剪应力和正应力正负相反的断层上发生极性反转。分层球对称地球模型下断层上的正应力几乎都为正值(图 6(a)), 因此随着等效摩擦系数的增加, 只有在个别剪应力为负值且数值十分大的断层部分库仑应力为负值之外, 其他断层大多为正值。三维不均匀性地球模型下, 西安达曼断层的南段部分正应力为负值(图 6(d)), 剪应力为负值或较小的正值, 由此得到的库仑应力为正值。在等效摩擦系数较大时, 除西安达曼断层的南段外, 其他断层的库仑应力在数值和极性上都基本相同(图 7(c)和(f))。
Fig. 7
 | Download: JPG larger image |
| (a)-(c)分层球对称模型；(d)-(f)三维横向不均匀性地球模型；(a)(d)等效摩擦系数为0.2; (b)(e)等效摩擦系数为0.4; (c)(f)等效摩擦系数为0.6。 图 7 苏门答腊地区及其周边断层15 km深度不同等效摩擦系数的库仑应力 Fig. 7 Coulomb stress of different frictional coefficients on the faults in Sumatra area and its vicinity at the depth of 15 km | |
3.3 余震位置处的库仑应力后续地震与主震引起的应力变化有关[42-45], 库仑应力增加的区域余震也会相应增加。由于大部分余震的震源机制解未知, 一般计算最优破裂面上的库仑应力。但最优破裂面不仅与地震造成的应力降有关, 还与背景应力场有关[41]，而背景应力场难以获得。苏门答腊地区不仅板块汇聚的机制等问题复杂, 而且物性参数在垂向和横向上变化都很大, 其背景应力场也比较复杂, 简单假定背景应力场很不合理。因此我们选用已知震源机制的余震：苏门答腊地震发生后到2005年3月27日之间的地震作为2004年苏门答腊地震的余震(图 8(a))。余震的破裂方向主要受控于地质构造和已有断裂走向, 而不是简单地与同震应力变化场和区域应力场联系[46]。震源机制解所提供的断层面有两个方向, 此处选用与主震破裂面走向相近的一组方向作为余震的破裂方向。
Fig. 8
 | Download: JPG larger image |
| (a)震级大于3级的余震震源机制; (b)分层球对称地球模型; (c)三维横向不均匀性地球模型。 (a)中标数字的细线黑色方框与图 9对应; 图中黑色和绿色方框同图 5。 图 8 余震发生处的库仑应力 Fig. 8 Coulomb stress in the site of aftershock | |
图 8(b)、8(c)与图 9均显示余震处的库仑应力, 前者为余震在地表上的投影; 后者为分别垂直于三段断层面走向上的投影, 以断层破裂面在地表的投影为横坐标的零点, 向东为正。图 9中, 右上角圆圈大小代表余震震级; 横轴坐标零点为断层初露地表位置, 纵坐标为深度; 灰色实线为断层面, 断层面从南向北变陡。使用分层球对称地球模型计算余震处的库仑应力, 67.6%的余震处的库仑应力为正值(图 8(b)); 使用三维横向不均匀性地球模型则有72.3%的余震处的库仑应力为正值。这大于缪淼和朱守彪[18]计算得到的49.8%的余震落在库仑应力增加区域。需要注意的是, 这里与他们的计算方法和数据不一样。在他们的计算中, 余震处的断层参数为最优破裂方向, 并且使用的余震的数量也大于我们所使用的余震数量, 得到只有不到一半的余震落在库仑应力增加的区域, 这个结果与通常认为的主震对余震的触发有很大的区别。
Fig. 9
 | Download: JPG larger image |
| (a)与(d)北段, 对应图 8(a)黑色方框1; (b)与(e)中段, 对应图 8(a)黑色方框2; (c)与(f)南段, 对应图 8(a)黑色方框3。(a)-(c)分层球对称地球模型; (d)-(f)三维横向不均匀性地球模型。(b)和(e)中的绿色短线代表西安达曼断层北段在剖面上的大致投影位置; (c)和(f)中的绿色短线代表西安达曼断层南段在剖面上的投影位置。 图 9 余震处库仑应力垂向剖面图 Fig. 9 Sections of Coulomb stress in the site of aftershock | |
尼克巴群岛以东的西安达曼走滑断层的北段附近沿断层展布发生了大量的余震, 且绝大多数为深度都不大的走滑型地震(图 8(a)), 主要分布在12~30 km深度(图 9b), 大部分余震处的库仑应力为正值, 并且附近的主要断层为走滑断层, 且这些断层上的库仑应力也为正值(图 6(c)和(f)), 故这些余震的发生应受到了附近断层的影响较大。由Chlieh等[6]给出的断层滑动模型显示, 在7°~11°范围内以逆冲为主(图 3(b)), 并且具有较大的滑动位移(图 3(a)), 释放了该地区的走滑断层部分正应力, 计算得到这些走滑型断层上的正应力为正值, 使得断层更容易滑动, 地震更容易发生。处于该位置处的西安达曼走滑断层北段的库仑应力(图 6)为较大的正值, 也有利于余震的发生。
在西安达曼断层的南段也有大量的余震发生, 但使用分层球对称地球模型计算得到的库仑应力在深度较大的余震处的库仑应力为负值(图 8(b), 图 9(b)与(c)), 三维横向不均匀性地球模型得到这些余震处的库仑应力为正值(图 8(c), 图 9(e)与(f))。在断层破裂面范围内靠近海沟一侧的余震处的库仑应力大多数都是负值, 且库仑应力的绝对值超过1 MPa。
图 9显示, 在发震断层破裂面的南段, 余震主要发生在断层的下盘; 在中段与北段, 余震主要发生在断层的上盘。余震的展布特征与主震断层破裂基本一致, 这印证了Chlieh等[6]所采用的断层几何模型的正确性(该几何模型由Ammon等[2]给出)。计算结果显示, 在深度小于35 km靠近断层破裂面的余震处的库仑应力大多数为负值, 在上盘余震处的库仑应力大多数为正值(图 9(a)-(b), (d)-(e))。在断层面的底部(图 9(b)-(e), (e)-(f))发生了大量地震, 在分层球对称地球模型下余震处库仑应力部分为较小的负值, 而使用三维横向不均匀性地球模型则都为正值, 且大部分超过1 MPa。
2005年3月28日在锡默卢岛与尼亚斯岛之间发生了Mw8.6级地震, 这里探讨其与2004年苏门答腊地震的关系。使用Konca等[35]的断层几何参数与滑动角数据, 等效摩擦系数为0.4, 计算2005年3月28日地震断层面上的库仑应力(图 10)。使用分层球对称地球模型和三位横向不均匀性地球模型得到的结果差异不是很大。2005年3月地震的震中附近接近一般认为的地震触发阈值0.1 bar (10 kPa)[44], 与Nalbant等[47]结果一致, 但小于Pollitz等[41]计算的0.25 bar (25 kPa)。在2005年3月28日地震震中的北西方向, 库仑应力为正值, 超过10 kPa; 在震中以北区域为负值, 超过-10 kPa; 在震中以南区域库仑应力较小。由此2004年苏门答腊地震在2005年3月28日地震的库仑应力在震中区域为为正值, 且在断层破裂面范围有部分区域的库仑应力超过地震触发阈值, 这里从库仑应力的角度得到2004年苏门答腊地震触发了2005年3月28日尼亚斯地震。
Fig. 10
 | Download: JPG larger image |
图 10 2004年苏门答腊地震对2005年3月28日地震库仑应力 Fig. 10 Coulomb stress induced by 2004 Sumatra earthquake on the fault plane of 2005 Nias earthquake | |
3.4 在中国西南地区产生的库仑应力2004年苏门答腊地震对中国西南地区也产生了很大的影响, 地震活动性显著增强, 也导致了地下水位的变化[14-15]。总体上, 苏门答腊地震引起的同震应力变化在中国西南地区较小(图 4), 15 km深度上中国西南地区主要断层上的库仑应力也相应较小(图 11)。图 11中断层缩写对应: HFS, 喜马拉雅断裂系; JLF, 嘉黎断裂; NJF, 怒江断裂; LCJF, 澜沧江断裂; RRF, 红河断裂; XSHF, 鲜水河断裂; XJF, 小江断裂; LMF, 龙门山断裂; HYF, 海原断裂; ALTF, 阿尔金断裂。YSF, 玉树断裂。小江断裂与红河断裂的库仑应力为正值, 在红河断裂的南段库仑应力值较大, 超过2 kPa。阿尔金断裂、嘉黎断裂、海原断裂、鲜水河断裂、龙门山断裂、澜沧江断裂与怒江断裂的库仑应力为负值, 但都小于500 Pa。这些断层上的库仑应力远小于触发余震的阈值10 kPa[44], 因此静态库仑应力的作用可以忽略[15]。但是2004年苏门答腊地震之后川滇地区却发生了大量的地震[14-15], 并且这些地震基本都是浅震[48]这些地震可能不是由于苏门答腊地震静态应力变化引起的, 而可能是由动态应力变化触发的[15]。
Fig. 11
 | Download: JPG larger image |
图 11 对中国西南地区的库仑应力 Fig. 11 Coulomb stress on the faults in Southwest China | |
4 讨论库仑应力变化与余震关系问题是一个复杂的问题。库仑应力计算中, 涉及两个要素:一个是接收断层的走向倾角; 一个是该处应力张量决定的潜在滑动方向。前者虽然也难以精确确定, 但地球物理和地质资料为大致确定提供了一些判断依据, 因此接收断层上正应力的变化计算也比较容易; 后者计算则比较困难, 因为应力张量一般均为未知, 潜在滑动方向并不容易确定, 也决定了断层上沿潜在滑动方向上的剪应力变化计算困难。对潜在错动方向的估算途径包括:利用知道震源机制的余震来确定滑动方向, 但震源机制中两个节面哪一个是断层面需要其他资料辅助确定; 参考主震错动模型外推, 但大地震不同段落错动不同、其他地方断层参数与应力状态与主震断层也未必相同; 利用应力实测资料进行推断, 但应力实测数据少而浅, 许多还仅仅提供二维平面应力, 有局限性; 根据多种方法(包括应力解除、水压破裂、震源机制、新构造等)绘制的世界和区域应力分布图中主压应力的方位推测潜在滑动方向, 但这类图件往往也只给出平面应力插值和平滑了的结果, 难以反映局部应力方向的变化和大小的起伏, 而许多余震与主震震源机制并不同, 说明这种插值平滑的结果未必反映局部现实。人们还考虑将特定走向、倾角断层上库仑应力最大的断层面, 即所谓“最优破裂面”上的库仑应力与余震活动比较, 但最优破裂面的确定也涉及主压应力方向或潜在滑动方向的难题。此外, 库仑应力计算还涉及孔隙流体压力的大小和变化等问题, 目前也没有得到定量深入研究。因此, 库仑应力的计算存在很多不确定因素。这样计算出的库仑应力大小正负也难以和余震或后续地震活动建立更加确定的关系。
从分层球对称地区模型到三维不均匀性地球模型, 余震处的库仑应力为正值的比例增加4.7%, 增加的部分大多数都位于断层破裂面的底部(图 9(b)和(e), 图 9(c)和(f)), 这显示考虑地球介质模型的三维不均匀性的影响的重要性。这是因为我们采用了更加真实的地球介质模型。如果有更为精细的地球模型, 这一比例有可能还会有小幅的增加。但其前提是要保证有限断层滑动模型的正确性。
在讨论地球介质模型的影响时, 需要针对某个确定的地震进行讨论。实际发生的地震由于其震级、震源深度、震源机制及地震所处的位置的不同, 地球介质模型的影响也不相同。本文所讨论的2004年苏门答腊地震为俯冲带地震, 所处位置岩石物性在垂向和横向上变化都很大, 同震效应受到岩石物性的不均匀性的影响相当大。本文所采用的三维不均匀性地球模型, 地壳部分为Crust1.0, 地幔部分为GyPSuM, 二者的精度都是1°×1°, 在苏门答腊地区在笛卡儿坐标系下相当于110 km×110 km的分辨率, 相对于海沟、岛弧等环境还是相当粗糙。
5 结论本文使用并行有限元方法, 使用已验证的与实际资料符合很好的Chlieh等[6]给出的断层滑动模型, 分别采用分层球对称地球模型与三维横向不均匀性地球模型, 计算苏门答腊地震引起的同震应力变化, 然后计算苏门答腊及周边地区断层与余震处的库仑应力。通过对计算结果的分析, 得出以下认识:
苏门答腊地震断层附近地壳与地幔三维横向不均匀性较强, 在同震破裂面附近, 介质横向不均匀性对同震应力变化的影响很大。考虑与不考虑介质横向不均匀性可导致南北方向同震正应力的最大差异可达147%, 东西方向同震正应力的最大差异为243%, 水平剪切应的最大差异为110%, 且在断层破裂面范围内同震应力变化各分量的差异超过1 MPa。
由分层球对称地球模型与三维横向不均匀性地球模型两种模型分别计算得到的苏门答腊及其附近断层上的库仑应力大致相近, 在西安达曼断层南段有较大的差异, 前者为正值, 后者为负值, 三维横向不均匀性地球模型更接近实际情况, 因此更加合理。这种差异也说明在地震断层破裂面附近有必要使用三维横向不均匀性地球模型。使用这两种模型计算得到的在67.6%与72.3%的余震处的库仑应力为正值, 这说明主震与余震有明显的相关性。
西安达曼断裂的北段和南段处的库仑应力与二者附近发生的余震, 均显示余震的发生与库仑应力为正值的主要断层有很好的对应关系, 余震受到其附近主要断层的控制作用。
2004年苏门答腊地震对2005年3月28日Mw8.6级地震震中的库仑应力约10 kPa, 接近一般认为的地震触发阈值0.1 bar。但在中国西南地区地库仑应力较小, 受苏门答腊地震的影响较小。
本文目前所使用的三维横向不均匀性地球模型的精度为1°×1°, 其精度还是比较粗糙, 介质属性对同震效应的具体影响还有待于进一步的研究和讨论。
感谢Cattin等[17]提供苏门答腊地区及其周边的断层数据。
参考文献
| [1] | Lay T, Kanamori H, Ammon C J, et al. The great Sumatra-Andaman earthquake of 26 December 2004[J].Science, 2005, 308(5725):1127–1133.DOI:10.1126/science.1112250 |
| [2] | Ammon C J, Ji C, Thio H, et al. Rupture process of the 2004 Sumatra-Andaman earthquake[J].Science, 2005, 308(5725):1133–1139.DOI:10.1126/science.1112260 |
| [3] | Banerjee F F P R. The size and duration of the Sumatra-Andaman earthquake from far-field static offsets[J].Secience, 2005, 308:1769–1772.DOI:10.1126/science.1113746 |
| [4] | Ni S, Kanamori H, Helmberger D. Seismology:Energy radiation from the Sumatra earthquake[J].Nature, 2005, 434(7033):582.DOI:10.1038/434582a |
| [5] | Vigny C, Simons W J, Abu S, et al. Insight into the 2004 Sumatra-Andaman earthquake from GPS measurements in southeast Asia[J].Nature, 2005, 436(7048):201–206.DOI:10.1038/nature03937 |
| [6] | Chlieh M, Avouac J, Hjorleifsdottir V, et al. Coseismic slip and afterslip of the great Mw 9.15 Sumatra-Andaman earthquake of 2004[J].Bulletin of the Seismological Society of America, 2007, 97(1A):S152–S173.DOI:10.1785/0120050631 |
| [7] | Rhie J, Dreger D, Rgmann R B U, et al. Slip of the 2004 Sumatra-Andaman Earthquake from joint inversion of long-period global seismic waveforms and GPS static offsets[J].Bulletin of the Seismological Society of America, 2007, 97(1A):S115–S127.DOI:10.1785/0120050620 |
| [8] | Shearer P, Rgmann R B U. Lessons learned from the 2004 Sumatra-Andaman megathrust rupture[J].Annual Review of Earth and Planetary Sciences, 2010, 38(1):103.DOI:10.1146/annurev-earth-040809-152537 |
| [9] | Bird P. An updated digital model of plate boundaries[J].Geochemistry, Geophysics, Geosystems, 2003, 4(3):1027. |
| [10] | Subarya C, Chlieh M, Prawirodirdjo L, et al. Plate-boundary deformation associated with the great Sumatra-Andaman earthquake[J].Nature, 2006, 440(7080):46–51.DOI:10.1038/nature04522 |
| [11] | Bock Y, Prawirodirdjo L, Genrich J F, et al. Crustal motion in Indonesia from global positioning system measurements[J].Journal of Geophysical Research:Solid Earth (1978-2012), 2003, 108(B8):2367.DOI:10.1029/2001JB000324 |
| [12] | Muller R D, Sdrolias M, Gaina C, et al. Age, spreading rates, and spreading asymmetry of the world's ocean crust[J].Geochemistry, Geophysics, Geosystems, 2008, 9(4). |
| [13] | 付广裕, 孙文科. 2004年苏门答腊地震引起的远场形变[J].大地测量与地球动力学, 2008, 28(2):1–7. |
| [14] | 张国民, 张晓东, 刘杰, 等. 印尼苏门答腊8.7级大震对中国陆区的影响[J].地震, 2006, 25(4):15–25. |
| [15] | Lei X, Xie C, Fu B. Remotely triggered seismicity in Yunnan, southwestern China, following the 2004 Mw9.3 Sumatra earthquake[J].Journal of Geophysical Research:Solid Earth (1978-2012), 2011, 116:B8303.DOI:10.1029/2011JB008245 |
| [16] | 马寅生, 史大年, 安美建, 等. 苏门答腊地震对柴达木地方震的触发作用[J].地质力学学报, 2005, 11(2):110–116. |
| [17] | Cattin R, Chamot-Rooke N, Pubellier M, et al. Stress change and effective friction coefficient along the Sumatra-Andaman-Sagaing fault system after the 26 December 2004(Mw=9.2) and the 28 March 2005(Mw=8.7) earthquakes[J].Geochemistry, Geophysics, Geosystems, 2009, 10:Q3011. |
| [18] | 缪淼, 朱守彪. 俯冲带上特大地震静态库仑应力变化对后续余震触发效果的研究[J].地球物理学报, 2012, 55(9):2982–2993. |
| [19] | Pollitz F F, Rgmann R B U, Banerjee P. Post-seismic relaxation following the great 2004 Sumatra-Andaman earthquake on a compressible self-gravitating Earth[J].Geophysical Journal International, 2006, 167(1):397–420.DOI:10.1111/gji.2006.167.issue-1 |
| [20] | Okada Y. Internal deformation due to shear and tensile faults in a half-space[J].Bulletin of the Seismological Society of America, 1992, 82(2):1018–1040. |
| [21] | Wang R, Mart I N F L, Roth F. Computation of deformation induced by earthquakes in a multi-layered elastic crust:FORTRAN programs EDGRN/EDCMP[J].Computers and Geosciences, 2003, 29(2):195–207.DOI:10.1016/S0098-3004(02)00111-5 |
| [22] | Fu G, Sun W. Effects of lateral inhomogeneity in a spherical earth on gravity earth tides[J].Journal of Geophysical Research:Solid Earth (1978-2012), 2007, 112:B6409.DOI:10.1029/2006JB004512 |
| [23] | 张贝, 张怀, 石耀霖. 有限元模拟弹性位错的等效体力方法[J].地球物理学报, 2015, 58(5):1666–1674. |
| [24] | Bangerth W, Rannacher R. Adaptive finite element methods for differential equations[M].Birkhauser: 2013. |
| [25] | Antonioli A, Piersanti A, Spada G. Stress diffusion following large strike-slip earthquakes:a comparison between spherical and flat-earth models[J].Geophysical Journal International, 1998, 133(1):85–90.DOI:10.1046/j.1365-246X.1998.1331490.x |
| [26] | Sun W, Okubo S. Effects of earth's spherical curvature and radial heterogeneity in dislocation studies:for a point dislocation[J].Geophysical research letters, 2002, 29(12):41–46. |
| [27] | Dong J, Sun W, Zhou X, et al. Effects of Earth's layered structure, gravity and curvature on coseismic deformation[J].Geophysical Journal International, 2014, 199(3):1442–1451.DOI:10.1093/gji/ggu342 |
| [28] | Qu W, Zhang B, Shi Y. Equivalent body force finite elements method and 3-D earth model applied in 2004 Sumatra earthquake[J].AGU Fall Meeting, 2015:NH51D–1921. |
| [29] | Dziewonski A M, Anderson D L. Preliminary reference earth model[J].Physics of the earth and planetary interiors, 1981, 25(4):297–356.DOI:10.1016/0031-9201(81)90046-7 |
| [30] | Laske G, Masters G, Ma Z, et al. Update on CRUST1.0-A 1-degree global model of earth's crust[J].Geophys Res Abstracts, 2013:EGU2013–2658. |
| [31] | Simmons N A, Forte A M, Boschi L, et al. GyPSuM:A joint tomographic model of mantle density and seismic wave speeds[J].Journal of Geophysical Research:Solid Earth (1978-2012), 2010, 115:B12310.DOI:10.1029/2010JB007631 |
| [32] | 瞿武林, 张贝, 黄禄渊, 等. 2004年苏门答腊地震的几个断层滑动模型的全球同震位移对比[J].地球物理学报, 2016, 59(8):2843–2858. |
| [33] | Sevilgen V, Stein R S, Pollitz F F. Stress imparted by the great 2004 Sumatra earthquake shut down transforms and activated rifts up to 400 km away in the Andaman Sea[J].Proceedings of the National Academy of Sciences, 2012, 109(38):15152–15156.DOI:10.1073/pnas.1208799109 |
| [34] | Mccaffrey R. The tectonic framework of the sumatran subduction zone[J].Annu Rev Earth Planet Sci, 2009, 37:345–366.DOI:10.1146/annurev.earth.031208.100212 |
| [35] | Konca A O, Hjorleifsdottir V, Song T A, et al. Rupture kinematics of the 2005 Mw8.6 Nias-Simeulue earthquake from the joint inversion of seismic and geodetic data[J].Bulletin of the Seismological Society of America, 2007, 97(1A):S307–S322.DOI:10.1785/0120050632 |
| [36] | King G C, Stein R S, Lin J. Static stress changes and the triggering of earthquakes[J].Bulletin of the Seismological Society of America, 1994, 84(3):935–953. |
| [37] | Stein R S. The role of stress transfer in earthquake occurrence[J].Nature, 1999, 402(6762):605–609.DOI:10.1038/45144 |
| [38] | Sumy D F, Cochran E S, Keranen K M, et al. Observations of static Coulomb stress triggering of the November 2011 M5.7 Oklahoma earthquake sequence[J].Journal of Geophysical Research:Solid Earth, 2014, 119(3):1904–1923.DOI:10.1002/2013JB010612 |
| [39] | Cotton F, Coutant O. Dynamic stress variations due to shear faults in a plane-layered medium[J].Geophysical Journal International, 1997, 128(3):676–688.DOI:10.1111/gji.1997.128.issue-3 |
| [40] | Pollitz F F, Banerjee P, Rgmann R B U, et al. Stress changes along the Sunda trench following the 26 December 2004 Sumatra-Andaman and 28 March 2005 Nias earthquakes[J].Geophysical research letters, 2006, 33:L6309.DOI:10.1029/2005GL024558 |
| [41] | 石耀霖, 曹建玲. 库仑应力计算及应用过程中若干问题的讨论:以汶川地震为例[J].地球物理学报, 2010, 53(1):102–110. |
| [42] | Das S, Scholz C H. Why large earthquakes do not nucleate at shallow depths[J]. 1983, 305:621-623. |
| [43] | Stein R S, Lisowski M. The 1979 Homestead valley earthquake sequence, California:control of aftershocks and postseismic deformation[J].Journal of Geophysical Research:Solid Earth, 1983, 88(B8):6477–6490.DOI:10.1029/JB088iB08p06477 |
| [44] | Harris R A. Introduction to special section:Stress triggers, stress shadows, and implications for seismic hazard[J].Journal of Geophysical Research:Solid Earth, 1998, 103(B10):24347–24358.DOI:10.1029/98JB01576 |
| [45] | Steacy S, Gomberg J, Cocco M. Introduction to special section:Stress transfer, earthquake triggering, and time-dependent seismic hazard[J].Journal of Geophysical Research:Solid Earth, 2005, 110:B05S01. |
| [46] | Mccloskey J, Nalbant S U L S, Steacy S, et al. Structural constraints on the spatial distribution of aftershocks[J].Geophysical research letters, 2003, 30(12):1610. |
| [47] | Nalbant S S, Steacy S, Sieh K, et al. Seismology:Earthquake risk on the Sunda trench[J].Nature, 2005, 435(7043):756–757.DOI:10.1038/nature435756a |
| [48] | Peng Y, Zhou S, Zhuang J, et al. An approach to detect the abnormal seismicity increase in Southwestern China triggered co-seismically by 2004 Sumatra Mw 9.2 earthquake[J].Geophysical Journal International, 2012, 189(3):1734–1740.DOI:10.1111/gji.2012.189.issue-3 |
