汶川8.0级地震对鄂尔多斯块体西南缘断裂体系影响的数值模拟

刘双庆,盛菊琴,曾宪伟,解朝娣

(1.天津市地震局,天津 300201;2.宁夏地震局,宁夏 银川 750001;3.云南大学,昆明 650091)

汶川8.0级地震对鄂尔多斯块体西南缘断裂体系影响的数值模拟

刘双庆1,盛菊琴2,曾宪伟2,解朝娣3

(1.天津市地震局,天津 300201;2.宁夏地震局,宁夏 银川 750001;3.云南大学,昆明 650091)

利用有限元方法及静态位移解理论,采用前人反演的汶川8.0级地震破裂过程时空参数,定量计算分析汶川地震对鄂尔多斯块体西南缘断裂体系的影响。数值计算表明,具逆冲右旋性质的汶川地震,在汶川—宝鸡一线以东产生顺时针的位移扰动场,以西产生逆时针的位移扰动场,这与该分界两侧附近的活断层整体走滑方向相一致。水平剪切应力场的局部极值区与震后的鄂尔多斯块体西南边界附近的断裂体系小震活动有对应性。震后数值模拟还表明汶川地震使得海原断裂带、香山—天景山断裂带的地震活动性明显增强。剪应力增强和地震暂缺的矛盾解释区也许是未来前兆观测更要密切关注的地区。

汶川地震;有限元模拟;静态位移解;鄂尔多斯块体

0 引言

近年来,关于大地震是如何影响后续地震发生时间和地点的力学机制以及其在地震危险性分析中的应用的定量研究问题,引起了国内外地震学界广泛的兴趣[1-6]。这些研究主要基于线弹性格林函数积分法及库仑应力变化以分析大震对远近断层的应力投影及其与后续地震的对应关系,其主要结论是处于正库仑应力区的断层,后续发震的几率将增加,具有增震的作用;反之则具有减震作用。2008年5月12日,在四川龙门山构造带发生的 Ms8.0级特大地震,造成了巨大的地面破坏和人员伤亡,产生的最大垂直位移5m,最大右旋走滑位移4.8m[7]。而汶川大地震后1年内,在宁夏固原附近发生了2次4级以上地震,香山—天景山等地震带小震活动明显增强,震中分布的带状特征较前7年显著(图7A,C处)。龙门山构造带以北的陕甘宁地区发育着西秦岭北缘断裂、海原断裂等多条大型活动断裂带,历史上均发生过强震或大地震,此次汶川地震对这些活动断裂今后地震活动的影响不应忽视。经地震地质及大地测量资料搜集分析,印度板块以50mm/a左右的速率向北推挤,与欧亚板块汇聚于喜马拉雅构造带一线,其东部受鄂尔多斯和华南刚性块体的阻挡,在青藏高原东缘形成了一系列的挤压推覆构造带[1]。作为阻碍体之一的鄂尔多斯西南缘,在青藏高原北缘东段的祁连山—河西走廊—海原—六盘山—宝鸡构造带上,1920—1932年的短短12年内,曾发生3次7.6级以上的大地震:1920年海原8.5级、1927年古浪8.0级和1932年昌马7.6级大地震[1]。它们曾被认为是一次典型的陆内大地震沿相关断裂带由东向西迁移的例子,或是它们之间相互作用触发的结果[8]。因此本文拟采用静态位移解理论和有限元方法模拟汶川8.0级地震对宁夏及其邻近地区构造应力场造成的影响,为探讨该地区今后地震活动提供构造应力场变化的背景依据。

1 模型建立与应力场计算理论

1.1 模型建立简介

本文应用Okada[9]公式计算的解析解模型为无限半空间介质模型,用3个断层单元模拟发震破裂,3个断层单元共分为134个子断层,各个子断层面上的破裂滑动方向和滑动量分布见文献[7]。半空间或层状介质空间静态位移解的详细理论可参见文献[10-11]。断层的坐标数据来自Mapsis数据库[12]。

