1976年唐山大震发生对华北地区各地块运动与变形影响的数值模拟研究

用非连续变形分析方法(DDA+FEM)数值模拟在华北地区各地块相互制约的构造环境中,发生1976年唐山地震的动力学过程;研究大震引起的华北地区各地块特别是鄂尔多斯地块运动变形及各地块边界断层上应力状态变化的特征。采用的1976年唐山地震复杂震源模型由3个子事件组成:北北东走向的主破裂,与主破裂北东端部相连接的北东走向子事件以及位于主破裂南东侧的北西走向的子事件,同前人的结果比较,这个模型与该震早期余震分布图像的位置更为相符。计算得到了地震释放的主应力场,最大剪应力变化等值线图,大震邻近地区地块滑移变形图,发震断层滑移距离随时间的变化以及最大走滑错距和应力降等震源参数,这些结果与该震的震源机制解,震源参数,宏观等震线,地表观测的水平位移矢量图基本一致。唐山地震的发生使得与其邻近的地块及边界断层受到不同程度的影响,发生了不同规模的移动和变形,鄂尔多斯地块似乎发生了整体性小规模的逆时针转动及变形,其中鄂尔多斯地块东北部与太行山,燕山,阴山地块交接带的相对运动变形较大,唐山地震的发生也引起了华北地区各地块边界断层上的应力状态发生了一定程度的变化,其中通过唐山震区的张家口-蓬莱北西走向断裂带以及鄂尔多斯地块北部边界断层的主要地带上剪应力增加而法向应力减小,使这些断裂带的地震危险度增加。这与唐山地震后鄂尔多斯地块北部边界断层上发生一系列6级强震的事实基本相符。     

用（DDA+FEM）方法数值模拟1975年海城、1999年岫岩地震发生的过程l

用非连续变形分析方法（DDA+FEM）， 数值模拟了华北地区各地块相互制约构造环境中发生的1975年海城地震及1999年岫岩地震的过程. 通过对岫岩前震及近期余震分布图象的研究， 以及前人对海城地震的研究， 提出了海城、 岫岩地震发震构造块体模型. 数值模拟结果给出了两次地震释放的主应力变化、 最大剪应力变化等值线图， 地震前后位移变化矢量图， 及发震断层滑移随时间的变化， 分别与相应地震的震源机制、 宏观等震线、 发震断层的走滑性质等的结果基本一致.      

震源动力学破裂过程数值模拟研究

首先阐述了震源动力学过程研究的重要意义, 在此基础上, 研究了复杂的断层几何形态及介质模型对动力学破裂过程的影响, 并对常用的有限元方法、 离散元方法、 有限差分方法和边界积分方程方法等进行了相应介绍. 讨论了这些数值模拟方法各自的优缺点, 建议在方法的选择上应视具体问题及计算的精度而定. 最后对动力学数值模拟的关键部分, 滑动摩擦准则进行了论述. 常用的滑动摩擦准则有滑动弱化准则、 速率弱化准则和速率-状态依赖摩擦准则. 在单纯考虑某个地震的动力学破裂传播过程时, 滑动弱化准则较为常用, 其中滑动弱化距离的选取至关重要. 但若考虑整个地震循环, 速率-状态依赖摩擦准则更为合适.      

1927年古浪8级大震破裂的三维数值理论模拟

在数值模拟计算的基础上,进一步从理论上论证了１９２７年古浪８级大震地表破裂带较短的原因．通过对地震烈度分布形态,地表破裂受阻的终止过程,理论模拟结果所显示的位移矢量的衰减过程,断裂受阻变形及最大剪应力分布特性等问题的研究,认为分布在主破裂带东西两端部的ＮＮＷ向断裂起到了“障碍构造”的作用,阻止了主破裂向东西方向的延伸扩展．     

2001年8.1级昆仑山大震破裂过程及对2008年汶川8.0级大震孕育发生影响的研究

