2. 中国地质科学院地质研究所, 北京 100037
2. Institute of Geology, Chinese Academy of Geological Sciences, Beijing 100037,China
2010年4月14日7时49分40.7秒在我国青海省玉树县(33.23°N,96.61°E)发生MS7.1 级地震,震源深度18km, 该地震发生在甘孜-玉树断裂带上,震源机制得出的断层面走向119°/倾角83°/滑动角-2°[1].这次地震破裂过程持续了约23s, 包括两次主要的子事件,第一次是在初始破裂后0~5.5s;第二次在5.5~23s左右[1].地表破裂主要集中在距震中东南约10~30km 处,破裂长度延续23km, 破裂错动以反扭走滑为主,错动量可达1.75 m[2].地震的破裂主要由北西向东南方向扩展,破裂滑移量最大的区域位于玉树城区附近.
玉树地震位于巴颜喀拉地块南边界,汶川地震位于巴颜喀拉地块东南边界.它们处于同一地块不同的边界上,具有相同的构造动力学背景.
地震造成严重破坏,人员伤亡和财产损失.震后第一时间很多专家赶赴现场进行考查,为重建,减灾救灾和地震科学研究提供了重要依据.在此基础上,我们对玉树地震发震应力场进行了分析,对非稳定发震过程进行了数值模拟,为地震机理的研究和预测预报提供依据.
已有很多文献做了大量地震模拟工作,为地震机理研究和地震危险区预测做出了重要贡献[3, 4].但是以往文献大多没有考虑应变软化的地震非稳定过程,因而不能给出地震产生的断层错动,应力降以及能量释放等重要参数.在以往工作的基础上,我们采用应变软化非稳定模型进行模拟.将围岩看成弹性介质,断层看成具有软化特性的弹塑性介质,采用非稳定断层滑动准则[5],模拟地震的非稳定过程.所用的软件是有限差分FLAC3D 软件[6].
2 玉树地区的应力场和速度场地震是在地应力场的作用下,应力逐渐积累,当应力达到断层的摩擦破坏强度时,断层发生软化,产生突然滑动,储存在围岩中的能量突然释放,形成地震[7~9].因此需要首先查明应力场.
根据地应力测量和震源机制解得到玉树地震区及附近地区地应力场最大水平主应力方向如图 1(图中应力方向据文献[10, 11]绘制,GPS速度场来自文献[12]).在玉树地震区的东南部沿龙门山断裂带最大水平主应力方向为北西到北西西向.在西部那曲附近为北北东到北东向.在北部张掖附近为北东到北东东向.在玉树附近东南方向为北东到北东东向.位于巴颜喀拉地块南边界的喀喇昆仑山口、玛尼、玉树三个地震,他们的破裂类型都表现为左旋走滑运动[13],这说明沿着这个边界主压应力方向为北东到近东西向.玉树MS7.1 级地震的震源机制解,其P轴方位为70°~74°[14].
根据P轴方位可以求出主应力方位.关于P轴方位与主应力方位之间的关系,文献[15]认为他们通常是不同的,他们之间的差别与岩石或断层的内摩擦等因素有关.破裂面与主应力方位之间的夹角θ可用下式表示[8]:
(1) |
式中φ 为内摩擦角.当存在内摩擦时,θ 角小于45°,当内摩擦角φ 为0时θ=45°,此时为P轴与破裂面的夹角,也就是P轴方位是当内摩擦角φ 为0时的主应力方位.考虑到内摩擦的存在,我们将由P轴得到的最大水平主应力的方位定为北东80°,与P轴方位相差6°~10°,这是比较合理的.这个方位与玉树附近的速度方位大体一致.我们将这个方位定为玉树地震发震应力场的方位.因此,在本模型中,取最大水平主应力为北东80°,与断层夹角为40°.
由浅部550 m 深度以内地应力测量结果[11],168个数据的统计得到的水平地应力的大小及其随深度的变化可表示为式(2)及(3):
(2) |
(3) |
可以看出水平地应力是随深度增加的.铅直应力σv可由上覆岩层重量计算,其表达式为
(4) |
式中H为深度(m).σH 为最大水平主应力,σh 为最小水平主应力.应力单位为MPa.三个主应力的关系为:σH >σv >σh, 符合该区走滑断层应力状态.
GPS 测量得到的速度场如图 1.可以看出,玉树附近的速度方向为北东70°到近90°.在玉树的西部速度快,东部速度慢,由西向东速度递减,有一个速度梯度.速度梯度可以定义为dv/dx,他表示在距离为dx的两个质点的速度变化.速度梯度等于应变速度[16].由于该区存在速度梯度,因而产生应变和应力.由玉树附近的西部地区和东部地区的速度差异可以估算出该区的速度梯度约为(1.7~4.2)×10-8 m/a.如果没有速度梯度,速度处处相同,速度梯度和应变速度为零,则只有物体的整体移动,而无变形.因而无应变及应力.
玉树地区的应力场和位移场是青藏高原应力场的一部分,是由于印度板块长期向北推挤,青藏高原向东南方向侧向挤压以及甘孜-玉树断层两侧速度差异的背景下形成的.该地区的应力场和速度场的特点,有利于甘孜-玉树断层产生左旋走向滑动,形成走滑地震.
3 计算模型,本构方程和参数 3.1 计算模型、边界条件和模型网格模型的位置取在玉树附近.断层产状为近直立的倾斜断层,走向120°,倾向210°,倾角80°.模型长度L为80km, 长轴方向北东80°,模型宽度D为15km, 深度H为27km, 如图 2.模型中的断层长度23km.
模型受到重力、地应力和孔隙压力的作用.最大水平主应力σH 的方向与模型长轴平行,与断层走向夹角为40°.最小水平主应力的方向与模型长轴垂直.最大水平主应力σH 的方向主要依据玉树地震的震源机制解并参考GPS和地应力测量结果选定,并认为应力方向不随深度变化.水平地应力大小据浅部地应力测量结果选定.最大水平主应力和最小水平主应力皆随深度变化,其表达式如式(2)及(3).铅直应力由上覆岩层重量确定.考虑孔隙压力p为静水压力,其表达式为
(5) |
式中孔隙压力单位为MPa.
铅直应力、两个水平主应力和孔隙压力均作为初始应力和初始孔隙压力施加在模型上.初始应力和初始孔隙压力由FLAC3D 软件自动施加.图 3给出了初始最大水平主应力的分布,可以看出模型内部初始最大水平主应力随深度均匀增加,在水平方向均匀不变.在地表最大水平主应力为5.5MPa, 在27km 深处为869.5 MPa.分布规律与式(2)相同.最小水平主应力、铅直应力以及孔隙压力的分布也是如此,他们分别与式(3),式(4),式(5)相同.
模型的边界条件只是在模型的左边界施加速度,右边界水平向约束.模型的前后边界为Y方向约束,底边界Z方向约束.
图 4是模型网格.对于围岩使用六面体网格,共划分1500个单元.断层使用不连续界面单元,使用三角形网格.
取库伦条件为断层的屈服条件[5]:
(6) |
式中c和μ 分别是内聚力和内摩擦系数,他们是k的函数;k为塑性体应变或塑性应变;τzx,τzy,σz为断层面的剪应力和法向应力.
当上式满足时,断层滑动.当f(σ,k)<0时,断层处于弹性状态.
3.2.2 断层的本构方程断层单元或节理单元首先由Goodman 等提出[17],Zienkiewicz等进行了发展[18].文献[5, 19, 20]在库伦屈服条件下,用增量理论进行了详细推导,建立了断层单元的本构方程,其表达式如下:
(7) |
(8) |
(9) |
式中kt,kn为断层的切向和法向刚度;τzx,τzy,τ 为沿断层面的剪应力,dτzx,dτzy,dσz为断层面的法向应力的增量;τ=
围岩为弹性体,他的本构方程为:
(10) |
式中De 为弹性矩阵,ε为应变.
3.3 断层非稳定滑动的能量准则断层和围岩组成的地质系统,在断层失稳滑动前的时刻,断层处在非稳定平衡状态,一个小的扰动,如远场位移,应力,孔隙压力等的扰动,系统的弹性应变能发生改变,释放的能量,除了驱动断层滑动外,还能产生地震波的动能.下式表示断层非稳定滑动的能量准则:
(11) |
当上式满足时,断层发生非稳定滑动.式中ΔU是系统弹性应变能的变化,ΔE是系统的耗散能.即断层非稳定滑动时,系统弹性应变能的改变大于断层破坏滑动过程所需的能量,其中的差值用于产生地震波的能量[5].
3.4 介质参数的选择我们参考文献[5, 21, 22]模拟地震时的摩擦系数和内聚力以及他们的软化情况进行了参数选择如表 1.
在初始水平地应力σH、σh、铅直应力σv、孔隙应力p和边界速度v的作用下,随着远场位移(边界位移)的增加,应力逐渐加大.当最大水平主应力的平均值σH 达到435.75 MPa时,断层突然错动(图 5),能量突然释放(图 6),应力σH 突然下降(图 7),此时发生地震.地震造成的应力降为1.0MPa, 释放的能量密度为8×104erg/cm3.能量释放和应力降的大小与断层位错量有关,断层错动量越大,能量释放和应力降也越大.
应该指出,这里所指的应力降与地震学中的应力降不同.地震学中的应力降定义为沿断层面的平均剪应力降.这里是最大水平主应力平均值的应力降.这样做是因为计算起来简单一些.
平均应力即最大水平主应力的平均值是指模型深度范围内的平均值,由下式得出:
(12) |
式中:H0 为最大深度,H0=27km.
按(12)式计算结果,平均应力相当于13.5km深处的最大水平主应力.
4.2 断层的错动图 8为地震时断层面的剪切错动,是水平错动和铅直错动的合矢量.最大错动量约为2m.断层两端错动量最小,靠近中心部位错动最大.在一定的深度范围内,错动量随深度变化不大,超过一定深度后,错动量随深度变小,最深处错动近为零.这种位移分布规律与文献[1]给出的相近.
图 9为地震时地表的同震位移.可以看出,断层为逆时针走向滑动.断层两端位移小,中间部分大.断层两盘的相对滑动量最大约为1.75m.与现场调查结果基本一致.地表破裂长度延续约23km.位移涉及的范围距断层约10km.
图 10为剖面上断层走向位移分布.剖面方向为北东80°,剖面位置为切割模型中央.由图 10可以看出,走向位移随深度变化不大,随离开断层距离的加大,走向位移变小,断层面附近位移最大.
除了走滑位移外,还有少量垂向位移如图 11.由图 11可以看出,地表垂向位移最大约为0.11m.在断层两端出现对称的隆起与下沉.如断层的北西端右侧出现鼓包,左侧出现下沉,断层的东南端左侧出现鼓包,右侧出现下沉.这反映了走滑断层的端部效应.在野外调查中也观察到这种相似的隆起现象[2].
图 12为剖面上沿断层倾向的位移,剖面方向与长边平行.同样可以看出,断层两端对称的隆起和沉降,但其幅度很小,约为9~10cm.
在当前的应力状态和地层参数的情况下,应力增加1.0MPa后,地应力又达到临界状态,即435.75 MPa, 可再发生7.1级地震.地震的复发周期与地应力的积累速度有关.地应力的积累速度越快,复发周期越短.地应力的积累速度与应变速率有关.可由图 7的应力位移曲线求出远场每单位位移引起的应力变化,而单位位移所需要的时间可由位移速度求出,速度可由速度梯度求出,进而求出地应力的积累速度.这里所指的速度是模型东西两侧的相对速度,即两侧的速度差.当速度梯度为(1.7~4.2)×10-8m/a时,求得的应力积累速度为(1.36~3.36)×10-3MPa/a.应变速率与速度梯度相同[11],故应变速率为(1.7~4.2)×10-8/a.这里可以看出,应力积累速度和应变积累速度是非常缓慢的.应力和应变的积累速度对地震预报的研究有重要意义.
在这样的应力积累速度下,应力增加1.0 MPa的时间约为300~700a.即玉树7.1 级地震复发周期约为300~700a, 如图 13 所示.文献[23]给出的玉树7.1级地震的复发周期为200a, 与本文给出的复发周期下限接近.
图 13为地震复发周期与应力积累的关系.由图 13可以看出,在应力积累期间,断层没有明显的位移错动,随着应力的积累,当应力达到临界状态,也就是静摩擦强度时,断层突然滑动,应力突然下降到动摩擦强度[24].然后应力重新积累,直到下一次断层滑动.
5 讨论与结论 5.1 讨论在非稳定地震过程的模拟计算中,要给定一个应力场,包括水平应力和铅直应力,作为初始应力场.应力场的方向根据浅层地应力测量和深部震源机制解等方法给出,并认为应力场的方向不随深度变化.但对深部地应力大小及其分布情况目前我们知道的甚少,建立深部地应力大小的分布比较困难.目前尚无比较好的方法.一种做法是用浅部地应力测量的结果外推到深部,如文献[25]在研究1955年神户地震应力场和断层强度时,曾利用浅部1.6km深度以内的地应力测量结果外推到14.3km 的深度.我们参考这种方法,将浅部地应力测量得到的分布规律外推到深部,作为模型的初始应力.外推后的应力场与本区断层活动类型,即与甘孜-玉树左旋走滑断层类型一致.从模拟结果来看,大体上也是可以被接受的.因此这种外推基本上是可行的.如果没有地应力测量资料作为参考,也要给定一个应力大小的分布.所以如何更合理地建立深部应力状态是需要进一步研究的一个重要问题.
地震的复发周期是由地应力的积累速度估算的.这里采用的应力积累速度由GPS观测的应变速率求出,反映短期的应力积累速度.如果能用地质时期长期应力积累速度进行估算可能更符合实际.
用本文估算地震复发周期的方法,曾估算了汶川MS8.0级地震的复发周期,其结果与其他方法所得结果相近[22].对于5 级或6 级地震,该方法也应该适用,但尚无实例验证.
本文的地震复发模型是一个简单的理想化规则模型.它只与三个参数有关,即应力(或应变)的积累速度、断层的静态摩擦强度和动态摩擦强度有关,没有考虑断层蠕变等因素.地震复发周期是在假定这三个参数不变的情况下做出的,它只适合于估算同样震级的复发周期.实际上,这些参数有可能随着时间改变.如果只有一个参数改变,如应力积累速度改变,则复发周期将延长或缩短.如果这三个参数都发生改变,地震复发周期和震级都将改变,这时,需用不规则模型.对于不同震级交替发生的周期问题,也需用不规则模型,但建立这种不规则模型是比较困难的.
如果能用多种方法,如历史地震法、地震矩率法以及地质学的方法[23, 26],对地震复发周期进行估算,并进行对比分析,可能会得到更符合实际的结果.
模型所施加的应力虽然随深度增加,但应力比值,如最大与最小水平主应力比值不随深度变化,各点基本一样,且断层强度各点也一样.因此断层面上各点是同时或基本同时达到屈服条件的.
5.2 结论(1) 玉树地震是印度板块持续向北挤压,青藏高原物质向东南方向流动以及甘孜-玉树断层两侧速度差异的背景下,在玉树地区形成近东西方向的应力场和速度场,在这样的应力场和速度场作用下,应力逐渐积累,再加上孔隙流体压力的作用,使甘孜-玉树断裂的应力达到其破裂强度,断层突然滑动,发生地震.
(2) 随着水平地应力的缓慢增加,当达到435.75 MPa时(平均值),断层失稳滑动,形成地震.能量释放密度为8×104erg/cm3,地震的应力降为1.0 MPa.断层以反扭走滑为主,挤压为辅.相对走向滑动量为1.75m.
(3) 在断层两端出现对称的隆起与下沉,断层前进端一侧出现隆起,另一端出现下降,反映了走滑断层的端部效应.隆起与下沉的幅度约为0.11m.
(4) 玉树地震区应力的积累速度为(1.36~3.36)×10-3 MPa/a.应变积累速度为(1.7~4.2)×10-8/a.应力积累速度和应变积累速度是缓慢的.应力和应变的积累速度对地震预报的研究有重要意义.
(5) 当玉树地震区的应力积累速度为(1.36~3.36)×10-3 MPa/a时.地震复发周期约为300~700a.当速度改变时,地震复发周期将会改变.
致谢感谢审稿人提出的宝贵意见!
[1] | 张 勇, 许力生, 陈运泰. 2010年4月12日青海玉树Ms7.1级地震. 2010, http://www.cea-igp.ac.cn Zhang Y, Xu L S, Chen Y T. Yushu Ms7.1. Earthquake in Qinghai in 2010. (in Chinese), 2010, http://www.cea-igp.ac.cn |
[2] | 马寅生, 张永双, 胡道功. 地质力学研究所专家组青海玉树地震断裂调查. 2010 http: //www. geomech. ac. cn/xinwen/2010/yushudizhen/yushu0423. htm Ma Y S, Zhang Y S, Hu D G. Qinghai Yushu Earthquake Fracture Investigation of Expert Group from Institute of Geomechanics. (in Chinese), 2010 http: //www. geomech. ac. cn/xinwen/2010/yushudizhen/yushu0423. htm |
[3] | 陈祖安, 林邦慧, 白武明, 等. 2008年汶川8.0级地震孕震机理研究. 地球物理学报 , 2009, 52(2): 408–417. Chen Z A, Lin B H, Bai W M, et al. The mechanism of generation of May 12, 2008 Ms8.0 Wenchuan earthquake. Chinese J. Geophys. (in Chinese) , 2009, 52(2): 408-417. |
[4] | 朱守彪, 张培震. 2008年汶川Ms8.0地震发生过程的动力学机制研究. 地球物理学报 , 2009, 52(2): 418–427. Zhu S B, Zhang P Z. A study on the dynamical mechanisms of the Wenchuan Ms8.0 earthquake. Chinses J. Geophys. (in Chinese) (in Chinese) , 2009, 52(2): 418-427. |
[5] | 殷有泉, 张宏. 模拟地震的应变软化的数学模型. 地球物理学报 , 1982, 25(5): 414–423. Yin Y Q, Zhang H. A mathematical model of strain softening in simulating earthquake. Chinese J.Geophys. (Acta Geophysica Sinica) (in Chinese) (in Chinese) , 1982, 25(5): 414-423. |
[6] | 陈育民, 徐鼎平. FLAC/FLAC3.D基础与工程实例. 北京: 中国水利出版社, 2009 . Chen Y M, Xu D P. FLAC/FLAC3D Case of Foundation and Engineering (in Chinese). Beijing: China Waterpower Press, 2009 . |
[7] | 陈颙. 地壳岩石的力学性能——理论基础与试验方法. 北京: 地震出版社, 1988 . Chen R. Mechanical Properties of Crest Rock—Theoretical Foundation and Experimental Methods (in Chinese). Beijing: Seismological Press, 1988 . |
[8] | 耶格J C, 库克NG W. 岩石力学基础. 中国科学院工程力学研究所译, 北京: 地震出版, 1981 . Jaeger J C, Cook N G W. Fundamentals of Rock Mechanics(in Chinese). Translated by Institute of Engineering Mechanics, Chinese Academy of Sciences (in Chinese). Beijing: Seismological Press, 1981 . |
[9] | Rice J R. Theory of precursory processes in the inception of earthquake rupture. Gerlands Beitr. Geophys. , 1979, 88: 91-127. |
[10] | 谭成轩, 马秀敏,杨 农等. 地质力学研究所地应力监测台站获得青海玉树Ms7.1级地震信息. http: //www. geomech. ac. cn/xinwen/2010/yushudizhen/yushujianbao. htm. 2010[2011-01-25] Tan C X, Ma X M, Yang N, et al. Qinghai Yushu Ms7.1 Earthquake Information Obtained by Crurst Stress Station in Institute of Geomechsnics (in Chinese). http: //www. geomech. ac. cn/xinwen/2010/yushudizhen/yushujianbao. htm. 2010[2011-01-25] |
[11] | 彭华, 马秀敏, 李金锁, 等. 南水北调西线一期工程地壳稳定性研究. 北京: 地震出版社, 2009 . Peng H, Ma X M, Li J S, et al. Study on Crustal Stability of West Route Engineering of South-to-North Water Transeer Project (in Chinese). Beijing: Seismological Press, 2009 . |
[12] | 赖锡安, 黄立文. 中国大陆现今地壳运动. 北京: 地震出版社, 2004 . Lai X A, Huang L W. Current Crustal Movement in Chinese Mainland (in Chinese). Beijing: Seismological Press, 2004 . |
[13] | 刁桂苓, 王晓山, 高国英, 等. 以震源机制类型划分汶川、玉树地震构造块体归属. 地球物理学报 , 2010, 58(3): 1778–1783. Diao G L, Wang X S, Gao G Y, et al. Tectonic block attribution of Wenchuan and Yushu earthquakes distinguished by focal mechanism type. Chinses J. Geophys. (in Chinese) (in Chinese) , 2010, 58(3): 1778-1783. |
[14] | Bureou of Geological Surver of America. Magnitude 6.9-Southern Qinghai, China. http://earthquake.usgs.gov. 2010 |
[15] | 王连捷, 王薇, 袁佳音. 中国地应力测量的进展. 1989 : 21 -32. Wang L J, Wang W, Yuan J Y. Progress of stress measurement in China. In: Selected Pupers on Geomechanics No.9 (in Chinese). Beijing: Geological Publishing House, 1989 : 21 -32. |
[16] | MeansW D. Means W D. 应力与应变. 北京: 煤炭工业出版社, 1980 . Means W D. Stress and Strain (in Chinese). .Beijing: Coal Mining Industry Press, 1980 . |
[17] | Goodman R E, Taylor R L, Brekke T L. A model for the mechanics of jointed rock. Soil Mech. Found. Div., ASCE , 1968, 94(M3): 637-659. |
[18] | Zienkiewicz O C, Taylor R L. The Finite Element Method (Fifth Edition), Volume 1: The Basis, Volume 2: Solid Mechanics, Volume 3: Fluid Dynamics, Xxford, Boston: Butterworth-Heineman, 2000 |
[19] | 殷有泉, 张宏. 断裂带内介质的软化特性和地震的非稳定模型. 地震学报 , 1984, 6(2): 135–145. Yin Y Q, Zhang H. The softening behaviour of fault zone medium and an instability model of earthquakes. Acta Seismologica Sinica (in Chinese) (in Chinese) , 1984, 6(2): 135-145. |
[20] | 殷有泉. 非线性有限元基础. 北京: 北京大学出版社, 2007 . Yin Y Q. Non-Linear Finite Element Method (in Chinese). Beijing: Peking University Press, 2007 . |
[21] | Stuart W D. Quasi-static earthquake mechanics. Reviews of Geophysics and Space Physics , 1979, 17(6): 1115-1120. DOI:10.1029/RG017i006p01115 |
[22] | 王连捷, 崔军文, 周春景, 等. 汶川5·12地震发震机理的数值模拟. 地质力学学报 , 2009, 15(12): 105–113. Wang L J, Cui J W, Zhou C J, et al. Numerical modeling for Wenchuan earthquake mechanism. Journal of Geomechanics (in Chinese) (in Chinese) , 2009, 15(12): 105-113. |
[23] | 任俊杰, 谢富仁, 刘冬英, 等. 2010年玉树地震的构造环境、历史地震活动及其复发周期估计. 震灾防御技术 , 2010, 5(2): 228–233. Ren J J, Xie F R, Liu D Y, et al. Study of tectonics, seismicity and recurrence interval of Yushu 2010 earthquake, Qinghai Province. Technology for Earthquake Disaster Prevention (in Chinese) (in Chinese) , 2010, 5(2): 228-233. |
[24] | 笠原庆一. 地震力学. 北京: 地震出版社, 1984 . Kasahara K. Earthquake Mechanics (in Chinese). Beijing: Seismological Press, 1984 . |
[25] | Yamashita F, Fukuyama E, Omura K. Estimation of fault strength: Reconstruction of stress before the 1995 Kobe Earthquake. Science , 2004, 306(5694): 261-263. DOI:10.1126/science.1101771 |
[26] | 李海兵, 王宗秀, 付小方, 等. 2008年5月12日汶川地震(Ms8.0)地表破裂带的分布特征. 中国地质 , 2008, 35(5): 803–813. Li H B, Wang Z X, Fu X F, et al. The surface rupture zone distribution of the Wenchuan earthquake (Ms8.0) happened on May 12th, 2008. Geology in China (in Chinese) (in Chinese) , 2008, 35(5): 803-813. |