地球物理学报  2017, Vol. 60 Issue (7): 2641-2651   PDF    
2016年11月13日新西兰凯库拉(Kaikoura)MW7.8地震同震位移和应力数值分析
程惠红,张贝,张怀,石耀霖    
中国科学院计算地球动力学重点实验室, 中国科学院大学, 北京 100049
摘要: 2016年11月13日新西兰南岛北端凯库拉(Kaikoura)发生了MW7.8大地震,造成了强烈的地表变形并引发大面积滑坡和海啸的发生.基于美国地质调查局(USGS)断层滑动模型,建立全球同震横向不均匀并行椭球型地球模型,计算了此次新西兰凯库拉大地震产生的同震形变和应力及库仑应力变化.初步计算结果表明:新西兰凯库拉MW7.8地震造成断层上盘东北向抬升,下盘西南俯冲;引起发震区域同震位移较大,从凯库拉到坎贝拉(Campbell)以及首都惠灵顿(Wellington)整体上东北向抬升,最大同震水平位移1.2 m,垂直位移1.1 m.此次大地震释放了发震断层上积累的压应力,但增加了发震断层两端的挤压力;同时,同震剪应力变化增加了NE-SW向断层发生右旋滑动的危险性;采用此次地震发震断层参数计算得出的最大库仑应力变化增加区域集中在发震断层两端,可达到MPa量级.当分别采用新西兰北岛Awatere断裂系和南岛Wellington断裂系参数计算库仑应力变化时,发现新西兰北岛和南岛震中以南区域的库仑应力均增加,可触发部分余震的发生.
关键词: 新西兰地震      同震变形      库仑应力变化      有限元模拟     
Calculation of the co-seismic deformation and stress changes of the Kaikoura MW7.8 earthquake, Nov 13, 2016
CHENG Hui-Hong, ZHANG Bei, ZHANG Huai, SHI Yao-Lin    
Key Laboratory of Computational Geodynamics of Chinese Academy of Sciences, University of Chinese Academy of Sciences, Beijing 100049, China
Abstract: On 13 November, 2016, an MW7.8 earthquake occurred around Kaikoura, northern South Island, New Zealand. The MW7.8 earthquake caused strong surface deformation, massive landslides and tsunami. Based on the fault slip model released by United States Geological Survey, the co-seismic deformation and stress changes by Kaikoura MW7.8 earthquake were computed with the global heterogeneous ellipsoid Earth model and high precision topography. The preliminary result show that the footwall of the MW7.8 earthquake uplifted to the northeast, while the hanging wall underthrusting southwest. The co-seismic deformation triggered by this earthquake up to several center meters from Kaikoura to Campbell and the capital city of Wellington. The maximum co-seismic horizontal displacement is about 1.2 m whereas the vertical is about 1.1 m. Although the accumulated strain along seismogenic faults has been released by this earthquake. The extrusion force at the both ends of the seismogenic faults also increased due to co-seismic stress changes that the maximum value of the co-seismic stress changes reaches the order of MPa. The Coulomb failure stress changes are also up to MPa that concentrates in the vicinity of epicenter. At the same time, the risk for NE-SW dextral slip faults increases due to the co-seismic shear stress. We also calculated the Coulomb Failure Stress changes in the North Island and South Island with local faults systems respectively. Our results show that Coulomb failure stress changes are positive in both two regions, which means subsequent earthquakes are possible.
Key words: Kaikoura earthquake      Co-seismic deformation      Coulomb Failure Stress      Finite element method     
1 引言