本文用三维流变非连续变形(块体边界)与有限元(块体内)相结合(DDA+FEM)的方法,在青藏高原及其东侧四川盆地,鄂尔多斯块体地区三维构造块体相互制约的大环境中,考虑了龙门山断裂带东西两侧地势、地壳厚度和分层的明显变化,及断裂带东侧四川盆地及鄂尔多斯块体坚硬地壳阻挡的影响,通过用GPS资料做位移速率边界约束和震源机制约束,计算得到研究区的速度场和应力场与该地区GPS测量结果和震源机制分布结果基本相似.在此基础上,数值模拟2001年昆仑山大震的破裂过程;研究大震引起各构造块体边界断层应力状态变化特征,特别是对2008年汶川大震发震断层的影响.结果表明:(1)数值模拟昆仑山大震发震断层发生左旋走滑错动,最大水平错距约4.5 m,最大应力降约18 MPa.计算获得大震释放的主压应力场图像,最大剪应力变化等值线图,大震发震断层垂直面上位错等值线图及大震引起垂直位移变化三维图分别与大震的震源机制,地表破裂带同震位移分布,GPS同震位移图及地震波反演和GPS反演的结果总体上均比较相近.(2)计算获得的最大剪应力变化等值线图分布具有不对称特点,大震发震断层南侧变化梯度明显大于北侧.(3)模拟计算得到大震引起汶川大震发震断层库仑破裂应力增加约0.016 MPa(上地壳层).昆仑山大震破裂过程是在东昆仑断裂带其发震断层上发生的左旋走滑错动,引起东昆仑断裂带南侧巴彦喀拉块体进一步东扩和一定规模的变形,并受到该块体东侧四川盆地较硬地壳的阻挡,使得块体东边界断层中低倾角的汶川大震发震断层库仑破裂应力增大,应变能积累增强.可以认为这一破裂过程对汶川大震发震断层发生逆冲型失稳起了促进作用.     

1997年玛尼地震对青藏川滇地区构造块体系统稳定性影响的三维DDA+FEM方法数值模拟

本文用三维非连续变形与有限元相结合（DDA+FEM）的方法,在青藏川滇地区三维构造块体相互制约的大背景中,通过用GPS资料做位移速率边界约束和震源机制约束,计算得到研究区的初始位移场和应力场与该地区GPS测量结果和震源机制分布结果基本一致.在此基础上进一步数值模拟1997年玛尼7.9级大震的发生过程,研究大震引起研究区各块体边界断层应力状态变化的特征.（1）发震断层两侧发生左旋走滑错动,最大水平位错大约7 m;（2）深部位错面上位错分布与用地震波资料震源反演的结果类似;（3）最大差应力变化等值线图与由星载D\|INSAR技术获取的地表形变场图像相似;（4）地表垂直位移表明地震断层面略向北逆冲.计算模拟得到了玛尼地震发生引起青藏川滇地区构造块体系统各边界断层上库仑破裂应力变化的分布,表明玛尼大震的发生除了使其发震断层的两端库仑破裂应力增大,应力进一步集中外,位于上地壳层上东昆仑断裂中段的2001年昆仑山8.1大震（H=11 km）发震断层段的库仑破裂应力增加约2 MPa,位于中地壳层上喀拉昆仑断裂带中的2008年改则6.9级地震（H=30 km）发震断层段的库仑破裂应力也增加约0.7 MPa,可见这两个已接近破裂强度地段的失稳对发生大震起了一定促进作用.研究结果也表明：作者发展的三维DDA+FEM方法能有效地用于大震活动与各构造块体相互作用关系的研究.     

含预制软弱带的岩石破裂过程的数值模拟

考虑到岩石脆性破坏过程中介质的不均匀性之特点，对岩石样品中含预制弱介质条带的岩石样品破坏过程进行二维有限元数值模拟，并对弱介质带的破坏贯通过程、新的断层的产生和有关的地震活动进行了研究．数值模型展示了岩石从变形、微观破坏到整体破坏的全过程以及微震活动的时空分布特征．应力、应变和微震活动的时空分布形象地描绘了岩石变形的局部化和时空迁移等现象，这与实际地壳中所观测到的现象是一致的．此外，模拟结果与实验观测结果也是一致或相似的．      

模拟震源动力学破裂过程的数值方法研究