应用有限元方法计算的模型为国外通用力学计算软件 Ansys10.0提供的壳体模型,分别是Shell63,shell143,它们都属6个自由度壳体单元模型,能计算相对复杂的弹性变形问题,其中shell143可以进行非线性大变形计算。壳体厚度分别给定8km、10km、12km进行计算,即计算地下某层位上下一定地壳厚度内的应力应变场[13]。介质参数见表1,数据主要参考王仁[13-14]、张国民[15]等人的结果。由于内陆板块层析反演的横向速度结构间断面一般不尖锐[16-17],所以本文在稳定块体之间不是直接用断层摩擦接触算法,而采用块体两侧节点耦合算法。整个模型除少量采用3节点单元外,绝大部分采用4节点单元。整个模型各单元尺寸约5km,节点总数37471个,单元总数37813个(图1)。其中x轴为 EW向,y轴为NS向,z轴向上为正,遵守右手螺旋法则。

图1 计算区域的有限元实体模型

对于数值解的位错源的设定,根据地震破裂前锋沿断层走向的扩展情况,可以估算出朝东北方向和朝西南方向的平均破裂速度分别约为3.4和2.2km/s[18],因此本文利用王卫民等人[7]反演的结果提取出沿破裂面的位移分布标量,经震源机制解参数进行三维空间矢量投影到截面各节点的Ux,Uy,Uz分量上,断层两侧的位错权重则利用张勇等人[18]反演的两侧破裂速度比值来求取。文中把整个地震过程分5个破裂子事件,既而设置5步加载法,子载荷步数设置为15,以保证非线性求解器迭代收敛。5个子事件分别对应王卫民等人[7]反演结果的0～10s、10～40s、40～65s、65～90s、90～110s 的 5个时间段。计算中开启非线性、大变形、步长预测器等稳定迭代过程的算法[19]。

对于数值解的边界设定,研究王琪等人[20]GPS约束的位置详见图1带“”所示边界,这些边界尽可能地从稳定块体中部剖开,将剖开线上的节点6个自由度设为固定约束。而各块体之间的边界采用节点自由度耦合算法。另外,本文通过边界调整、修改材料参数、变动网格尺寸及考察辅助断层内数值变化情况,来分析边界效应的影响,同时分析边界约束的合理性。

1.2 应力场计算理论简介

反映全(半无限)空间某点的点力或位错在空间另一点上产生的力或位移的函数表达式一般称为格林函数。对于多层水平层状介质,通过Laplace及Hankel变换,构建出多层界面间物理量传递的Thomson-Haskell传递矩阵,从而实现一个反映多层介质中空间某点的点力在空间另一点上产生的力或位移的广义格林函数。应力场核心公式是 F=κ·,其中:力 F是位错的函数,κ为空间广义弹性模量(反映在格林函数中)。完整刻画线弹性空间上一点的应力状态,需要二阶应力张量。对于地震,讨论震源处的位错在地表或地下某点产生的影响,通常用索米扬那公式表示。1985年Okada[9]在索米扬那公式的基础上推导出断层位错在半空间表面产生的静态位移和应变变化的解析表达式。本文中使用的静态位移解方法按Okada的公式计算。

对于有限元法,不是去寻找求解域内的解析的格林函数,而是将原问题研究区域进行离散细分,对每一个细分块体进行原问题的等效积分,利用形函数来分段逼近原问题可能存在的格林函数,实现每个单元的有限元方程。对所有单元按空间点位进行集成,构成求解的整个基本模型,即用分片光滑的“形函数”来逼近广义弹性模量κ。有限元理论中的整体刚度矩阵与解析解法的格林函数具有相似的数理意义。观测场结果可以发现计算区域年最大形变率在10cm/a以内,扣除该区域的 GPS均值,则各观测点之间的相对形变速率更小。而大地震一般会快速释放应力能,在短暂时间(秒的量级)就产生强烈的地表位移,地震引起的形变速率远大于 GPS观测场的速率(与形变率的比值约109量级)。因此本文认为,对模型施加的边界条件,采用固定约束或采用GPS约束没有明显的差异。本文设定边界约束条件如下:(1)鄂尔多斯块体北部边界、东部边界固定约束;(2)阿拉善块体北部边界、西部边界固定约束;(3)青藏块体西部边界固定约束;(4)扬子地台东部边界、南部边界固定约束;(5)其余边界自由。块体

