日本MW9.0级地震同震位移场和应力场的有限元模拟

王连捷, 王红才, 王 薇, 孙东生, 崔军文, 乔子江, 赵卫华, 郭彬彬

1)中国地质科学院地质力学研究所, 北京 100081; 2)中国地质科学院地质研究所, 北京 100037

日本MW9.0级地震同震位移场和应力场的有限元模拟

王连捷1), 王红才1), 王 薇1), 孙东生1), 崔军文2), 乔子江1), 赵卫华1), 郭彬彬1)

1)中国地质科学院地质力学研究所, 北京 100081; 2)中国地质科学院地质研究所, 北京 100037

根据同震位移GPS观测数据, 利用有限元法反演了2011年3月11日本MW9.0级地震的断层滑移模式。在此基础上, 计算了日本MW9.0级地震引起的同震位移场和应力场, 给出了位移和应力的分布, 分析了他们的变化规律并与实测结果进行了对比。计算结果表明: 日本MW9.0级地震的静态断层滑移量最大可达25 m。地震引起断层上盘向东位移, 最大位移在震中附近, 可达24.25 m, 日本东北地区向东位移最大可达6 m。震后地表隆起, 隆起幅度可达5.6 m, 隆起的最高点也在震中附近。日本东北地区东海岸附近有一下沉带, 下沉量可达0.8 m。同震地表位移的计算值与GPS测量结果基本一致。地震引起应力变化, 导致震后应力下降。应力变化是不均匀的, 在震中附近约为9.9 MPa, 在深处可达32 MPa, 在日本东北地区地表应力变化小于4.4 MPa。地震引起的应力变化主要是水平应力, 垂直应力基本不变。

日本地震; 断层滑移; 地应力场; 同震位移; 有限元模拟

2011年3月11日, 日本发生了MW9.0级地震。震中位于日本本州东海岸附近海域, 北纬38.1°, 东经142.6°, 震源深度 20 km。 距最近的海岸城市仙台130 km。震源机制解的P轴方位115°, 倾角–4.5°,破裂长度700 km, 断层走向195°, 倾向北西西, 倾角14°(USGS, 2011a, b)。主震前两天发生过4次大于6级的地震, 震后发生了几十次6级以上余震。地震引发海啸, 席卷了沿岸许多地区。福岛核电站受损, 引起核泄漏, 造成严重核污染, 大量人员伤亡和巨额财产损失。

震后一些作者对地表同震位移进行了研究, 给出了日本岛上同震位移GPS观测资料, 并进行了断层滑移分布的反演和同震位移的计算(陈运泰等, 2011; Pollitz et al., 2011; Sato et al., 2011; Ozawa et al., 2011; Ide et al., 2011; Simons et al., 2011; Hayes, 2011)。一些作者进行了库伦应力的计算(Stein et al., 2011)。但对地震引起的主应力的变化(本文称为同震应力)却研究得很少。而同震应力场对地震机理的研究和震前应力场的重建是有意义的(Yamashita et al., 2004)。因此, 本文在已有研究的基础上, 用弹性力学有限元法进行了断层滑移分布的反演, 同时计算了同震位移场和应力场的分布。为震前应力场的重建和地震机理的研究提供了基础。

1 构造动力背景

2011年3月11日在日本东北地区附近海域发生的Mw9.0级地震, 是由太平洋板块和欧亚板块(也称北美板块)之间的俯冲带(日本海沟)附近的逆冲断层引起的(USGS, 2011a, b)。太平洋板块以大约80 mm/年的速度向西, 朝欧亚板块移动如图1(孙文斌等, 2004; 彭华等, 2011), 引起欧亚板块变形并形成北西西方向的压应力场(Lin et al., 2011)。随着太平洋板块向欧亚板块的推挤, 板块边界(日本海沟)附近应力逐渐增加, 当应力积累达到断层的摩擦强度时,断层突然滑动, 能量突然释放, 形成这次地震。

2 计算模型与计算方法

图1的矩形框为计算模型的范围。图2为有限元计算模型, 模型的Y轴(长边)方向为北北东方向,其长度为1100 km, 模型的X轴(短边)方向为北西西方向, 其长度为542 km, 深度55 km。模型为均匀各向同性弹性介质, 参照文献(殷有泉等, 1982)将弹性模量取为70000 MPa。泊松比为0.22。约束条件为:模型的底部为Z(铅直)方向约束, 左边界X方向(水平方向)约束, 前后边界为Y方向约束。使用矩形网格, 共划分网格5049个, 节点6256个, 断层面上节点328个。

