2. 中国地震局兰州地震研究所, 兰州 730000
2. Lanzhou Institute of Seismology, China Earthquake Administration, Lanzhou 730000, China
由于青藏高原的地形地貌、深部结构和地质构造极具特殊性,其动力学和运动学过程相当复杂,它的演化过程广泛引起地学家的关注.自有限元方法发展到地学领域以来,在青藏地区开展了不少有限元分析和研究,获得了一些较有价值的研究成果(Tapponnier and Molnar, 1976; 傅容珊等, 2000; 陆师阔和蔡永恩, 2004; 邓宗策等, 1990; 戴黎明等, 2011).有限元方法已成为研究和探索该特殊地区的动力学过程的重要方法之一.
在应用有限元法对青藏高原进行的构造应力场模拟方面,地学家们认为印度板块以一定的速度持续向北挤入,这是青藏高原隆升的主要动力.傅容珊等的研究认为挤压过程主要是受岩石的力学特性、边界条件以及剥蚀作用等因素的影响,并认为高原隆升过程是一个非均匀的演化过程(傅容珊,2000;郑洪伟等, 2006).陈开平等的研究认为印度大陆和欧亚大陆碰撞是川滇菱形块体向东南运动的动力源(陈开平和马瑾, 1995).杨立强等(2006)对青藏高原壳幔形变模拟分析认为,青藏高原3D壳幔形变主要特征表现为主碰撞带在纬向上存在距离上的差异,在经向上地壳物质存在“逃逸”现象.张东宁等(1995)用黏弹性幂律蠕变本构关系探究了高原内部张性构造应力分布状态的可能机制.在此基础上,以现代构造活动为主,建立了青藏地区构造形变场有限元模型(张东宁等, 2007).王辉等(2006)利用二维有限元模型,以GPS资料为约束条件,探讨了青藏地区活动块体在应变积累过程中的变形与运动状态.
在应用有限元法对青藏高原进行的地震机制与地震预测模拟方面,曹雪峰等(2003)研究显示南北地震带上的地震存在呼应特征,认为这可能是区域应力场平衡的结果,是地块边界力(外载荷)作用的整体表现.陈祖安等(2008)模拟分析了1997年玛尼7.5级地震对青藏地区地块系统稳定性的影响.获得了该地震引起的地块边界带上库仑应力的变化,其结果表明除了发震断层两端库仑破裂应力增大和更加集中外,2001年昆仑山8.1级地震发震断层段(东昆仑断裂中段)上的库仑破裂应力增加了约2 MPa,2008年改则6.9级地震发震断层段(喀拉昆仑断裂带)上增加了约0.7 MPa,其应力的加载使得这2次地震的发震断层(接近破裂强度)处于失稳状态,加速了地震的发生.陈连旺等(2008)研究认为,在一次强震发生之后,发震断层本身是处于快速卸载的过程,而在其他潜在地震孕育活动断裂带上库仑破裂应力处于加载的过程,可能加载作用起到主要作用,其相互作用主要表现为应力的加载,一次强震的发生可能对下一个强震的孕育过程起到加速的作用,引起强震或多次地震的发生,直到区域应变能达到低值水平状态,区域地震将进入一个新的平静时段.
以上两方面的大尺度研究,其有限元模型是平均化意义上的三维模型,一般把地块当成均匀实体来处理,也就是说同一地块是用一组力学参数来描述的,但实际上地块内部间的非均匀性是明显存在的,使用更为精细的模型进行模拟计算,其结果无疑更加接近实际.目前能够反映整个模型区域三维结构的模型是全球地壳CRUST 1.0模型,该模型主要把地壳分为上、中、下三层,水平分辨尺度为1°×1°.CRUST 1.0全球地壳模型是本研究模型建立的重要数据源.本文以2001年11月14日昆仑山口8.1级、2008年5月12日汶川8.0级和2015年4月25日尼泊尔8.1级3次特大地震为例,应用有限元数值模拟方法,研究震级相近的大地震在不同地区和不同构造条件下产生的物理场特征及其差异性.
2 建立有限元黏弹模型 2.1 模型区域与实体模型的建立模型区域包括青藏地块区的6个地块、塔里木地块南部、阿拉善地块东南部、鄂尔多斯地块西部、华南地块西部、滇缅地块区的滇南地块和滇西地块,共包括12个地块和15条边界带(张培震等, 2003; 张国民等, 2005),其中Ⅰ级地块区边界有8条,次级地块边界有7条.区域东西长约2800 km, 南北宽约1800 km.根据全球地壳模型CRUST 1.0数据,可以获得模型区上、中、下地壳1°×1°网格的P波速度、S波速度、介质密度和几何参数等数据,利用这些参数可计算网格中的杨氏模量和泊松比.上地壳为完全弹性模型.中、下地壳为黏弹模型,其黏滞系数分别设置为1019和1020 Pa·s,各层黏性是均匀的.地壳上、中、下三层模型由1521个实体拼接而成,完成了实体模型的建立.
2.2 模型网格划分ANSYS系统中可提供使用的单元类型有上百种,我们所分析的模型是不规则的,希望或要求模型在施加外力后,其相应的变形量小.为此可以优先选用solid187单元,其单元为具有6个中间节点的高阶三维固体结构的单元,具有二阶位移插值精度,可以更好地表示不规则模型,支持大变形和大应变等模型.为了描述剪切松弛核函数G(t)和体积松弛核函数K(t),我们采用的是Prony级数形式.在网格划分时,采用网格单元平均长度为20 km的四面体单元进行模型非结构化网格自动划分,共划分为488707个单元,节点数为740841个.
2.3 模型断层参数设定断裂带的描述有2种形式,一种是用曲线线条数表示,即一条断裂是用多条线段连接而成的线条;另一种是用一系列的四边形区域组成的一条不等宽度区带来表示,这是一种接近断裂带实际宽度的描述.第一种形式的断裂可根据张国民(2005)的对地块边界带的划分宽度进行给定.对位于这些断裂带或地块边界带内所有单元的材料参数进行弱化(陈连旺等, 2011).材料参数弱化规则为:上地壳弹性层的断层参数:杨氏模量=原杨氏模量/3, 泊松比=原泊松比+0.02;中、下地壳黏弹层的断层参数:杨氏模量=原杨氏模量/3, 泊松比=原泊松比+0.02,黏性系数=原黏性系数/2.
2.4 模型边界条件根据甘卫军(私人交流)提供的GPS数据,对模型区进行数据挑选,对稀疏的区域进行必要的数据插值,在喜马拉雅以南地区GPS测点非常少,仅有几个点,应用这几个测点的GPS数据,在平面上增补了几个控制点.在此基础上,对GPS数据进行二维平面插值,获得要求密度点位的年位移,可对模型侧边界进行位移边界约束,在侧边界竖直方向位移条件是等同的.模型表层为自由面,约定下地壳底界面垂直于地表方向不产生位移,其分向位移为零.在模型构造结处进行水平位移约束,即水平位移分量为零.在模型的南边界向北推挤的最大年位移量可接近50 mm,模型侧边界不同区面位移边界约束大小和方向都存在一定差异.
2.5 地震断层位错加载根据昆仑山口8.1级(震中90.90°E, 36.40°N)、汶川8.0级(震中103.40°E, 30.95°N)和尼泊尔8.1级(震中84.70°E, 28.20°N)3次大地震(图 1a)的实际野外考察和同震资料,昆仑山口8.1级地震断层为左旋走滑断层,倾角陡(80°),极震区平均水平位移滑动量为6 m,断裂长度为350 km(陈文彬等, 2001);汶川8.0级地震断层为右旋逆走滑(徐锡伟等, 2008),倾角较陡(60°),全长260 km,南端三江—映秀段,水平位移为0~3 m,垂直位移0~3.5 m.映秀—平驿段,平均水平位移为3~3.5 m,平均垂直位移3.5~3.7 m.平驿—青川,水平位移为3.5~0 m,垂直位移3.7~0 m.尼泊尔8.1级地震断层为逆断层,倾角为11°,全长约120 km,断层上盘向南的最大位移量为3.5 m(张贝等, 2015).由于这些地震都发生在地块的边界带上,边界带的宽度不小于60 km,满足断裂在深部也位于相应的边界带内的情况.地震位错断层加载宽度60 km(断层两边各30 km),深度为25 km(假定震源深度为15 km),对于走滑位错量,两盘各加载一半,方向相反.对于逆冲或正断位错量,只对上盘进行加载.震源深度至上地壳顶部为均匀加载,其他方向为线性衰减加载(最小加载位错量为0.1 m).
由于模式为黏弹性模型,在计算时采用的是时间步计算,参考时间为1900年或其他时间,参考时间之前以1万年为起算年,从参考时间起按黏弹性模型进行计算,每个地震需要计算两步,地震发生之前为第一步,地震位错加载后为第二步,按地震发生时间依次进行时间步计算,并对目标地震位错加载计算结果与其地震加载之前的计算结果进行差值处理,可以获得目标地震产生的位移场和应力应变场,下面重点对上地壳中的各种物理场进行分析.
3.1 位移分布特征2001年11月14日昆仑山口8.1级地震位于巴颜喀拉地块与柴达木地块的边界东昆仑带的中段,野外考察极震区平均水平位移滑动量为6 m,实际计算上地壳顶界面破裂面极震区的平均位移为5.9 m,水平位移量大于0.4 m的分布区域位于地震破裂带的南北两侧,存在一定的斜对称性,在地块内部存在明显的非均匀性.1 cm以上位移量的区域基本上遍及青藏高原中部和北部地区,以及塔里木地块的东南边缘(图 1b),地块边界为强衰减带.位移方向基本以发震断层为界,北部向西,南部向东,与同震GPS观测一致(万永革等, 2004),但在较远区域位移方向存在一定的变化.2008年5月12日汶川8.0级地震位于巴颜喀拉地块与华南地块的边界岷江—龙门山带上,野外考察极震区平均水平位移为3~3.5 m,平均垂直位移为3.5~3.7 m.实际计算上地壳顶界面破裂面极震区的平均位移量为3.3 m,水平位移量大于0.3 m的分布区域位于地震断层两侧,有一定的对称性.1 cm以上位移量的区域也位于断层带两侧(图 1c),西部区域范围大,巴颜喀拉地块北边界为明显的衰减带.东部区域同等位移量分布范围小.位移方向为西部区域向东、北部区域为东南、东部区域方向为NWW,与GPS同震位移观测结果基本一致(徐韶光等, 2010),南部区域以安宁河—小江带为界位移方向反向.与昆仑山口8.1级地震类似,地块边界为强衰减带,位移大范围分布呈现出非均匀性和不对称性.2015年4月25日尼泊尔8.1级地震位于喜马拉雅带的中段,发震断层为逆冲断层,最大位错量为3.5 m,实际计算接近3.6 m,位移量大于0.3 m的分布区域主要位于地震破裂带以北的NNE区域的拉萨地块中部地域,呈口字型分布,每经过一条边界带,其位移衰减是明显的.位移方向表现为极震区以南为NNE,北部地区向南或指向震中(图 1d).这次地震主要造成青藏高原向南运动.
3.2 应变场分布特征昆仑山口8.1级地震引起的加卸载区各2个,卸载区域位于破裂带的南北两侧,等效应变和最大剪应变为负值,等效应变减少0.01×10-5~0.12×10-5,最大剪应变减少0.01×10-5~0.07×10-5.加载区域位于破裂带两端,等效应变和最大剪应变为正值,等效应变增加0.01×10-5~0.12×10-5,最大剪应变增加0.01×10-5~0.07×10-5 (图 2a—2b).这次地震之后的5年时间里,青藏高原70%的5级以上地震和约90%的6级以上地震发生在加载区或加卸载区交界部位.汶川8.0级地震引起的加卸载区各有2个,卸载区分别位于四川盆地、巴颜喀拉地块东部地区和川滇地块最北部地区,等效应变减少0.01×10-5~0.07×10-5,最大剪应变减少0.01×10-5~0.045×10-5.加载区分别是地震破裂带的两端、川滇地块中部地区、以及北部和西北部地区,即岷江—龙门山带中北段(编号10)、安宁河—小江带北段(编号11),等效应变增加0.01×10-5~0.07×10-5,最大剪应变增加0.01×10-5~0.045×10-5 (图 2c—2d).2010年的玉树7.1级地震、2013年的芦山7.0级地震和岷县漳县6.6级地震以及强余震都发生在加载区.尼泊尔8.1级地震引起的加卸载区各1个,极震区北部地区大范围为卸载区,等效应变减少0.01×10-5~0.12×10-5,最大剪应变减少0.01×10-5~0.07×10-5.加载区主要分布于喜马拉雅带内,等效应变增加0.01×10-5~0.09×10-5,最大剪应变增加0.01×10-5~0.045×10-5(图 2e—2f),7级以上强余震都发生在加载区.
一次大地震的发生,其发震断层两盘会产生相对运动,运动的位移称为地震位错.我们可以通过野外考察的地震位错数据来模拟计算深部地震时的应力场状态,以及地震产生的应变场分布,帮助我们了解需要多大的应力才能产生实际大小的位错.也就是说发生相等大小的地震需要释放多大的应力和应变,用地震时的应力和应变分别与地震前的应力和应变相减,在断层处的值为正值;如果用地震前的值减去地震后的值则结果为负值,这个量值我们可以理解为是地震释放的相关量值.这3次大地震所产生的等效应力和等效应变、最大剪应力和最大剪应变如图 3和图 4所示,3次地震的震源体大小存在区别,尤其在不同深度上,应力和应变分布差异明显.汶川地震的应力应变分布状态在不同深度上差异不大,这一点可以说明震级相当的这3次大地震在水平方向产生位移场、应力和应变场存在较大的变化,汶川8.0级地震所产生的应力应变场分布于整个地壳,而昆仑山口大地震和尼泊尔大地震所产生的应力主要分布于中、上地壳,应变场在下地壳明显减弱.这可能与地震断层性质和所处的空间介质结构有关.一次特大地震发生,释放的应力可达10 MPa以上,产生的最大剪应力为6 MPa以上,产生的等效应变接近3.0×10-4,最大剪应变可大于1.0×10-4(图 3和图 4).
我们应用有限元方法对青藏高原及周缘的三次特大地震进行数值模拟和分析,可得到如下的结论.
(1) 这3次特大地震产生的位移场、应力应变场在空间分布上存在异同,昆仑山地震断层为纯走滑型,其位移场、应力应变场呈4象限分布,与左旋走滑断层模型吻合;尼泊尔地震断层为逆冲型,其位移场、应力应变场呈2象限分布,与断层模型吻合;汶川地震断层为逆走滑型,其位移场、应力应变场呈4象限分布,但分布特征较昆仑山和尼泊尔地震复杂.
(2) 一次大地震的发生,会对周围应力应变场环境引起改变,形成明显的加卸载区(正负值相间的区域),应力应变加载区是强余震和未来几年内的强烈地震主要活动地区.
(3) 一次大地震的发生,在极震区可产生10 MPa以上的应力降,应变的改变量可达2.5×10-4以上.不同地震断层性质和不同深部结构对地震产生的位移场、应力应变场在空间分布上存在明显的差异.
在应用有限元方法对地学模型进行数值模拟时,我们应该看到,计算结果依赖于地球内部结构、边界约束条件和地震位错等数据的合理性和精度.随着地球内部的精细探测以及高密度和高精度地壳运动观测,可以获得精度更高的地壳结构模型和地壳运动的时间位移模型,同时可以应用现代观测手段和测量技术得到更高精度的地震断层位错数据.以此为基础,可对每次大地震进行模拟计算,获得应力应变加卸载分布区域,为未来地震发生研究提供区域依据.
Cao X F, Yin J Y, Yang L M. 2003. Study of simulation for the evolution of stress field on north-south seismic belt in China. Northwestern Seismological Journal , 25 (3) : 221-225. | |
Chen K P, Ma J. 1995. Numerical analysis of tectonic deformation by continental collision between India and Eurasia. Seismology and Geology , 17 (3) : 277-284. | |
Chen L W, Zhang P Z, Lu Y Z, et al. 2008. Numerical simulation of loading/unloading effect on Coulomb failure stress among strong earthquakes in Sichuan-Yunnan area. Chinese J. Geophys. , 51 (5) : 1411-1421. | |
Chen L W, Zhan Z M, Ye J Y, et al. 2011. Numerically modeling the influence of rheological properties on tectonic deformation of Tibet plateau. Journal of Geodesy and Geodynamics , 31 (3) : 8-14. | |
Chen W B, Xu X W, Zhang Z J, et al. 2001. Preparatory field investigation on the surface deformation zone of the Qinghai-Xinjiang MS8.1 earthquake on Nov. 14, 2001. Northwestern Seismological Journal , 23 (4) : 313-317. | |
Chen Z A, Lin B H, Bai W M. 2008. 3-D numerical simulation on influence of 1997 Mani earthquake occurrence to stability of tectonic blocks system in Qingzang and Chuandian zone using DDA+FEM method. Chinese J. Geophys. , 51 (5) : 1422-1430. | |
Dai L M, Li S Z, Tao C H, et al. 2011. Three-dimensional numberical modeling of activity of the Longmenshan fault zone driven by the India Plate. Progress in Geophys. , 26 (1) : 41-51. DOI:10.3969/j.issn.1004-2903.2011.01.004 | |
Deng Z C, Cui J W, Feng X F, et al. 1990. The analysis of tectonic dynamics and deformation mechanics of Qinghai-Xizang plateau by elastic finite element method. Bulletin of the Chinese Academy of Geological Sciences , 21 : 65-67. | |
Fu R S, Xu Y M, Huang J H, et al. 2000. Numerical simulation of the compression uplift of the Qinghai-Xizang plateau. Chinese J. Geophys. , 43 (3) : 346-355. | |
Lu S K, Cai Y E. 2004. Review of numerical simulation on the dynamics of Qinghai-Xizang plateau. Acta Seismologica Sinica , 26 (5) : 547-559. | |
Tapponnier P, Molnar P. 1976. Slip-line field theory and large-scale continental tectonics. Nature , 264 (5584) : 319-324. DOI:10.1038/264319a0 | |
Wan Y G, Wang M, Shen Z K, et al. 2004. Co-seismic slip distribution of the 2001 west of Kunlun mountain pass earthquake inverted by GPS and leveling data. Seismology and Geology , 26 (3) : 393-404. | |
Wang H, Zhang G M, Shi Y L, et al. 2006. Numerical simulation of movement and deformation of Qinghai-Tibet plateau. Journal of Geodesy and Geodynamics , 26 (2) : 15-23. | |
Xu S G, Xiong Y L, Liao H, et al. 2010. Static and kinematic analysis of coseismic deformation of Wenchuan MS8.0 earthquake. Journal of Geodesy and Geodynamics , 30 (3) : 27-30. | |
Xu X W, Wen X Z, Ye J Q, et al. 2008. The MS8.0 Wenchuan earthquake surface ruptures and its seismogenic structure. Seismology and Geology , 30 (3) : 597-629. | |
Yang L Q, Du J, Chen Y. 2006. Numerical modelling of the crust/mantle deformation in the Tibetan plateau. Earth Science Frontiers , 13 (5) : 360-373. | |
Zhang B, Cheng H H, Shi Y L. 2015. Calculation of the co-seismic effect of MS8.1 earthquake, April 25, 2015, Nepal. Chinese J. Geophys. , 58 (5) : 1794-1803. DOI:10.6038/cjg20150529 | |
Zhang D N, Xu Z H. 1995. An alternative interpret ation of earthquake activity caused by normal fault of upper crust in southern Qinghai-Xizang (Tibetan) plateau. Acta Seismologica Sinica , 17 (2) : 188-195. | |
Zhang D N, Yuan S Y, Shen Z K. 2007. Numerical simulation of the recent crust movement and the fault activities in Tibetan Plateau. Chinese J. Geophys. , 50 (1) : 153-162. | |
Zhang G M, Ma H S, Wang H, et al. 2005. Boundaries between active-tectonic blocks and strong earthquakes in the China mainland. Chinese J. Geophys. , 48 (3) : 602-610. DOI:10.1002/cjg2.693 | |
Zhang P Z, Deng Q D, Zhang G M, et al. 2003. Active tectonic blocks and strong earthquakes in the continent of China. Science in China Series D:Earth Sciences , 46 (Supppl. 2) : 13-24. | |
Zheng H W, Li Y D, Gao R, et al. 2006. The advance of numerical simulation in geodynamics. Progress in Geophys. , 21 (2) : 360-369. | |
曹雪峰, 尹京苑, 杨立明. 2003. 南北地震带应力场演变的数值模拟研究. 西北地震学报 , 25 (3) : 221–225. | |
陈开平, 马瑾. 1995. 印度与欧亚大陆碰撞构造变形数值分析. 地震地质 , 17 (3) : 277–284. | |
陈连旺, 张培震, 陆远忠, 等. 2008. 川滇地区强震序列库仑破裂应力加卸载效应的数值模拟. 地球物理学报 , 51 (5) : 1411–1421. | |
陈连旺, 詹自敏, 叶际阳, 等. 2011. 流变特性对青藏高原构造变形影响的数值模拟. 大地测量与地球动力学 , 31 (3) : 8–14. | |
陈文彬, 徐锡伟, 张志坚, 等. 2001. 2001年11月14日青新交界MS8.1地震地表破裂带的初步调查. 西北地震学报 , 23 (4) : 313–317. | |
陈祖安, 林邦慧, 白武明. 2008. 1997年玛尼地震对青藏川滇地区构造块体系统稳定性影响的三维DDA+FEM方法数值模拟. 地球物理学报 , 51 (5) : 1422–1430. | |
戴黎明, 李三忠, 陶春辉, 等. 2011. 印度板块挤压驱动龙门山断裂带活动的三维数值模型. 地球物理学进展 , 26 (1) : 41–51. DOI:10.3969/j.issn.1004-2903.2011.01.004 | |
邓宗策, 崔军文, 冯晓枫, 等. 1990. 青藏高原构造动力及变形机制的弹性有限元分析. 中国地质科学院院报 , 21 : 65–67. | |
傅容珊, 徐耀民, 黄建华, 等. 2000. 青藏高原挤压隆升过程的数值模拟. 地球物理学报 , 43 (3) : 346–355. | |
陆师阔, 蔡永恩. 2004. 青藏高原动力学数值模拟方法与研究进展. 地震学报 , 26 (5) : 547–559. | |
万永革, 王敏, 沈正康, 等. 2004. 利用GPS和水准测量资料反演2001年昆仑山口西8.1级地震的同震活动分布. 地震地质 , 26 (3) : 393–404. | |
王辉, 张国民, 石耀霖, 等. 2006. 青藏活动地块区运动与变形特征的数值模拟. 大地测量与地球动力学 , 26 (2) : 15–23. | |
徐韶光, 熊永良, 廖华, 等. 2010. 汶川地震同震形变的静态和动态分析. 大地测量与地球动力学 , 30 (3) : 27–30. | |
徐锡伟, 闻学泽, 叶建青, 等. 2008. 汶川MS8.0地震地表破裂带及其发震构造. 地震地质 , 30 (3) : 597–629. | |
杨立强, 邓军, 陈赟. 2006. 青藏高原壳幔形变数值模拟研究. 地学前缘 , 13 (5) : 360–373. | |
张贝, 程惠红, 石耀霖. 2015. 2015年4月25日尼泊尔MS8.1大地震的同震效应. 地球物理学报 , 58 (5) : 1794–1803. DOI:10.6038/cjg20150529 | |
张东宁, 徐忠淮. 1995. 青藏高原南部上地壳正断层地震活动的一种可能解释. 地震学报 , 17 (2) : 188–195. | |
张东宁, 袁松涌, 沈正康. 2007. 青藏高原现代地壳运动与活动断裂带关系的模拟实验. 地球物理学报 , 50 (1) : 153–162. | |
张国民, 马宏生, 王辉, 等. 2005. 中国大陆活动地块边界带与强震活动. 地球物理学报 , 48 (3) : 602–610. | |
张培震, 邓启东, 张国民, 等. 2003. 中国大陆的强震活动与活动地块. 中国科学 , 33 (增刊) : 12–20. | |
郑洪伟, 李延栋, 高锐, 等. 2006. 数值模拟在地球动力学中的研究进展. 地球物理学进展 , 21 (2) : 360–369. | |