表1 模型分块及相应的参数设置[10,13,15,19]

2 计算结果与比较

图2 无限半空间σxy解析结果

图3 数值模拟水平剪切应力场结果

图2为半无限空间格林函数法计算得到的σxy结果,结果显示了汶川8.0级地震的右旋特征。其最大应力水平在0.1MPa的量级,这与 Toda[21]研究南极板块1998年 Mw8.1地震对余震的触发典型值0.1～0.2MPa的结果相符。图3为有限元法计算的水平剪切应力场,其在破裂面上的值稍偏高于0.1MPa的量级,明显看出在以汶川县—宝鸡市分界,汶川地震对分界线以东产生顺时针的剪应力扰动场,以西产生逆时针的剪应力扰动场。剪应力正负区间的分布可反映出地震的右旋特征,应力场在青川附近增强,且Ux4与Ux3之间的变化明显比Ux4与Ux5之间大(图4),这可能与地震破裂带由逆冲转向走滑,同时撕开较稳定的扬子地块边缘有关,这种结果符合生成新界面的应力分布特征[10]。无限半空间解析解与有限元法计算出的水平剪应力场比较,可以发现自由边界约束的解析结果的应力分布要规则得多,延伸更远。而有限元法数值模拟的结果更能刻画因介质物理属性差异,断层存在与否、边界是否受到约束等因素的影响。由于有限元法使用二维壳体模型,对王卫民等人[7]反演的三维结果进行了数据均值处理,因此计算的应力场需要与Okada三维解析方法计算的结果进行对比。从图2、图3仔细分析可看出二者在震源区附近具有一定的相似性,表现在破裂带两端的正负区间分布一致且沿破裂带两侧50km内的应力场正负分布有明显的对应性,表明数值法简化的震源模型是较可靠的。

图45 步加载法计算的Ux,Uz值分布

图5 震后1年内川甘宁地震分布

图6 汶川地震在关注断层上产生的旋转场

图7 2001年5月11日—2009年5月11日各年段LogES图像

图4为采用有限元方法,将震源发震过程分5个破裂事件而设置5步加载法得到的位移x分量各载荷步计算结果和第5步时的垂向位移Uz结果。该结果与前人远震波形反演[18]及灾害现场调查[22]的结果相符合。其中水平位移向量场分布与汶川烈度图[23]基本一致。

为了检验数值结果的可靠性,将模型的壳单元厚度设为8km、10km和12km,将居于中部的陇西地块[24]的弹性模量设为60Gpa和50Gpa,对整个模型的西、南边界沿边界法向进行约100km的向内退缩(修改边界上的关键点坐标)。按以上调整的参数分别计算,结果表明壳体厚度设置越厚,Uz分量值越小,但差值不超过5%。陇西地块的弹性模量降低后,Uz值增大,差值不超过2%,而且与前算例相比,等值线朝北侵过西秦岭断裂的范围更大。西、南边界调整,整个计算区域相对缩小,Uz值有所减小,差值不超过3%。且在边界选择上,本文设置的边界控制范围远大于相关文献[13]的建议。总体上,所有算例Uz的计算结果分布特征是基本一致的。Uz,Ux,Uy三者矢量和表现出逆冲右旋走滑特征。