用有限元法进行位移场和应力场分析时, 通常有两种方法: 一是有限元叠加法。这种方法需要先求出断层面上各个节点的单位载荷在GPS测点处引起的位移分量, 根据所有观测点的观测数据, 通过求解联立方程得出断层面上各个节点的载荷, 进而求出断层面的滑移分布(李志明等, 1997; 王连捷, 1991; 王薇, 2004; Wang, 1983)。第二种方法是给定断层面上节点的倾向和走向位移分布, 进行试算并将结果与GPS观测值进行对比。如不满意, 调整位移分布, 重新计算。经过多次试算, 直到得出满意的结果(王仁等, 1982)。

我们选定断层产状为: 走向195°, 倾向北西西75°, 倾角14°(Pollitz et al., 2011)如图2, 并假定只有倾向滑动分量。用第二种方法对断层面上的滑移分布进行了反演。得出了较为满意的结果, 并同时给出了位移场和应力场的分布。

图1 板块运动示意图(彭华等, 2011)Fig. 1 Sketch diagram showing plate motion (after PENG et al., 2011)

3 计算结果

图2 计算模型及网格Fig. 2 Model with mesh for calculation

3.1 断层滑移分布

用上述方法得到的上盘断层面的滑移分布如图3(为作图方便, 我们将计算结果只显示北纬34°到41°之间的区域)。由图3可以看出沿上盘断层面的静态位移可达25 m, 最大静态位移在震源附近。滑动方向沿断层面向上。

3.2 地表同震位移

图4为模拟计算得到的地表同震水平位移等值线云图。可以看出水平位移的分布情况, 水平位移最大的地区在震中附近, 最大水平位移可达24.25 m。远离震中位移逐渐减小, 东海岸一带位移为4～6 m, 向西位移变小, 到西海岸附近位移小于2 m。水平位移方向总体上向南东东方向位移(图6)。

图5为地表同震铅直位移等值云图。可以看出,地表同震铅直位移同样在震中区最大,铅直位移造成地表隆起和下沉, 隆起幅度可达5.6 m。在隆起的边缘地区, 日本东海岸附近有一个下降区, 下降幅度可达0.8 m。

图3 由GPS反演的上盘断层滑移等值线及位移矢量分布在地表的投影Fig. 3 Fault slip model projected onto the surface inverted from GPS data by FEM

图4 地表同震水平位移(计算值)Fig. 4 Co-seismic surface horizontal displacement(calculated value)

将模拟结果与GPS观测结果进行了对比如图6和图7。图6a和图6c为GPS观察值, 图6b和图6d为模拟计算的结果, 图7为海底同震位移观测值与计算值。

由对比可以看出观测值与计算值基本上是一致的。

同震位移的观察值和计算值的一致性, 说明断层面上的滑移分布基本上是符合实际的。

3.3 同震应力场

图8是模拟计算得到的地表同震应力场主应力等值线云图。它不是地应力的绝对值, 是地震引起的应力的变化值, 是增量。由图可见, 同震应力分布是不均匀的, 靠近震中地区应力值最高, 可达9.9 MPa(张应力)。远离震中区应力变小, 日本东海岸附近应力为3.3～4.4 MPa。西海岸附近应力小于1 MPa。从剖面上看(图9), 在深处断层面附近有一个高应力区, 应力可达32 MPa。由路径图(图10)可以看出, 地表应力由西向东(由B到B1)应力的变化情况。西部应力较小, 约0.5 MPa, 且变化平缓。在250 km处, 应力开始较快增加, 在震中附近主应力S1达到9.9 MPa, S2达到2.8 MPa, S3基本不变。

同震主应力S1及S2近水平, S3近于铅直。主应力S1的方位随地区位置而变化如图11。在日本东北地区北部, 主应力S1方位为北西向, 南部为北东向,在海区主应力方位为近东西向。

4 讨论与结论

图5 地表同震铅直位移Fig. 5 Co-seismic surface vertical displacement