2016年11月13日新西兰南岛北端凯库拉(Kaikoura)发生了MW7.8大地震,震中(经度173.1°E;纬度42.8°S)距离基督城(Christchurch)约95 km(https://en.wikipedia.org/),地震向NE方向扩展至库克(Cook)海峡,震源深度25.0 km,震源机制在南部以逆冲断层为主,北部为右旋走滑为主(http://earthquake.usgs.gov/).此次地震浅海破裂引起海啸(在凯库拉附近海啸波浪幅度达到1.5 m高,首都惠灵顿(Wellington)波浪高度为0.5 m左右),且导致沿海的山体发生大面积滑坡,是自1931年2月3日Hawke港湾M7.8地震之后新西兰发生的最大地震(http://www.iris.edu/[2016-11-13]).MW7.8大地震发生后余震不息且沿着滑动破裂面两侧分布,截止到11月19日,已监测到3000多个余震,400多个M4.0级以上的余震(http://earthquake.usgs.gov/[2016-11-19]),见图 1.新西兰位于太平洋板块和澳大利亚板块的交界位置,地质背景复杂,是构造运动和地震活动极为活跃的地区.大平洋板块在新西兰北岛和南岛北端向澳大利亚板块之下斜俯冲,而在南岛南端则是澳大利亚板块向太平洋板块之下俯冲(石耀霖和范桃园, 2001).南岛北部的陆壳俯冲带被大洋板块拖曳,冷而高密度的大洋俯冲板块携带大陆组分的Chatham隆起及附近窄条低密度陆壳岩石俯冲到新西兰南岛北段之下,复杂的板块边界条件导致诸多地震发生(石耀霖和范桃园, 2002),例如,2010年9月3日和2011年2月21日基督城附近相继发生了M7.1和M6.3级地震.凯库拉MW7.8地震的发生距离这两次地震约100 km,位于大洋俯冲带沿走向过渡为大陆俯冲-碰撞带(李忠海和石耀霖, 2016),该区域发育了一系列近北东走向的断裂体系(以右旋走滑为主)和次级连接断层(Langridge et al., 2003).此次MW7.8地震主震持续约2 min,从Hope断层南部首先破裂(石许华,个人通讯),进而向北东方向传递并引发一系列其他断层上的地震及地表破裂,破坏长度可达200 km(新西兰地质学和核科学研究所Institute of Geological and Nuclear Sciences (GNS) https://www.gns.cri.nz/[2016-11-21]).

图 1 新西兰及其邻近地质构造背景和地形高程 红-白色震源球表示此次MW7.8地震震源位置(下方来源于USGS,上方来源于CMT),黑-白色震源球为1976-01-01—2016-11-13的震级>6.0的历史地震分布(http://www.globalcmt.org[2016-11-13]),蓝-白色震源球为震后震级>5.0的余震;紫色圆点为MW7.8地震震后震级>2.5的余震和区域小震分布(截止2017年01月10日(http://earthquake.usgs.gov/).红色断层为俯冲带,蓝线表示新西兰南岛北部Awatere断裂系和北岛近N-E向Wellington断裂系.左上图为新西兰南岛北部大洋俯冲带沿走向过渡为大陆俯冲-碰撞带示意图 Fig. 1 The geological background and elevation in New Zealand and its adjacent area The focal sphere with red-white stands for the main MW7.8 earthquake and focal spheres with blue-white stand for aftershocks with magnitude >5.0. The purple squares stand for aftershocks in the figure before Jan 10, 2017 (http://earthquake.usgs.gov/); Focal spheres with black-white mean historical earthquakes with magnitude larger than 6.0 since 1976 (http://www.globalcmt.org). Red faults mean subduction zone; Blue lines represent Awatere fault system in the South Island and Wellington fault system in the Northern Island. Schematic diagram in the upper left subgraph shows transition zone from oceanic subduction to continental collision around New Zealand.

地应力长期累积增加而超过岩体强度时往往会发生地震.一次大地震的演化过程一般分成孕育、发生和震后调整三个阶段,然后进入下一个大地震的孕育阶段.其中,地震发生时产生的同震形变和应力变化最为迅速和巨大,且对区域后续地应力演变和地震活动有重要影响.因此,通过计算地震发生时同震位移和应力及库仑应力变化对于后续地震危险性评估和地震触发等问题的研究有重要的意义.本文应用全球同震横向不均匀并行椭球型地球模型计算了新西兰MW7.8大地震同震位移和应力变化,同时计算了地震引起周围区域库仑应力变化,以期估计大地震后新西兰地区地震危险性和应力环境的影响.

2 计算方法与断层模型

地球作为一个复杂的物理系统,在受到力的作用时短期内可表现为弹性体.因此,可将地球作为弹性体处理来研究分析地震发生短期内的同震效应.据此,20世纪50年代以来研究学者将位错理论引入地震学(Steketee, 1958),并发展了这套理论(Chinnery, 1961, 1963; Press, 1965; Okada, 1985).Okada给出了计算平面半无限空间地球的解析公式(Okada, 1985, 1992);汪荣江等给出了计算平面分层地球的解析公式(Wang et al., 2003);孙文科等给出计算球状分层地球的解析公式(Sun et al., 2009).然而,真实的地球既不是半无限空间,也不是简单的同心球层结构,而是一个地表有地形起伏,地壳底部有Moho面起伏,具有椭率的球体,且在同一圈层内物质的密度、弹性模量等参数也不尽相同.虽然解析解有着计算快捷简便的优点,但是只适用于比较简单的地球模型,无法计算上述这些因素对结果造成的影响.Ben-Menahem和Israel(1970)通过均质无自重球形模型计算,得出地球的曲率影响在震源矩20°以内可忽略不计,而对于一个震级8.5的浅源地震在一定的纬度上最大有30%的误差.Sun和Okubo(2002)认为地球曲率和分层对同震垂直位移的影响不可忽略:假设地球自重影响小于地球曲率影响的情况下,随着震源深度的增加地球曲率的影响增大,且引张震源的曲率效应比剪切震源的大.林晓光和孙文科(2014)采用二维有限元方法计算了2011年日本MW9.0大地震得出地形效应对计算结果影响为34%.然而,采用高性能数值计算方法可以很好地避免这些方面的困难,可以灵巧地将地球椭率、地形和Moho起伏等影响因素融入计算模型,可得出更加切合实际的计算结果,进而进行详细的分析研究.

张贝等(2015)解决了传统有限元方法数值计算同震位错断层不连续问题,提出了采用弹性位错等效体力有限元方法,大幅降低了数值建模的难度,减少了整个计算的耗时.在实际使用中,计算等效体力项所需的代码量小,同时,添加体力项也是各种主流有限元计算程序都支持的功能,对计算网格更是没有任何额外的要求,并可以处理任意复杂形状的断层,十分方便.基于张贝等(2015)等效体力方法的基础上,我们将全球球型分层计算模型进一步改进,添加地形起伏、地球椭率和介质的横向不均匀性.

在构建含椭率和地形起伏计算网格模型时,首先生成没有地形的球形地球,然后选取合适的地形数据(例如,Etopo1),计算出地表每个节点需要挪动的距离,将这些节点挪动到目标位置处.为了避免挪动节点造成的网格畸变,对地球内部的节点施加适当的挪动量.具体方法是将地表节点的挪动量看做位移,将这个位移作为地球模型的边界条件,而地表以下的节点按比例移动,以地球模型为求解区域,求解一个Poisson方程(见式(1)),得到地球内部各节点需要的位移量.Poisson方程解的光滑性保证了使用这种方式添加地形不会造成太大的网格畸变.

(1)

对于横向不均匀性问题处理时不采用移动网格的方式,而是在组建刚度矩阵的时候,根据所采用的介质参数数据,针对不同的Gauss积分点位置选取对应的参数进行计算.

图 2a给出了采用等效体力处理断层位错的凯库拉MW7.8地震同震计算模型.模型在使用结构化网格的基础上,断层附近采取自适应加密网格方法,且将地球地表高程和Moho起伏及横向不均匀性效应加入到并行计算模型上.边界条件的设置为:上表面为自由表面,下表面为核幔边界,处理为弹簧边界条件,弹簧刚度取决于地核和地幔物质密度差.该模型优点:① 采用网格自适应加密技术,自动将断层处网格加密,在保证计算精度的同时合理控制计算量;② 使用完整的地球参与计算,无需设定侧边界条件,避免不合理的边条件带来的误差;③ 等效体力代替位错源,避免在计算网格中显式地包含断层,使计算中最耗时的建模过程大大简化;④ 可灵活设置断层形态和模型参数; ⑤ 使用球坐标可以避免直角坐标计算应力应变的系统误差;⑥ 将地表高程和Moho面起伏及地球的椭率融入到计算模型中.

图 2 全球同震横向不均匀椭球型地球模型和新西兰凯库拉MW7.8断层滑动模型 (a)计算模型和网格;模型灰色系色标表示椭球型地球高程,彩色标显示了计算中采用的不均匀密度值;在发震区网格加密;(b)美国地质调查局提供的滑动模型(http://earthquake.usgs.gov);(c)数值计算中采用等效体力方法拟合的滑动模型. Fig. 2 The computing ellipsoidal Earth model with topography, the curvature of Earth and heterogeneity of the earth interior and heterogeneity of the earth interior (a) The ellipsoidal Earth model. Right sub-picture shows encryption grids with adaptive mesh method; (b) The slip model from USGS (http://earthquake.usgs.gov); (c) The slid model in our numerical model with equivalent body force approach.

新西兰凯库拉MW7.8地震之后,美国地质调查局(http://earthquake.usgs.gov[2016-11-13])迅速给出断层滑动模型图 2b,发震断层参数:走向219°,倾角38°,滑动128°;地震矩7.04×1020N·m,滑动面上最大滑动位移至少4 m.据此模型,我们采用等效体力方法给出模型在断层面上对应的等效体力,见图 2c.计算网格53万个单元,65万个节点.

3 计算结果 3.1 新西兰凯库拉MW7.8地震产生的同震位移场

图 3给出了新西兰凯库拉MW7.8大地震产生的地表同震位移模拟计算结果.由于此次MW7.8地震向北东方向扩展且在地震开始之后释放了最大能量,因此,同震位移变化明显的区域从凯库拉到坎贝拉(Campbell)以及首都惠灵顿,水平位移整体上是东北向抬升.但是,最大同震位移并未出现在震中,而是沿着断层在震中的东北50~60 km处产生较大同震位移.在此次数值模拟计算中,我们选取USGS的断层滑动模型,计算得出最大同震水平位移约为1.2 m;地震造成南岛区域东西方向上的水平位移相对南北方向较大,新西兰南岛基督城和惠灵顿之间区域东向和北向运动至少3 cm以上:新西兰南岛基督城及邻近区域东向位移至少5 cm;凯库拉以北地区东向运动则至少15 cm;首都惠灵顿东向运动大于30 cm;新西兰东部海域区域则西向运动约5 cm.另外,该地震造成沿海海岸线和海底有一定的抬升,据滑动模型计算得出的最大位移可达2 m.由图 3(d, e)地震地表同震水平位移矢量图,可以看出此次新西兰MW7.8地震是一个逆冲地震,上盘向东北抬升,下盘则向西南俯冲;地震造成周围100 km范围内的区域同震位移0.15 m以上,引起垂直于发震断层500 km范围内的新西兰周围海域变化约5 cm.

图 3 新西兰MW7.8地震引起其及邻近区域同震位移变化 东西向以东为正,南北向以北为正,垂向上以上为正. (a)—(c)分别表示同震水平位移南北向和东西向及垂直位移. (d)—(e)地震引起的地表水平位移矢量图. Fig. 3 Co-seismic displacements caused by the MW7.8 earthquake with positive in East, North and Upward (a)—(c) show co-seismic displacements along N-S horizontal and E-W horizontal directions, and vertical direction, respectively; (d)—(e) present vector maps of surface co-seismic horizontal displacement.

凯库拉MW7.8地震后,新西兰地质学和核科学研究所(GNS)给出了整个新西兰地区的GPS观测结果(http://info.geonet.org.nz/display/quake/2016/11/18/GPS+allowed+rapid+detection+of+land+movements+due+to+M7.8+earthquake),其观测显示此次地震造成了整个国家移动,凯库拉向北东方向移动了近70 cm,凯库拉地区约为2 m,首都惠灵顿和Kapiti地区整体向北移动2~6 cm,南岛部分西海海岸向东移动约10 cm,而北岛北部和南岛南部仅仅移动了几个mm.进一步同GPS观测结果相对比,我们的计算结果可以较好地吻合上,同时,也验证了我们计算程序的可靠性.

3.2 新西兰凯库拉MW7.8地震同震应力场变化

地震的发生是一个区域应力调整和应变能释放的过程.为了进一步探讨新西兰凯库拉MW7.8地震对周围区域应力场的影响,图 4给出了此次地震造成的震源深度处同震应力场变化.因垂直方向应力变化较小,在此仅绘出水平应力变化图.可以看出MW7.8地震造成同震应力最大值在MPa量级,从凯库拉到坎贝拉以及首都惠灵顿同震应力变化均比较大,量级大于10.0 kPa.地震造成发震断层两侧的东西向正应力σφφ和南北向正应力σθθ释放,压应力显著降低(弹性力学定义张应力为正,正应力增加即表示压应力减小);而使得断层两端的压应力增加.同时,此次地震引起的水平剪应力σθφ变化明显,余震均落在这些区域;正的剪应力变化将会引起南北向断层发生右旋滑动和东西向断层左旋滑动.

图 4 新西兰MW7.8地震引起其及邻近区域同震应力变化 (a)—(c)分别为东西向水平正应力σφφ、南北向水平正应力σθθ和水平剪应力σθφ变化. Fig. 4 Horizontal stress changes caused by the MW7.8 earthquake (a)—(c) E-W normal stress σφφ, N-S normal stress σθθ, and horizontal shear stress σθφ.

近年来的研究表明,在背景应力场不清楚的情况下,可用主震产生的库仑应力(Coulomb Failure Stress, CFS)变化来判断余震触发和地震危险性(石耀霖和曹建玲, 2010).Jacques等(1996)研究1978年北非Asal裂谷地震群,计算了库仑应力变化值,解析了为期六个星期地震群分布特征,且得出+0.03 MPa的库仑应力变化是可以触发地震发生的结论.Stein和Lisowski等(1983)考察了1979年加州地震造成沿破裂面及断层周边上应力的变化.Lomnitz(1996)研究了90年世界7级以上的浅源地震,并得出地震间的触发发生的距离在300~1000 km.King等(1994)计算了1992年Lers7.4级地震造成破裂面及周围断层上的库仑应力变化.石耀霖(2001)提出在大量的小地震、慢地震和无震滑动对应力调整、应力触发和应力大小的影响.邵志刚等(2010)计算了汶川地震引起的库仑应力动态演化过程,分析了对周边断层活动的影响.金笔凯等(2011)以2010年新西兰Canterbury MW7.1主震和2011年Christchurch MW6.3余震为例,应用汪荣江(Wang et al., 2003)PSGRN/PSCMP程序计算了同震应力变化,用库仑应力变化判断准则研究主震对余震的触发作用.

库仑应力变化表达式为

(2)

式中Δσn是正应力变化,Δτ是沿断层面剪应力沿滑动方向的变化,μ为断层面内摩擦系数,Δp为孔隙压力变化.本文中应力正负遵循弹性力学定义.公式(2) 表明了当断层面上的剪应力和孔隙压力增加、正应力减小时,库仑破裂应力为正,说明应力调整增加了断层面滑动的危险性;相反,则阻碍了断层破裂滑动.设地震发生应力场变化了Δσij,由震源机制解可得到发震断层几何参数(走向、倾角和滑动角),则可求出断层面法向量(n1, n2, n3)及正应力和剪应力变化量:

(3)

其中Δσij是断层面上正应力变化,Δτ是沿断层面剪应力沿滑动方向的变化,ni是断层面的法向量.经张量计算可以求出地震发生造成断层面上的库仑应力变化ΔCFS.库仑应力变化考虑了正应力、剪应力和孔隙压力变化的力学效应.通过库仑应力变化计算,可以了解哪些地方、哪种走向倾角和哪种性质的断层相对变得更危险或更安全.本文在此暂不考虑孔隙水压的作用.由公式(2) 可知,摩擦系数μ对库仑应力变化的影响较大,据此,我们分别采取μ值为0.2、0.6和0.8来讨论摩擦系数对计算结果的影响.另外,水的化学效应,例如水可以使断层弱化、强度减低,由于目前还缺乏摩擦系数变化Δμ与水的定量关系表达式,在数值模拟中还没能纳入.

图 5(ac)给出了采用不同摩擦系数和此次MW7.8震源机制参数情况下周围区域库仑应力变化,可以看出MW7.8地震释放了一定的能量,最大库仑应力变化量级在MPa范围.另外,当采用凯库拉MW7.8地震发震断层参数计算时,新西兰南岛西北包括纳尔逊(Nelson)区域及西北侧和震区东侧海岸线附近、东侧太平洋、首都惠灵顿的库仑应力变化为负值,使得这些区域走向和倾角接近发震断层的断层变得相对安全,但同时也将导致南岛震中以南的区域和北岛西部和东部地区这种类型的断裂危险性增加.同时,进一步对比图 5(ac),可以看出采用不同的摩擦系数计算得出的库仑应力变化正负范围有一定的影响,例如,摩擦系数的增加,北岛北部库仑应力变化负值增加,区域类似于发震断层的安全性提高;南岛南部类似发震断层的断裂库仑应力变化值也相对减小一些.在图 5a子图给出了5个M>5.0(http://www.globalcmt.org[2016-11-17])和330个M>2.5(http://earthquake.usgs.gov/[2016-11-17])在库仑应力变化图的分布情况.可以看出大部分的余震并未分布在库仑应力变化为正的区域.由公式(3) 可知,库仑应力变化与计算时候采用的断层参数有密切的关系,而震后野外地质考察认为此次MW7.8震有多条断层发生破裂(http://info.geonet.org.nz/display/quake/2016/11/18/Viewing+the+M7.8+Kaikoura+earthquake+from+space).邓园浩等(2017)对2016年3月2日苏门答腊MS7.8地震同震计算分析时发现有效摩擦系数的变化和区域构造应力场的耦合作用以及研究区域介质结构和断层的复杂性对库仑应力变化计算值有一定的影响.利用改进的库仑应力变化计算方法和最优破裂方法计算可以显示库仑应力触发理论更好地解释余震分布,因此,对于此次MW7.8地震的震后余震的库仑应力变化需要进一步结合与余震发生的断层参数来做详细分析.

图 5 新西兰MW7.8地震引起其及邻近区域库仑应力变化 (a)—(c)采用此次USGS确定的主震错动面(走向221.2°,倾角38.0°和滑动角145°)的断层参数,摩擦系数分别为0.6、0.2和0.8计算的库仑应力变化;(d)—(e)分别为采用新西兰北岛和南端断层参数(走向235°、倾角45°和滑动角-52°)和(走向57°、倾角60°和滑动角136°),摩擦系数为0.6计算的库仑应力变化.图a子图中红-白色震源球分别为USGS和哈佛大学CMT给出的凯库拉MW7.8地震震源位置,蓝-白色震源球为MW7.8地震后5个M>5.0的余震分布(来源http://www.globalcmt.org);白色圆点为330个M>2.5余震分布(来源http://earthquake.usgs.gov/). Fig. 5 The changes of Coulomb failure stress caused by the MW7.8 earthquake (a)—(c) show the ΔCFS caused by the MW7.8 earthquake using the USGS slip model with the friction coefficient 0.6、0.2 and 0.8; (d)—(e) show the ΔCFS caused by the MW7.8 earthquake using the parameters from mechanisms of historical earthquakes of North Island and South land, respectively. The focal sphere with red-white stands for the main MW7.8 earthquake and focal spheres with blue-white stand for 5 aftershocks with magnitude >5.0. The purple squares stand for 330 aftershocks in the figure (http://earthquake.usgs.gov/).

图 5(de)分别为结合北岛和南岛发生的历史地震和该区域断裂组参数来计算区域库仑应力变化.利用北岛的断层参数(走向235°、倾角45°和滑动角-52°)、南岛的断层参数(走向57°、倾角60°和滑动角136°)来分别计算库仑应力变化.由图可知,区域库仑应力变化分布差异较大,2016年新西兰凯库拉MW7.8地震引起整个北岛和南岛的同震库仑应力变化为正值,最大可达到MPa量级,但随着离震源位置距离的增加,库仑应力增量逐渐衰减,在新西兰南侧沿岸地区库仑应力增加为kPa量级.

4 结论与讨论

新西兰地处太平洋板块和印度洋板块之间的活动边缘上,新构造运动十分强烈,地震活动频繁.此次南岛北部凯库拉MW7.8地震,距2011年新西兰基督城(Christchurch)MW6.3地震约150 km,是新西兰东部地震活动性高的区域(Shi et al., 2017).早期的地质学者曾经根据构造运动和活断层的分布和走向对新西兰地区进行了构造规划和地震危险性区划,且预测了该地区在一百年内有多次破坏性地震的发生(Alder et al., 2016; 陈时军等,2003).长时间的构造应力积累使得大地震发生,地震的发生会引发一系列物理场变化,其中力学场的变化包括反映地球弹性的同震(Co-seismic)位移、应力的变化和反映地球黏弹性响应的震后(Post-seismic)位移、应力变化等(石耀霖,2013).这些物理场的变化,进一步对周围的介质产生作用.本文据美国地质调查局给出的此次地震滑动模型,搭建全球同震横向不均匀并行椭球型地球模型,对新西兰凯库拉MW7.8地震进行了同震位移和应力的计算.计算结果获得以下认识:(1) 新西兰凯库拉MW7.8地震是一个海洋板块俯冲大陆板块的逆冲兼一定走滑的地震,断层上盘东北向抬升,下盘西南俯冲.地震造成的同震位移较大,计算采用USGS滑动模型得出的最大水平位移1.2 m,相比地震造成东西向水平位移和垂直方向的位移比南北方向上水平位移较大;(2) 此次MW7.8地震释放了断层两侧一定的正应力,但断层两端的挤压力增加;造成同震应力最大值在MPa量级;正的剪应力变化将会引起南北向断层发生右旋滑动和东西向断层左旋滑动.由此次震源机制计算得出的库仑应力变化显示地震造成新西兰南岛震中以南地区的库仑应力变化为正,地震危险性依然较大;同时,也将导致北岛东西海岸区域的库仑应力增大;摩擦系数的减小将增大库仑应力变化为正的区域;(3) 分别以北岛和南岛断层发震机制来计算库仑应力变化,结果显示此次MW7.8地震造成北岛和南岛大部分区域库仑应力增加,量级在kPa范围.

本文采用全球同震横向不均匀并行椭球型地球模型,消除了地球曲率、地形和均质的影响.但采用的线弹性本构,只计算地震发生瞬间的同震响应,探讨震后较短时间内对地震活动性的影响,未考虑地震发生后黏弹性的松弛作用的黏弹效应.Casarotti等指出,震后松弛作用造成的应力变化不仅可以影响到近场断层,而且对远场的断层稳定性也会造成影响(Casarotti et al., 2001; Shao et al., 2016);并能影响到一定区域内的地震活动性等.这些研究表明震后应力对断层稳定性以及地震活动性有不可忽略的作用.因此,后续将在此模型的基础之上发展成为黏弹性位错模型,对平衡方程使用Laplace变换,将包含黏弹性本构关系的平衡方程转化成和弹性本构的平衡方程同样的形式,使用弹性位错的方法进行求解,再对得到的解使用Laplace逆变换,最终获得黏弹性位错问题的解.同时,断层形态、摩擦系数等参数有着较大的不确定性,在计算的震前构造应力场时较为简单,且没有考虑动态应力触发和流体作用等影响,在一定程度上这些影响因素对计算结果有着很大的影响,特别是余震的分布是否落在库仑应力变化为正的区域,是否库仑应力触发理论可更好地解释余震分布,这些方面今后将继续深入研究.此外,同震-震后位错数值计算模型主要采用断层破裂滑动模型来作为初始条件,采用不同研究机构给出的滑动模型得出的计算结果也不同,甚至相反.瞿武林等(2016)对2004年苏门答腊MW9.1~9.3大地震4个代表性断层滑动模型全球同震位移研究发现,在利用同震位移观测资料反演断层滑动模型时由于采用的介质模型不同而得出不同的断层滑动模型,且这些滑动模型差别颇大,特别是对于大地震远场位移分布的影响.因此,对于大地震的同震-震后的形变数值计算可选取不同研究机构的断层滑动模型,据全球同震位移观测同计算结果对比,进而挑选出一个更为合理的滑动模型来对研究区域做详细的分析.

致谢

感谢新加坡南洋理工大学石许华博士的交流和宝贵意见.感谢本文编辑和两位审稿专家对本文提出的建设性意见和建议.

参考文献
Alder S, Smith S A F, Scott J M. 2016. Fault-zone structure and weakening processes in basin-scale reverse faults:The Moonlight Fault Zone, South Island, New Zealand. Journal of Structural Geology, 91: 177-194. DOI:10.1016/j.jsg.2016.09.001
Ben-Menahem A, Israel M. 1970. Effects of major seismic events on the rotation of the Earth. Geophysical Journal International, 19(4): 367-393. DOI:10.1111/gji.1970.19.issue-4
Casarotti E, Piersanti A, Lucente F P, et al. 2001. Global postseismic stress diffusion and fault interaction at long distances. Earth and Planetary Science Letters, 191(1-2): 75-84. DOI:10.1016/S0012-821X(01)00404-6
Chinnery M A. 1961. The deformation of ground around surface faults. Bulletin of the Seismological Society of America, 51(3): 355-372.
Chinnery M A. 1963. The stress changes that accompany strike-slip faulting. Bulletin of the Seismological Society of America, 53(5): 921-932.
Deng Y H, Cheng H H, Zhang H, et al. 2017. Numerical study on the co-seismic deformation and stress changes of the MS7.8 Sumatra earthquake, March 2, 2016. Chinese Journal of Geophysics, 60(1): 174-186. DOI:10.6038/cjg20170115
Jacques E, King G C P, Tapponnier P, et al. 1996. Seismic activity triggered by stress changes after the 1978 events in the Asal Rift, Djibouti. Geophysical Research Letters, 23(18): 2481-2484. DOI:10.1029/96GL02261
Jin B K, Ni S D, Zhan Z W, et al. 2011. Coulomb stress change sensitivity due to variability in mainshock models and receiving fault parameters:a case study of the 2010-2011 Chritchurch, New Zealand Earthquakes.//The Chinese Geophysics 2011 (in Chinese). Beijing:Chinese Geophysical Society, 364.
King G C P, Stein R S, Lin J. 1994. Static stress changes and the triggering of earthquakes. Bulletin of the Seismological Society of America, 84(3): 935-953.
Langridge R J, Campbell J, Hill N, et al. 2003. Paleoseismology and slip rate of the Conway Segment of the Hope Fault at Greenburn stream, South Island, New Zealand. Annals of Geophysics, 46(5): 1119-1139.
Li Z H, Shi Y L. 2016. Constrains of 3-D plate geometry on the dynamics of continental deep subduction. Chinese Journal of Geophysics, 59(8): 2806-2817. DOI:10.6038/cjg20160808
Lin X G, Sun W K. 2014. Effects of topography and local geological structure on computing co-seismic deformation——a case study of the 2011 Japan Tohoku earthquake (MW9.0). Chinese Journal of Geophysics, 57(8): 2530-2540. DOI:10.6038/cjg20140814
Lomnitz C. 1996. Search of a worldwide catalog for earthquakes triggered at intermediate distances. Bulletin of the Seismological Society of America, 86(6): 293-298.
Okada Y. 1985. Surface deformation due to shear and tensile faults in a half-space. Bulletin of the Seismological Society of America, 75(4): 1135-1154.
Okada Y. 1992. Internal deformation due to shear and tensile faults in a half-space. Bulletin of the Seismological Society of America, 82(2): 1018-1040.
Press F, Briggs P. 1975. Chandler Wobble, earthquakes, rotation, and geomagnetic changes. Nature Letter, 256(5515): 270-273. DOI:10.1038/256270a0
Qu W L, Zhang B, Huang L Y, et al. 2016. Comparisons of global coseismic displacements from several fault slip models for the 2004 Sumatra earthquake. Chinese Journal of Geophysics, 59(8): 2843-2858d. DOI:10.6038/cjg20160811
Shao Z G, Zhou L Q, Jiang C S, et al. 2010. The impact of Wenchuan MS8.0 earthquake on the seismic activity of surrounding faults. Chinese J. Geophys, 53(8): 1784-1795. DOI:10.3969/j.issn.0001-5733.2010.08.004
Shao Z G, Zhan W, Zhang L P, et al. 2016. Analysis of the far-field co-seismic and post-seismic responses caused by the 2011 MW9.0 Tohoku-Oki Earthquake. Pure and Applied Geophysics, 173(2): 411-424. DOI:10.1007/s00024-015-1131-9
Shi X, Wang Y, Jing L Z, et al. 2017. How complex is the 2016 Mw7.8 Kaikoura earthquake, South Island, New Zealand?. Science Bulletin, in Press.
Shi Y L. 2001. Stress triggers and stress shadows:How to apply these concepts to earthquake prediction. Earthquake, 21(3): 1-7.
Shi Y L, Fan T Y. 2001. Analysis on maximum scale of continental crust sliver driven by subducting oceanic slab:a case study of south island of New Zealand and UHPM at Dabie, China. Chinese J. Geophys., 44(6): 754-760.
Shi Y L, Fan T Y. 2002. Continental sliver subduction driven by oceanic slab:Tectonic settings for large scale Ultra-high pressure Metamorphism. Acta Geologica Sinica, 76(1): 77-82.
Shi Y L, Cao J L. 2010. Some aspects in static stress change calculation-case study on Wenchuan earthquake. Chinese J. Geophys., 53(1): 102-110. DOI:10.3969/j.issn.0001-5733.2010.01.011
Steketee J A. 1958. On Volterra's dislocations in a semi-infinite elastic medium. Canadian Journal of Physics, 36(2): 192-205. DOI:10.1139/p58-024
Stein R S, Lisowski M. 1983. The 1979 Homestead valley earthquake sequence, California:control of aftershocks and postseismic deformation. J. Geophys. Res., 88(B8): 6477-6490. DOI:10.1029/JB088iB08p06477
Sun W K, Okubo S. 2002. Effects of the earth's spherical curvature and radial heterogeneity in dislocation studies-for a point dislocation. Geophys. Res. Lett., 29(12): 46-1-46-4.
Sun W K, Okubo S, Fu G Y, et al. 2009. General formulations of global co-seismic deformations caused by an arbitrary dislocation in a spherically symmetric earth model-Applicableto deformed earth surface and space-fixed point. Geophys. J. Int., 177(3): 817-833. DOI:10.1111/gji.2009.177.issue-3
Wang R J, Martín F L, Roth F. 2003. Computation of deformation induced by earthquakes in a multi-layered elastic crust-Fortran programs EDGRN/EDCMP. Comput. Geosci., 29(2): 195-207. DOI:10.1016/S0098-3004(02)00111-5
Zhang B, Zhang H, Shi Y L. 2015. Equivalent-bodyforce approach on modeling elastic dislocation problem using finite element method. Chinese Journal of Geophysics, 58(5): 1666-1674. DOI:10.6038/cjg20150518
陈时军, David, Harte, 等. 2003. 新西兰地区地震活动时空分布的多重分形特征研究. 地震学报, 25(3): 298–307.
邓园浩, 程惠红, 张怀, 等. 2017. 2016年3月2日苏门答腊MS7.8地震同震位移和应力场数值模拟研究. 地球物理学报, 60(1): 174–186. DOI:10.6038/cjg20170115
金笔凯, 倪四道, 詹中文等. 2011. 主震模型和接受断层参数的选取对库仑应力变化的影响: 2011年新西兰Christchurch地震研究. //中国地球物理2011. 北京: 中国地球物理学会, 364.
李忠海, 石耀霖. 2016. 三维板块几何形态对大陆深俯冲动力学的制约. 地球物理学报, 59(8): 2806–2817. DOI:10.6038/cjg20160808
林晓光, 孙文科. 2014. 地形效应和局部地质构造对计算同震形变的影响——以2011年日本东北大地震(Mw9. 0) 为例. 地球物理学报, 57(8): 2530–2540. DOI:10.6038/cjg20140814
瞿武林, 张贝, 黄禄渊, 等. 2016. 2004年苏门答腊地震的几个断层滑动模型的全球同震位移对比. 地球物理学报, 59(8): 2843–2858. DOI:10.6038/cjg20160811
邵志刚, 周龙泉, 蒋长胜, 等. 2010. 2008年汶川MS8.0地震对周边断层地震活动的影响. 地球物理学报, 53(8): 1784–1795. DOI:10.3969/j.issn.0001-5733.2010.08.004
石耀霖. 2001. 关于应力触发和应力影概念在地震预报中应用的一些思考. 地震, 21(3): 1–7.
石耀霖, 范桃园. 2001. 大洋岩石层拖曳窄条陆壳俯冲的极限尺度分析——以新西兰南岛和大别山超高压变质带为例. 地球物理学报, 44(6): 754–760.
石耀霖, 范桃园. 2002. 大洋板块拖曳窄条陆壳俯冲——大规模超高压变质形成的构造条件. 地质学报, 76(1): 77–82.
石耀霖, 曹建玲. 2010. 库仑应力计算及应用过程中若干问题的讨论——以汶川地震为例. 地球物理学报, 53(1): 102–110. DOI:10.3969/j.issn.0001-5733.2010.01.011
石耀霖, 张贝, 张斯奇, 等. 2013. 地震数值预报. 物理, 42(4): 237–255.
张贝, 张怀, 石耀霖. 2015. 有限元模拟弹性位错的等效体力方法. 地球物理学报, 58(5): 1666–1674. DOI:10.6038/cjg20150518