图5是汶川8.0级大地震后1年内28°～40°N,100°～110°E区域内的2级以上地震震中分布。考察北川以北区域可以发现,地震空间分布有一定的成带性及丛集性,表现明显的区域有:西秦岭断裂带西段,古浪—兰州NNW向断裂系中段,海原断裂带东南段,西安东北近郊,银川—吴忠—同心一线,香山—天景山断裂局部区段,其中宝鸡北部的地震震中丛集性最为明显。进一步比较2001年5月11日到2009年5月11日的每年段地震活动分布图(图7),发现图7上标注A,B,C,D的箭头位置附近在震后一年内较前些年都有地震活动增强现象。而点E处则显示活动性减弱。图7做法是先将底图按0.05°×0.05°的网格布置零值基准面,将每一个年段的所有地震震级按Gutenberg1945年拟合的公式LogES=1.5MS+4.8转换成有物理意义的地震能量。然后按地震经纬度将能量对数值加密到零值网格中,再进行局部多项式拟合以进行成像。这种处理方式,能够在数值上连续地反映地震震级大小,反映地震空间分布的关联性,突出地震丛特点,同时弱化空间分布较独立的中强震的影响。

剪切应力场在海原断裂带东南段呈左旋性,这与该区域的断层左旋性质对应,而该区域的地震活动水平在震后明显增强,反映一定的增震效果(图7B处)。在香山—天景山断裂、海原断裂带中段上,地震空间分布不连续,计算的旋转分量场也显示其不连续性,而且局部极值(图3、图6)与地震丛集区(图7A、C、F处)有一定的对应性。而西秦岭断裂中部偏东断裂上的地震活动性不高,数值计算反映在西秦岭断裂上由西往东,旋转分量的方向几乎改变了180°,而且中东段以压性为主(第1主应力和第3主应力的结果篇幅限制不再列出)并略带右旋,这可能与该区域的原有背景应力场反向。在中卫以南附近,有旋转分量的极值,但是地震空间分布围绕中卫却不聚合或沿带状分布(图5),这种现象尚待进一步解释。而古浪—兰州NNW向断裂系西北段,数值模拟结果与小震多发性对应不强,这可能与未引入该带本底应力场有关。

综合以上分析,给出汶川地震对鄂尔多斯块体西南缘断裂体系的影响如表2。

表2 汶川大地震对关注断层的影响

3 结论与讨论

解析法和有限元数值法的计算结果都反映出了汶川8.0级地震的右旋特性及震源的复杂性,二者在场源附近具有较一致性,在宝鸡附近及鄂尔多斯西缘断裂带上都产生逆时针的旋转场。但由于介质模型及边界约束的不同,使得有限元法能反映计算区域某些局部重要特征。数值解的一些局部极值区与震后小震活跃区有对应性,而解析解未能有效反映。这可能由于震源区与关注断裂体系的空间关系过于复杂、解析模型介质简单,以及设定自由化边界有关系。如汶川大地震后,不少学者利用传统Okada方法计算库仑应力场并发表了一系列文章,石耀霖针对这种现象专门剖析了传统Okada方法的缺陷,指出地震过程中主应力的变化或边界条件变化会影响最终的同震应力场分布[25]。而地壳体内部的断层存在与否也影响着局部应力的分布,如在特例下,含隧道的山体在隧道附近及非隧道区的应力差异是非常明显的[26]。因此本文利用格林函数法和FEM方法将计算结果相互映证,同时又进一步揭示出地壳断层的引入与边界约束条件的变化对计算结果的影响。

数值结果表明,汶川地震对西秦岭断裂带的东段呈减震作用,对其西段有一定的增震作用。对西秦岭断裂中西段,在鄂尔多斯西缘断裂带南部拐角处(图7F),汶川地震导致该处附近出现部分小震,数值结果显示旋转分量很强;而西秦岭断裂东段进一步贯通,在其本体区域应力场的作用下,可能产生新的地震孕震区。海原断裂带上受到的汶川应力影响并没有过多释放(震前海原断裂带是否处于松弛状态待考虑),数值解显示高值,但汶川地震后该带地震活动性不强。在香山—天景山断裂带的局部区域造成旋转分量增强;而其它断裂带的计算结果,未显示较大量级的扰动值(表2)。在中卫附近,海原断裂带和香山—天景山断裂带的解释矛盾区值得进一步分析,而古浪—兰州NNW向断裂系西北段的地震活动可能受本底应力场影响更大。