用有限元进行应力场模拟时, 模型的范围要足够大, 以使得边界对变形或位移的影响尽量减小。本文计算模型的边界范围较适当, 边界对研究区域内的计算结果影响较小。

图6 同震位移的模拟结果与GPS观测值的比较Fig. 6 Comparison of co-seismic displacement of observation with calculation

由GPS数据用有限元法可反演断层滑移模式,进而计算同震位移场和应力场。反演得到的日本MW9.0级地震的静态断层滑移量最大约为25 m。反演得到的地表同震水平位移在震中附近最大, 可达24.25 m。离开震中位移逐渐减小。日本东北地区东海岸一带水平位移为4～6 m, 西海岸附近小于2 m。水平位移方向总体向东, 向太平洋板块方向运动,但不同地区位移方向略有不同, 有些地区位移方向大体上指向震中。

地震引起上盘总体隆升和局部下沉。隆升幅度同样是震中附近最大, 可达5.6 m。在日本东北东海岸附近有一个下降区带, 下降幅度可达0.8 m。模拟得到的同震位移与GPS观测结果很一致。

图7 海底同震位移GPS观测值(Sato et al., 2011)与计算值的比较Fig. 7 Co-seismic horizontal displacement on the seafloor

地震引起的地应力的变化是不均匀的, 地表靠近震中附近地区主应力S1最大, 其数值约9.9 MPa(张应力, 相对于震前应力下降), 远离震中区应力变小, 日本东海岸附近应力为3.3～4.4 MPa。西海岸附近应力小于1 MPa。深处应力较高, 可达32 MPa。主应力方位随地区位置有所变化, 总体上S1方位为近东西向, 北部和南部分别为北西向和北东向。

图8 同震地表主应力S1的等值云图Fig. 8 Distribution of co-seismic principal stress S1on the surface

图9 沿图8中B-B1剖面同震主应力S1的分布Fig. 9 Distribution of co-seismic principal stress S1in profile of Fig. 8

图10 沿图8的B-B1路径同震应力的变化(S1, S2, S3, 为主应力)Fig. 10 Change of co-seismic stress along Path B-B1in Fig. 8(S1, S2, S3shows principal stress)

三个主应力S1、S2、S3中, S1、S2近水平, S3近于铅直。S1幅度最大, S2较小, S3基本为0。

地应力测量表明, 2008年汶川MW8.0地震后,表层地应力下降5.6 MPa(郭启良等, 2009), 2001年昆仑山MW8.1地震后地应力下降4.6～9.4 MPa(Liao et al., 2003)。本文计算的日本MW9.0 地震引起的应力降, 在表层可达9.9 MPa, 这与上述地区的实测应力降在数量级上是一致的。

将模型范围扩大延伸到中国东部地区, 可以研究日本MW9.0地震对中国东部地区的断层和地震活动性的影响。这将是下一步的工作。

图11 同震主应力S1的方位Fig. 11 Azimuth of principal stress S1

致谢:对日本高知岩心研究所, 日本南海海洋地球科学技术局(JAMSTEC)林为人教授的帮助, 审稿人提出的宝贵意见, 表示感谢！

陈运泰, 张勇, 许力生. 2011. 2011年3月11日日本附近海域Ms8. 6地震[OL/EB]. [2011-08-06]. http://www.cea-igp.ac.cn/.

郭启良, 王成虎, 马洪生, 王崇艮. 2009. 汶川Ms8. 0级大震前后的水压致裂原地应力测量[J]. 地球物理学报, 52(5): 1395-1401.

李志明, 张金珠. 1997. 地应力与油气勘探开发[M]. 北京: 石油工业出版社.

彭华, 马秀敏, 姜景捷, 彭立国. 2011. WFSD地应力台应变特征及其同震效应分析[J]. 地质力学学报, 17(1): 1-14.

孙文斌, 和跃时. 2004. 西太平洋Benioff带的形态及其应力状态[J]. 地球物理学报, 47(3): 433-440.

王连捷. 1991. 叠加法反演山西某铜矿现今地应力场[C]//第二届全国岩石力学数值计算与模型实验学术研讨会论文集. 上海: 同济大学出版社.

王仁, 黄杰藩, 孙旬英, 安欧, 郭世风, 宋惠珍, 潘善德, 黄庆华,沈淑敏, 杨建华. 1982. 华北地震构造应力场的模拟[J]. 中国科学(B辑), (4): 337-344.