为模拟震源的动力学破裂过程,本文讨论一种求解剪切裂纹动力学扩展问题的有限差分法。在研究二维反平面破裂的基础上,我们把所用方法推广到三维问题。研究了许多简单震源模型,通过对比同一问题的解析解和其它数值解,对方法的正确性和计算精度进行了检验,结果表明我们的数值方法是可行的。最后,作为初步应用的例子,我们研究了非均匀断层的自发破裂问题。     

模拟震源动力学破裂过程的数值方法研究

为模拟震源的动力学破裂过程,本文讨论一种求解剪切裂纹动力学扩展问题的有限差分法。在研究二维反平面破裂的基础上,我们把所用方法推广到三维问题。研究了许多简单震源模型,通过对比同一问题的解析解和其它数值解,对方法的正确性和计算精度进行了检验,结果表明我们的数值方法是可行的。最后,作为初步应用的例子,我们研究了非均匀断层的自发破裂问题。     

地震剪切破裂的三维断裂力学有限元分析—1966年邢台地震的数值模拟

以往用断裂力学方法研究地震的发生，绝大多数是针对二维问题或三维椭圆盘断层，而真实的断层是三维的。本文运用线性断裂力学理论，以Ｋｏｓｔｒｏｖ的Ｇｒｉｆｆｉｔｈ静肪表面能为断裂准则，对地震剪切破裂进行了三维有限元分析。在断裂面前缘，我们采用了由２０点等参单元退化而成的１５点元，非奇异单元用２０点元，证明了１５点元能反映应力的奇异性，推导了局部会标系下的应力强度因子的有限元解，并编制了大型的三维线弹性断     

1927年古浪8级地震的破裂习性及破裂机制的数值模拟

1927年古浪8级地震是继1920年海原8．5级地震之后发生在海原一祁连山断裂带上的另一次特大地震。对这次地震的破裂习性，前人的研究结果存在较大分歧。本文在多探槽揭露和地面追踪调查基础上，结合相关的历史资料分析后，认为本次地震是由构成古浪推覆体的天桥沟一黄羊川断裂、皇城双塔断裂冬青顶断裂段以及武威一天祝隐伏断裂等共同作用的结果。对古浪推覆体平面和剖面变形机制的有限元数值模拟结果显示，其应力和应变的集中区主要分布在天桥沟一黄羊川断裂西段、皇城一双塔断裂冬青顶一带以及武威一天祝隐伏断裂和古浪一双塔断裂所在的古浪峡一带，这与地震地表破裂带的展布是一致的，同时也进一步说明了1927年古浪8级地震是该推覆体整体活动的结果。     

唐山大地震孕育过程的数值模拟研究

根据1976年唐山7.8级地震的前兆资料，结合岩石破裂物理、震源物理、构造物理、地震前兆、地震力源、地壳结构等方面的研究结果，采用含有菱形硬包体的饱和水孔隙介质模型，数值模拟了唐山大地震的孕育过程．首先导出了考虑岩石非弹性体积膨胀和应变软化、孔隙水渗透的固-液两相介质的孕震动力学方程组，提出了用增量形式求解该方程组的有限差分法，然后计算模拟了硬包体内含有断层软弱带时的唐山大震的孕育过程，及其时空变化特征．数值模拟结果表明，唐山大地震的孕育主要经历了弹性积累、早期非弹性体积膨胀、断层第1次加速软化(加速破坏和蠕滑)、天津-宁河一带第2次非弹性膨胀、断层第2次加速软化(临震加速破坏和蠕滑)和发震等过程．这里，与断层的加速破坏相对应，孕震介质的某些部分出现了非弹性膨胀的弹性恢复和二次非弹性膨胀．与前兆观测资料比较分析可以发现，上述模拟结果可以比较好地解释实际观测资料的变化特征．同时，数值模拟还进一步揭示了唐山地震孕育场的时空变化规律和复杂图象，为地震前兆场变化的物理机制的研究奠定了基础．      

1976年唐山地震震源动力过程的数值模拟