从模型数值计算的结果来看,对于二维壳体模型Ux,Uy分量与地震波反演结果及地震地质灾害调查结果非常接近;但是Uz分量大1～2倍(不同学者[7,18,22]的结果也有1～2倍的差距),其原因一是二维模型只是三维模型的近似,特别是纵向分量上的简化;二是纵向分量除在稳定边界进行固定外,其余各点都自由,没有施加重力作用,也没有块体下边界的耦合。这两个原因都将使计算结果比实际结果偏大。壳体厚度增大可以减小Uz值。本模型的结果更适合理解成在震源深度水平面上的薄壳应力应变特征。同时,本文的研究范围限于地震释放的应力降在研究区域的作用,即扰动作用场,而不是背景场与扰动场的和。由于以地震应力降作为力源来计算,所以计算模型本底(初始)应力场—如重力场、历史形变场—可以与之分离。但现场观测的结果中,重力场实际上是起作用的,所以计算结果的垂向分量会比观测结果偏大。不过文中关注的断裂带以走滑特征为主,并与计算结果的水平分量场有较好的对应性,Z向分量则侧重模型稳定性分析及辅助解释。和前人的三维计算模型一样,本文模型对板块介质参数的设定也是尽可能利用前人地应力模拟的参数和地震层析速度结果。

在几百公里的震源附近区域,垂向形变的地质调查结果基本都在15m以内,而且除少数区域的位错场是奇异之外,其他都是较光滑的,因此利用壳体模型有其实用性。即使采用三维计算模型,由于块体内部物理参数的各向异性和不确定性,特别是青藏高原叠瓦形的奇异构造,给模型结果的准确性和可靠性带来了很多的人为因素,同时亦不易于调整参数和解释结果。在今后的深入工作中,利用GPS和历史地震构造经验性初始应力场,再在这个基础上加载大震扰动场,也许能进一步减少理论计算结果与实测结果的差异。

致谢:本文完成过程中得到赵卫明研究员和四川省地震局梁明剑工程师的多次建议和帮助,在此非常感谢。

[1] 傅征祥,刘桂萍,陈棋福.青藏高原北缘海原、古浪、昌马大地震间相互作用的动力学分析[J].地震地质,2001,23(1):35-42.

[2] 韩竹军,董绍鹏,谢富仁,等.南北地震带北部5次(1561～1920年)M≥7级地震触发关系研究[J].地球物理学报,2008,51(6):1776-1784.

[3] Harris,R.A.stress triggers,stress shadows and implications[J].Geophys Res,1998,103(B10):24347-24358.

[4] Hill,D.P.,Reasenberg,P.A.,Michael,A.,et al.Seismicity remotely triggered by the magnitude 7.3Landers,California,Earthquake method[J].Science,1993,260:1617-1623.

[5] Anderson,J.G.,Brune,J.N.,Louie,J.N.,et al.Seismicity in the Western Great basin apparently triggered by the Landers,California,earthquake,28June 1992[J].Bull.Seismol.Soc.A m.,1994,84:863-891.

[6] Bodin,P.and Gomberg,J.Triggered seismicity and deformation between the Landers,California,and Little Skull mountain,Nevada,internation[J].Bull.Seismol.Soc.A m.,1994,84:835-843.

[7] 王卫民,赵连峰,李娟,等.四川汶川8.0级地震震源过程[J].地球物理学报,2008,51(5):1403-1410.

[8] 国家地震局兰州地震研究所,宁夏回族自治区地震队.一九二○年海原大地震[M].北京:地震出版社,1980.

[9] Okada,Y.Surface deformation due to shear and tensile faults in a half space[J].Bull.Seismol.Soc.A m.,1985,75:1135-1154.

[10] 许金泉.界面力学[M].北京:科学出版社,2008.

[11] Rongjiang,Wang,et al.Computation of deformation induced by earthquakes in a multi-layered elastic crust—FORTRAN programs EDGRN/EDCMP[J].Computers&Geosciences,2003,29:195-207.