王薇, 王连捷, 乔子江, 王红才. 2004. 三维地应力场的有限元模拟及其在隧道设计中的应用[J]. 地球学报, 25(5): 587-591.

殷有泉, 张宏. 1982. 模拟地震的应变软化的数学模型[J]. 地球物理学报, 25(5): 414-423.

References:

CHEN Yun-tai, ZANG Yong, XU Li-sheng. 2011. The 2011. 3. 11 Japan sea areaMs8. 6 Earthquake[OL/EB]. [2012-08-06]. http: //www.cea-igp.ac.cn/(in Chinese).

HAYES G, USGS. 2011. Result of the Mar 11, 2011 Mw 9. 0 Earthquake Offshore Honshu, Japan[OL/EB]. [2012-08-06]. http://earthquake.usgs.gov/earthquakes/eqinthenews/2011/usc0 001xgp/finite_fault.php.

GEONET. 2011. Hashimoto[OL/EB]. [2012-08-06]. http://supersites. Earth observations.org/honshu.php.

GUO Qi-liang, WANG Cheng-hu, MA Hong-sheng, WANG Chong-gen. 2009. In-situ hydro-fracture stress measurement before and after the Wenchuan Ms8. 0 earthquake of China[J]. Chinese Journal of Geophysics, 52(5): 1395-1401(in Chinese with English abstract).

IDE S, BALTAY A, BEROZA G C. 2011. Shallow Dynamic Overshoot and Energetic Deep Rupture in the 2011 Mw 9. 0 Tohoku-Oki Earthquake[J]. Science, 332(6036): 1426-1429.

LI Zhi-ming, ZHANG Jin-zhu. 1997. In-situ stress and petroleum exploration & development[M]. Beijing: Petroleum Industry Press(in Chinese).

LIAO Chun-ting, ZHANG Chun-shan, WU Man-lu, MA Yin-sheng, OU Ming-yi. 2003. Stress change near the Kunlun fault before and after the Ms 8. 1 Kunlun earthquake[J]. Geophysical Research Letters, 30(20): 2027, doi: 10.1029/2003GL018106.

LIN Wei-ren, SAITO S, SANADA Y. 2011. Principal horizontal stress orientations prior to the 2011 Mw 9. 0 Tohoku-Oki, Japan, earthquake in its source area[J]. Geophysical Research Letters, 38, L00G10, doi: 10. 1029/2011GL049097.

OZAWA S, NISHIMURA T, SUITO H. KOBAYASHI T, TOBITA M, IMAKIIRE T. 2011. Coseismic and postseismic slip of the 2011 magnitude-9 Tohoku-Oki earthquake[J]. Nature, 475: 373-378.

PENG Hua, MA Xiu-min, JIANG Jing-jie, PENG Li-guo. 2011. Strain characteristics of WFSD stress stations and ITS CO-Seismic Effects analysis: A case study upon MW9. 0 E arthquake in Japan[J]. Journal of Geomechanics, 17(1): 1-14(in Chinese with English abstract).

POLLITZ F F, BÜRGMANN R, BANERJEE P. 2011. Geodetic slip model of the 2011 M9. 0 Tohoku earthquake[J]. Geophysical Research Letters, 38: L00G08, doi: 10.1029/2011GL048632.

STEIN R, SEVILGEN, TODA S. 2011. Coulomb stress imparted by the Mw9. 0 earthkuake[OL/EB]. [2011-08-06]. http://supersites.earthobservations.org/honshu.php.

SATO M, ISHIKAWA T, UJIHARA N, YOSHIDA S, FUJITAM, MOCHIZUKI M, ASADA A. 2011. Displacement above the hypocenter of the 2011 Tohoku-Oki earthquake[J]. Science, 332(6036): 1395.

SIMONS M, MINSON S S, SLADEN A, ORTEGA F, JIANG J, OWEN S E, MENG L, AMPUERO A P, WEI S J, CHU R S, HELMBERGER D V, KANAMORI H, HETLAND E, MOORE A W, WEBB F H. 2011. The 2011 Magnitude 9. 0 Tohoku-Oki Earthquake: Mosaicking the Megathrust from Seconds to Centuries[J]. Science, 332(6036): 1424-1425.