用新LDDA(Lagrangian DiscontinuousDeformation Analysis)方法模拟了唐山地震断层的破裂、错动和应力释放的整个动力过程．模拟结果表明,唐山地震的震源滑动过程在发震断层上各处不一样．近场位移受断层的曲率影响,断层凹侧的位移大于断层凸侧的位移．滑动过冲现象在震中处最大,并向断层两端衰减．我们发现,唐山地震断层的破裂速度和应力降与断层上的初始剪应力大小有关．唐山发震断层的最大动态、准静态位错量和剪应力降均发生在中间部位,分别是7.1 m、6.2 m和8.1 MPa、5.4 MPa,发震断层的平均准静态位错量和剪应力降分别为4.5 m和3.3 MPa,断层破裂的传播速度从震中向东南和西北方向分别为3.08 km/s和1.18 km/s．      

三维数值模拟在华北地区现今构造变形分析中的应用

基于并行版ANSYS数值模拟软件平台,根据华北地区活动地块划分及活动断裂分布,结合GPS资料等,确定模型的几何边界范围是99.8°~121.4°E,27.9°~42.3°N,模型的边界范围包括华北活动地块区的绝大部分以及周边青藏、西域、南华和东北亚地块的部分地区。模型单元平均边长为25km,模型共划分为单元416 582个,节点582 392个。主要研究成果为:1)地壳运动速度模拟与分析对比:华北地区地壳运动速率的整体分布呈现由东向西逐渐减小、自北向南逐渐增大的特点,数值模拟结果与GPS观测运动速度大部分较吻合;2)断层运动性质模拟与分析对比:华北地区有限元模型考虑了华北地区已知所有晚更新世以来的活动断裂,模拟结果所得的断层运动性质与活动断裂地质调查的运动性质符合得较好;3)应变场模拟与分析对比:通过数值模拟,计算了1999—2004年和2004—2007年华北地区地表第1主应变和第3主应变(最大压应变)的大小和方向,与前人利用震源机制结果、GPS测量结果等资料反演得到的华北地区主应变场较一致。     

1906年玛纳斯7.7级地震孕震过程的数值模拟研究

根据玛纳斯 7 7级强震区的构造格架 ,利用三维有限元方法对 190 6年玛纳斯 7 7级地震的强震区进行了应力模拟研究。结果表明 ,玛纳斯地震的震源区不仅位于应力集中的部位 ,而且该部位集中的应力还会随着时间的增长而增大。数值模拟较好地再现了 190 6年玛纳斯 7 7级地震震源区的应力积累过程     

松潘—甘孜地块地壳和上地幔顶部现今变形模式数值模拟研究

松潘—甘孜地块至四川盆地是整个青藏高原地形梯度最大的地区,自西向东30~50 km范围内存在约4 km的显著地形高差,这可能与松潘—甘孜地块埋深20~25 km处的低阻低速层、地壳至上地幔顶部水平向运动速率随深度增加以及地幔对流拖曳力均有关联.为获得松潘—甘孜地块地壳和上地幔顶部现今变形模式,定量评价低阻低速层、中下地壳和上地幔顶部变形对川西高原隆升的影响.本文依据跨龙门山断裂带探测剖面,构建了长300 km、宽104 km的二维有限元接触模型,以1991—2016年GPS观测结果为初始约束,考虑岩石在长时间载荷蠕变作用的前提下,通过改变低阻低速层蠕变参数、模型西边界和底边界加载条件,共计开展了11项数值模拟实验.将模拟结果与综合多学科、不同时间尺度获得的地表隆升速率进行对比,各项模拟结果横向对比,对松潘—甘孜地块地壳至上地幔顶部现今变形及物质运移模式展开讨论,主要结论如下：（1）松潘—甘孜地块长期稳定、更为客观的最大隆升速率可能为1.4~1.5 mm·a-1;（2）低阻低速层可以作为中下地壳和上地幔物质运移的一个滑移面,构成上地壳与中下地壳的"解耦"带,促进松潘—甘孜块体整体隆升,但影响有限,不足以形成如此规模的强烈的地表抬升;（3）同时考虑低阻低速层作用及模型西边界加载速率随深度线性增加1.8倍的前提下,模拟的地表隆升速率与综合多学科、不同时间尺度获得的隆升速率最为吻合,这一动力条件下,中下地壳和上地幔物质受四川盆地强烈阻挡,速率快速减小分布于松潘—甘孜地块内部（模型50~150 km范围）,松潘—甘孜地块东缘和龙门山断裂带速率减小已不显著,故地表隆升中心位于川西高原而非龙门山断裂带;（4）青藏高原东缘物质汇聚到一定程度后,会迫使部分物质向两侧运动,实现物质运移动态平衡.     