[12] 屈春燕.最新1/400万中国活动构造空间数据库的建立[J].地震地质,2008,30(1):298-304.

[13] 王仁,何国琦,殷有泉,等.华北地区地震迁移规律的数学模拟[J].地震学报,1980,2(1):32-42.

[14] 王仁,孙荀英,蔡永恩.华北地区近700年地震序列的数学模拟[J].中国科学,1982,(8):745-754.

[15] 张国民,汪素云,李丽,等.中国大陆地震震源深度及其构造意义[J].科学通报,2002,47(9):663-668.

[16] 汪素云,T.M.Hearn,许忠淮,等.中国大陆上地幔顶部Pn速度结构[C]//大陆强震机理与预测论文集,2003,(I):287-292.

[17] 李松林,张先康,任青芳,等.西吉—中卫地震测深剖面及其解释[J].地震地质,2001,23(1):86-92.

[18] 张勇,冯万鹏,许力生,等.2008年汶川大地震的时空破裂过程[J].中国科学,2008,38(10):1186-1194.

[19] 任重.ANSYS实用分析教程[M].北京:北京大学出版社,2003.

[20] 王琪,张培震,马宗晋.中国大陆现今构造变形 GPS观测数据与速度场[J].地学前缘,2002,9(2):415-429.

[21] Toda,S.and Stein,R.S.Did stress triggering cause the large off-fault aftershocks of the 25March 1998Mw=8.1Antarctic plate earthquake[J].Geophys Res lett,2000,27:2301-2304.

[22] 徐锡伟,闻学泽,叶建青,等.汶川 Ms8.0地震地表破裂带及其发震构造[J].地震地质,2008,(3):23-55.

[23] 国家汶川地震专家委员会.汶川地震灾区地震地质灾害图集[M].北京:中国地图出版社,2008.

[24] 范俊喜.鄂尔多斯地块运动特征研究[D].北京:中国地震局地质研究所,2002.

[25] 石耀霖,曹建玲.库仑应力计算及应用过程中若干问题的讨论——以汶川地震为例[J].地球物理学报,2010,53(1):102-110.

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

Numerical Simulation of Influence of the WenchuanMS8.0Earthquake upon Faults in the Southwest Margin of the Ordos Block

LIU Shuang-qing1,SHENGJu-qin2,ZENG Xian-wei2,XIE Zhao-di3
(1.Earthquake Administration of Tianjin Municipality,Tianjin 300201,China;
2.Seismological Bureau of Ningxia Hui Autonomous Region,Yinchuan 750001,China;
3.Yunnan University,Kunming 650091,China)

Based on inversed parameters of WenchuanMS8.0earthquake,using the FEM and Static Displacement Analytical Solution Theory,this paper analyzes the influence of the great earthquake on faults in southwestern margin of Ordos block.The numerical result shows that the WenchuanMS8.0earthquake,which is the type of thrust with right-lateral component,leads to clockwise disturbance displacement field to the east of the line of Wenchuan-Baoji and anticlockwise disturbance displacement field to the west,which is consist with slide direction of the active faults near the line.The distribution of horizontal shear stress extreme values is consisted with seismicity in the fault zones to the southwest of Ordos block during the year after the great earthquake.The numerical result also indicates that the Wenchuan earthquake may enhance the seismicity in the Haiyuan and Xiangshan-Tianjingshan fault zones.The region with higher shear stress value and lower seismicity may be the area that we should pay close attention to in the future.

Wenchuan earthquake;FEM;Static Displacement Analytical Solution;Ordos block

P315.0

A

1003-1375(2011)02-0057-08

2010-08-06

宁夏回族自治区科技厅攻关协作项目,浙江省科技厅协作项目(2008C13043)

刘双庆(1982-),男(汉族),广西桂林人,助理工程师,硕士研究生,现主要从事地震震源,地震烈度速报与震害评估研究.E-mail:goodmorningabc@163.com