SUN Wen-bin, HE Yao-shi. 2004. Characteristics of the subduction zone in the western Pacific and its stress state[J]. Chinese Journal of Geophysics, 47(3): 433-440(in Chinese with English abstract).

The ARIA team at JPL and Caltech. 2011. http://supersites. earthobservations.org/honshu.php.

USGS. 2011a. Magnitude 9.0- Near the East coast of Honshu, Japan[OL/EB]. [2011-08-06]. http://earthquake.usgs.gov/ Earthquakes/eqinthenews/2011/ usc0001xgp.

USGS. 2011b. Magnitude 9. 0 near the east coast of Honshu, Japan 2011 March 11[OL/EB]. [2011-08-06]. http://earthquake.usgs. gov/earthquakes/eqinthenews/2011/usc0001xgp/#summary.

WANG Lian-jie. 1991. Inversion of present geo-stress field by finite element method in a mine in Shanxi Province // Proceedings of Symposium of Second National Rock Mechanics numerical calculation and model experiments[C]. Shanghai: Shanghai Tongji University Press(in Chinese).

WANG Ren, HUANG Jie-fan, SUN Xun-ying, AN Ou, GUO Shi-feng, SONG Hui-zhen, PAN Shan-de, HUANG Qing-hua, SHEN Shu-min, YANG Jian-hua. 1982. Simulation of the earthquake tectonic stress field of North China[J]. Science in China (Series B), (4): 337-344(in Chinese with English abstract).

WANG Ren. 1983. A short note on the inversion of tectonic stress fields[J]. Tectonophysics, 100(1-3): 405-411.

WANG Wei, WANG Lian-jie, QIAO Zi-jiang, WANG Hong-Cai. 2004. Modeling of Three Dimensional Crustal Stress Field and Its Application in Tunnel Design[J]. Acta Geoscientica Sinica, 25(5): 587-591(in Chinese with English abstract).

YAMASHITA F, FUKUYAMA E, OUKUYAMA K. 2004. Estimation of Fault Strength:Reconstruction of Stress before the 1995 Kobe Earthquake[J]. Science, 306(5694): 261-263.

YIN You-quan, ZHANG Hong. 1982. A mathematical model of strain softening in simulating earthquake[J]. Chinese Journal of Geophysics, 25(5): 414-423(in Chinese with English abstract).

Modeling of the Co-seismic Displacement Field and Stress Field of Mw9.0 Earthquake in Japan by Finite Element Method

WANG Lian-jie1), WANG Hong-cai1), WANG Wei1), SUN Dong-sheng1), CUI Jun-wen2), QIAO Zi-jiang1), ZHAO Wei-Hua1), GUO Bin-bin1)
1) Institute of Geomechanics, Chinese Academy of Geological Sciences, Beijing 100081; 2) Institute of Geology, Chinese Academy of Geological Sciences, Beijing 100037

According to the GPS data, the fault-slip model for MW9.0 earthquake in March 2011 was inversed by using the finite element method. On the basis of the inversion, co-seismic displacement and stress fields were calculated, and the distribution of displacement and stress was given. The results show that fault-slip was up to 25 m. Northeast Japan was moved eastward by 1～6 meters, and the maximum displacement in the epicentral area was up to 24.25 m. After the earthquake, the surface was uplifted by about 5.6 meters near the epicenter. There was a depression area of 0.8 meters in the east coast of Northeast Japan. The calculated co-seismic surface displacement is consistent with GPS measurements. Stress was changed by the earthquake, resulting in the decrease of the post-earthquake stress. Stress change was about 9.9 MPa near the epicenter, 32 MPa in the depths, and less than 4.4 MPa in Northeast Japan. Earthquake-induced stress changes were mainly horizontal stress, whereas vertical stress changes were very small.

earthquake in Japan; fault slid; stress field; co-seismic displacement; finite element modeling

P315.725; P315.8

A

10.3975/cagsb.2013.06.03

本文由国土资源部“深部探测技术与试验研究”专项(编号: S inoProbe-07)和国家科技支撑“汶川地震断裂带科学钻探井中探测”项目(编号: 2008wfsd-03)联合资助。

2012-12-10; 改回日期: 2013-03-13。责任编辑: 闫立娟。

王连捷, 男, 1933年生。研究员。主要从事地质灾害、地应力测量、应力场数值分析及应用研究。E-mail: wanglj01@sina.com。