地幔柱与岩石圈相互作用过程的数值模拟

地幔柱的研究是地球科学研究的热点之一.本文主要集中研究地幔柱与岩石圈的相互作用过程.基于质量守恒方程、动量守恒方程和能量守恒方程,通过有限元数值方法可以计算得到地幔柱与岩石圈相互作用的温度场、速度场和有效黏度等的时空图.本文的流变本构模型主要基于非牛顿流体的有效黏度模型,通过数值模拟计算分析了地幔柱与岩石圈相互作用过程,着重讨论了地壳流变结构对此过程的影响.数值模拟结果显示,地幔柱与岩石圈的相互作用分为三个阶段:地幔柱上升期,时间持续到0.2 Ma,平均速度为2.75 m·a-1,地幔柱顶部地形开始向上隆起;地幔柱与岩石圈纵向作用期,时间从0.2 Ma到0.26 Ma,地幔柱上升的平均速度为0.83 m·a-1,地表地形隆升达到最大值;地幔柱与岩石圈横向作用期,0.26 Ma以后,岩石圈开始剪切变形,地幔柱水平运动速度为0.47 m·a-1,当剪切变形达到一定程度,岩石圈底部开始出现拆沉作用.当下地壳流变强度比较小时,上地壳的流变结构控制着地幔柱顶部地表地形隆起程度,流变强度越大,隆升高度越小;而下地壳的流变结构控制着地幔柱两侧地表地形的下沉幅度,下地壳流变强度越小,下沉幅度越大.最后,讨论了数值模拟对峨眉山大火成岩省地幔柱发展演化的应用.     

唐山地震断层破坏及其力学过程的数值模拟

本文用数值分析方法探讨了地震从孕震到发震的力学过程,以76年的唐山地震为实例,用非线性有限元方法,分析唐山发震断层破坏的力学机制,用增量分析的方法,近似地按时间序列模拟发震断层和周围岩石从孕震到发震的力学过程,并与实测资料进行分析比较     

2008年汶川8.0级地震孕震机理研究

用三维流变非连续变形与有限元相结合（DDA+FEM）的方法,在青藏高原及邻近地区三维构造块体相互制约的大背景中,考虑了龙门山断裂带东西两侧地势、地壳厚度和分层的明显变化,及断裂带东侧四川盆地及鄂尔多斯块体坚硬地壳阻挡的影响,通过用GPS资料做位移速率边界约束和震源机制约束,计算得到研究区的速度场和应力场与该地区GPS测量结果和震源机制分布结果基本一致.在此基础上,模拟计算现今构造块体边界断层上表征剪应力及法向应力等综合影响的危险度分布.结果表明,上、中地壳层危险度分布中危险度较高的地段多数与近几十年来发生的七级以上大震区域基本一致.包括2008年汶川8.0级等大震的发震断层.通过分别对龙门山断裂带东西两侧的两种不同构造格局进行试算表明,龙门山断裂带东西两侧地势、地壳厚度、分层与物性明显变化对汶川大震的孕育发生均起了关键性作用.计算得到的应变率强度分布图可见,高原东部整个边缘地带均接近应变率强度的陡变带.其中以龙门山断裂带上的陡变最为明显,西侧应变率强度是东侧的近4倍,而且断裂带东侧应变率强度等值线衰减比西侧快.反映了汶川大震逆冲型发震断层地区独特的特征.此外,由计算得到的应变能密度分布图可见,龙门山断裂带在上、中地壳层中均位于宽度相同、其走向与龙门山断裂带走向一致的高应变能密度带中,在上地壳层这个带的东西两侧则是应变能密度较低的地区,而在中地壳层,其强度在断裂带东侧逐渐向东衰减,西侧应变能密度高,而东侧应变能密度较低.表明在印度板块强烈推挤作用和高原各构造块体相互制约及龙门山断裂带东西两侧特殊构造环境中,高原地壳物质向东水平运动,受到龙门山断裂带东侧介质刚性强度较大的四川盆地阻挡,使得汶川大震发震断层在大震前已积累了相当水平的应变能,并同时处于力学上的不稳定状态.