龙门山断裂带中段及邻区构造变形特征及应力场演化的数值模拟研究*

付 平，张新辉，尹健民，韩晓玉，周春华

(长江科学院 水利部岩土力学与工程重点实验室，湖北 武汉 430010)

0 引言

岩石圈上部应力场是上覆岩层重力、大尺度力源(板块边界力、地球动力过程产生的力)和小尺度力源(岩体介质各向异性、侵入、开挖卸荷等产生的非均匀力)共同作用的结果(Zobacketal，1989)。影响地应力场的因素众多，有板块相互作用、地幔热对流、地球内应力、地心引力、地球旋转、岩浆浸入和地壳非均匀扩容等，另外，温度和水压梯度、地表剥蚀及其他物理化学变化也可引起相应的应力场变化，但其中只有重力和构造应力能达到较高量级(刘允芳等，2014；姚瑞等，2017)，即自重和地质构造作用是影响地应力场的主要因素。

龙门山是青藏高原巴颜喀拉块体和中国东部华南地块体的边界构造带，地质构造十分复杂。目前，对于龙门山地区的变形场及应力场已有不少研究成果：如白玉柱和周庆(2013)通过龙门山断裂映秀—青川段的有限元模型探讨了在垂直于断层走向的应力作用下逆冲构造运动造成的地表变形和断层附近的应力场变化；柳畅等(2012)和万永魁等(2017)通过构建岩石圈的三维粘弹性有限元模型研究了龙门山断裂带及其周缘区域的形变运动和应力累积状态；陈棋福等(2015)和尹力等(2018)分别基于粘弹性接触和二维粘弹塑性计算模拟了强震复发过程中龙门山断裂带周围地区地壳地表的变形分布与演化。上述有限元模型大多建立在断裂面是平面或为单个断裂面的假定情况下，主要研究地震复发周期内的时间演化历程，对龙门山地区叠瓦状多断裂面、内铲型逆冲运动形成的变形及应力演化的研究较少，而这些是历史持续构造运动的直接体现，值得深入研究。本文综合考虑龙门山断裂带中段及邻区的地形地貌、地壳厚度差异、深部岩体物理力学性质以及板块运动特征，基于最新地质构造剖面，建立龙门山断裂带中段的二维摩擦接触有限元模型，模拟其在构造挤压作用下的变形以及应力累积过程，探讨板块运动与构造应力之间的关系，研究断裂运动方式和地应力场演化规律。

1 研究区概况

龙门山断裂带南起泸定、天全，向东北经灌县、茂汶、北川、广元后进入陕西勉县一带，总体呈NE—SW向展布，倾向NW，长约500 km，宽30～50 km，主要为北西侧的汶川—茂县断裂和青川断裂、中央的映秀—北川断裂、南东侧的灌县—安县断裂4条近平行断裂及其控制的构造岩片组成的前展式逆冲推覆构造带(岳冲等，2017；李海兵等，2018)，其中映秀—北川断裂也称主中央断裂，如图1所示。该区域处于中国西部地质、地貌和气候的陡变带，地形起伏剧烈，高差达(3 500±500)m，莫霍面急剧变化，地壳厚度由西北方青藏高原东缘的近70 km减至东南方四川盆地的近40 km(王椿镛等，2003；Huangetal，2002)。相对稳定的华南地块阻挡了巴颜喀拉块体SE向的逃逸，使得在巴颜喀拉块体东边界产生NW—SE向或WNW—ESE向的水平挤压，造成应力在龙门山推覆构造带上高度累积，不仅导致青藏高原的持续隆升，还使其周缘地区产生强烈的地壳变形和频繁的地震活动(朱介寿等，2017)，仅2008年以来，在龙门山断裂带上就发生了2次7级以上强震：2008年5月12日汶川MS8.0地震和2013年4月20日芦山MS7.0地震。

图1 龙门山地区构造背景示意图Fig.1 Simplified tectonic setting in the Longmenshan areas

2 龙门山断裂带中段计算模型

2.1 有限元模型

巴颜喀拉块体沿先存主干断裂向SE方向逃逸，在龙门山地区受到四川盆地阻挡，其应力在断裂带上不断地累积和调整，因此，认为在这一地球动力学过程中平面应变假定近似成立。汶川地震震源机制解(郑勇等，2009；丰成君等，2018)揭示龙门山断裂带除北段表现出明显的走滑性质外，其余部位均表现为纯逆冲性质，动力源为NWW—SEE向的水平挤压，所以本文将龙门山中段构造环境简化成青藏高原向SE向推挤、受四川盆地阻挡的二维平面应变模型。选取如图1所示的垂直于龙门山断裂带走向的剖面AA′为研究对象，该剖面跨越巴颜喀拉块体、龙门山断裂带和四川盆地3个主要地质构造单元，剖面沿程地形地貌变化见图2。图2显示地形起伏集中出现在龙日坝与龙门山断裂带之间；映秀—北川断裂为龙门山地形陡变带的边界；青藏高原东缘的龙日坝断裂带具有右旋走滑兼逆冲性质，在一定程度上吸收了巴颜喀拉块体向四川盆地的挤压作用，龙泉山背斜是川西前陆盆地沉积边缘，被认为是龙门山造山带对前陆的最大影响范围。

为避免龙日坝断裂带活动的影响并体现四川盆地的阻挡作用，笔者综合上述区域地质背景及地球物理研究给出的最新深部构造特征(张新彦等，2017)，构建了如图3所示的数值计算模型，模型长度设定为160 km、深度约40 km。龙门山断裂带出露地表处位于模型中部，汶川—茂县断裂(F1)、映秀—北川断裂(F2)和灌县—安县断裂(F3)以叠瓦式分布，断裂倾角由浅至深逐渐放缓，F3和F2在深度约16 km处合并。F2向深部延伸至约23 km处与F1合并，采用四节点平面应变单元进行有限元网格划分，模型共有16 213个节点、15 709个单元，在断裂部位进行网格加密。

岩石介质的本构关系采用摩尔-库仑屈服准则描述，物理力学参数包括弹性模量、内聚力、内摩擦角、剪胀角和泊松比。龙门山地区主要发育宝兴杂岩和彭灌杂岩，该岩体抵抗破裂的强度大，能够积累较高的弹性应变能密度。计算模型中的弹性模量和密度参考花岗闪长岩的岩石力学特性试验参数(姚琪等，2012；李玉江等，2013)，并根据各单元质心位置坐标按照均匀变化的梯度形式进行赋值，以体现研究区地壳速度结构从浅部到深部表现出的“低速—高速”层状分布特点(柳畅等，2014)，其它参数为定值。将活动断裂上下两盘的相对运动看作是摩擦接触运动，在断裂面两侧设置接触对，进行断层的滑移模拟，接触对主要参数包括摩擦系数、滑移系数、法向接触刚度等。本文模型参数的设置参考了前人研究成果(朱广安等，2016；王红才等，2012；龙小刚，朱守彪，2015)，岩体和断层具体参数见表1。

图2 龙门山断裂带中段剖面AA′地形变化图Fig.2 Topography of the profile AA′ across the middle segment of Longmenshan fault zone

图3 龙门山断裂带中段及邻区有限元模型示意图Fig.3 Schematic of finite element model of the middle segment of Longmenshan fault zone and its adjacent areas

表1 模型岩体和断层参数设置Tab.1 Parameters setting of rocks and fault in the lithosphere

2.2 动力学边界条件

将龙门山断裂带所处的构造环境简化成位移约束，作为驱动计算模型的动力学边界条件。考虑到垂直于断裂带在SE方向上表现为挤压连续变形特征，选取包含剖面AA′在内的、长约600 km、宽约200 km的长方形区域，如图4a所示，将该区域内的GPS观测速度结果(Zhengetal，2017)转换为以华南地块为基准，并根据龙门山断裂带的走向分解为垂直于和平行于断裂的分量。由图4b可见，在龙日坝和龙门山断裂带附近存在2处较为明显的速度降，这是由于龙日坝和龙门山断裂先后吸收了巴颜喀拉块体的部分推挤作用，向SE跨越龙日坝断裂带后的构造推挤速度降为约2.1 mm/a。已有研究结果(Densmoreetal，2007；Zhangetal，2004)亦表明，在华南参考框架下横跨龙门山断裂带中段NW—SE方向上的水平缩短速率约1.0～4.0 mm/a。因此，计算模型的边界条件设置如下：模型右侧水平方向位移固定、垂直方向位移自由，底部水平方向位移自由、垂直方向位移固定，顶部完全自由。考虑到中下地壳与地表运动速度的差异(胡幸平等，2011)，并根据Li等(2011)利用汶川地震前数字地震台网记录的波形资料估算得出龙门山断裂带周围深部滑动速率约为GPS和地质等浅表观测滑动速率值的2倍这一结论，这里假定构造推挤速度沿深度线性变化，如图3所示，在模型左侧施加由地表2.1 mm/a增加到底部4.2 mm/a的水平向相对速度荷载。

图4 研究区相对于华南地块的GPS速度场(a)与垂直于龙门山断裂带走向的速度分量(b)Fig.4 GPS velocities relative to south China block in the study area(a)and the components perpendicular to the strike of Longmenshan fault zone(b)

3 龙门山断裂带中段及邻区变形特征

合理地给定初始状态是地球动力学数值模拟的关键前提条件(祝爱玉等，2016)，本文模型的初始应力场为静岩应力场，模型中任意一点各方向的正应力等于其上覆岩体的重力，且无剪应力，具体如下：施加重力后，再经自重应力场平衡算法消除掉因自重产生的节点位移，在这个应力场近似为自重应力、位移场近似为零的初始状态下对模型进行速度加载。从应力计算结果看，历经约67万年的持续运动加载，应力变化不显著，基本趋于稳定，模型进入到一个相对稳定的构造状态。

速度加载至67万年，垂直位移分布如图5所示，汶川—茂县断裂、映秀—北川断裂、灌县—安县断裂出露地表处上下两盘垂直方向上的相对位移演化曲线如图6所示，F1，F2，F3的垂直方向相对滑移量分别为72.5，271.2，186.4 m，由此得出其平均逆冲速率分别为0.11，0.40，0.28 mm/a。GPS观测和地质年代学研究(马保起等，2005；李勇等，2006；张培震等，2008)显示龙门山断裂带的逆冲速率如下：汶川—茂县断裂为0.3～0.8 mm/a、映秀—北川断裂为0.4～1.2 mm/a、灌县—安县断裂为0.3 mm/a。该结果表明，数值模拟得到的运动速率与实际观测结果基本一致。

图7显示这3条活动断裂的上盘抬升速率均大于下盘，汶川—茂县断裂左侧的巴颜喀拉块体和龙门山断裂带区域呈水平缩短和不断抬升趋势，灌县—安县断裂右侧的四川盆地在垂直向则保持相对稳定。这3条断裂的逆冲速率关系为vF2＞vF3＞vF1，汶川震后调查表明(何宏林等，2008；张培震，2008)，沿汶川—茂县没有发生破裂，而在映秀—北川—平通—南坝一带产生了严重的地表破裂和错动，沿灌县—安县也形成了较长的地表破裂，可以看到，本文模拟结果支持了上述地震地质研究结论。

已有研究认为龙门山地区自西向东形成系列的逆冲断层和推覆构造，逆冲断层将青藏高原东部地壳向东逃逸转换为隆升(Xuetal，2009)，中央的映秀—北川断裂的强逆冲性使此断裂长期发震从而产生较大的垂直位移累积量，尽管与同震累积垂直位移量和地层剥蚀量相比，持续逆冲作用所产生的垂直位移量较小，但仍然表明映秀—北川断裂对龙门山地区的隆升起主要贡献。根据Slemmons和Depolo(1986)提出的活动断裂分类方案，龙门山断裂属于中低活动速率断裂带，而该区域地震活动频繁且强度大，这也从侧面说明：虽然滑动速率是衡量断裂活动性的一个重要参数，但仅依靠地表调查和震前GPS观测得出的活动速率评估活动断裂区的地震危险程度具有一定的局限性，还需结合断裂深部动力学过程和区域地应力状态开展系统分析研究。

图5 龙门山断裂带中段及邻区垂直位移云图

Fig.5 Vertical displacement distribution of the middle segment of Longmenshan fault zone and its adjacent areas

图6 龙门山断裂带中段逆冲位移演化曲线Fig.6 Thrust displacement curve of the middle segment of the Longmenshan fault zone

图7 地表水平向和垂直向变形速率Fig.7 Horizontal and vertical deformation rates of the surface

4 构造运动过程中应力场演化规律

研究人员在龙门山断裂带区域开展了包括应力解除、水压致裂等方法在内的大量原地应力测试工作(安其美等，2004；陈群策等，2012；秦向辉等，2013)。受汶川地震同震应力效应影响，汶川地震后龙门山断裂带实测地应力沿断裂走向具有明显分段特征，本文选取断裂带中段剖面AA′附近的地应力实测点(测试钻孔位置见图1)，整理其中岩心相对完整的深部测试段获得的主应力数据列于表2。为探讨区域构造运动对地应力场的影响，在有限元模型F1，F2，F3这3条断裂上下两盘的不同深度位置上设置若干组主应力监测点，提取各监测点主应力量值列于表3，其中，中间主应力为面外主应力。为检验模拟结果的可靠性，将监测点的特征值主应力比λ=σ1/σ3与实测点的特征值kHv=σH/σv进行对比，λ和kHv均表征最大水平主应力作用强度。由表2，3可以看到，模型监测点λ的范围为1.88～3.17，且基本上呈现深度越深比值越小的特点；实测点kHv的范围为1.39～3.03，同样呈现上述规律，但数值离散性较强，原因在于浅部应力易受岩体结构的影响，总体上看，模拟计算得出的特征值与实测结果大致相符。

表2 龙门山断裂带中段及邻区地应力实测结果Tab.2 Results of in-situ stress measurement in the middle segment of Longmenshan fault zone and its adjacent areas

表3 模型各监测点主应力量值Tab.3 Principal stresses of monitoring points in the simulation model

利用2008年汶川MS8.0主震及余震序列龙门山中段地区的震源机制解(Global CMT，2013；崔效锋等，2011)，对该地区的构造应力场进行反演，得到了主应力方向和应力形因子R的最优解(图8)。图8显示，最大主应力轴方位角100.7°，倾角2.6°；中间主应力轴方位角9.8°，倾角20.8°；最小主应力轴方位角196.5°，倾角70.6°，R为0.72，应力结构为逆冲型，显示该地区主要受近NWW向的水平挤压应力场控制。本文所构建二维平面模型得到的最大主应力方向为剖面的取向NWW向，这里统计了3条断裂监测点处最大主应力倾角的平均值，F1上下盘分别为3.2°和7.3°，F2上下盘分别为11.1°和8.8°，F3上下盘分别为12.3°和2.7°，均为小角度近水平，可以看到数值模拟得到的主应力方向和倾角与震源机制解反演结果较接近。应力形因子R=(σ1-σ2)/(σ1-σ3)反映了3个主应力的相对大小，表3中各监测点的R值大于反演得到的最佳值，原因在于本文对构造环境的平面应变简化使得沿断裂带走向的变形受到强约束，未能考虑NNE向中间主应力产生的压张效应。

图8 汶川地震序列龙门山中段地区震源机制解反演结果

Fig.8 Inversion results of focal mechanism solutions of Wenchuan earthquake sequence in the middle segment of Longmenshan fault zone

由图9所示的断裂带附近监测点主应力随时间演化曲线可以看出，随着持续的区域构造运动，主应力基本呈不断增大趋势，符合挤压逆冲型构造环境，主应力量值的变化过程大致可分为3个阶段：较低水平缓慢增加的低态稳定阶段、持续增加的应力累积阶段和达到临界点后小幅度变化的临界稳定阶段。历经10～15万年的低态稳定阶段构造运动后，各监测点主应力量值开始显著增加，但进一步达到临界状态所需时间不同，深部要晚于浅部进入这个临界相对稳定阶段。F3断层6.5 km、F1断层6.7 km、F2断层8.1 km处监测点达到临界稳定状态的时间约为45～50万年，F1断层15.2 km处监测点达到临界稳定状态的时间约为60万年，F1断层20.4 km处监测点则看不到明显的稳定阶段，60万年后主应力量值仍表现出降低—升高—降低的波动性，说明当活动断裂带系统的构造应力达到一个临界值后，应力将处于不断调整与重新分配的过程，应力的累积和释放会形成一个相对稳定的自平衡状态，且深度越深，这种自发式的应力调整活动越频繁、越剧烈。临界稳定阶段出现的应力波动被认为是震前的自发式调整，是发震的必要前提条件，若该阶段岩体在无法承受更高水平的应力累积时快速破裂，应力迅速降低，这种长期累积的应变能突然释放即预示着地震的发生。已有研究(闻学泽，2018；孙云强，罗纲，2018)认为断裂构造带在经过长时间的应力累积期后，将进入能量释放的地震丛集期，而在发震后应力又会以一定速度继续累积，2008年汶川、2013年芦山以及2017年九寨沟地震事件也初步证明龙门山地区处于应力累积—地震发生—应力再累积的地震周期。由此看出，数值模拟得出的低态稳定—应力累积—临界稳定的阶段式构造应力演化过程与该地区现实中的地震周期性活动特征相呼应。

图9 不同深度处断裂带附近主应力演化曲线Fig.9 Evolution of principal stress near faults at different depths

5 结论与讨论

本文通过对龙门山断裂带中段及邻区构造变形以及应力场演化进行限元数值模拟，并结合已有研究成果，得出如下结论：

(1)龙门山断裂带属于挤压逆冲型断裂，汶川—茂县、映秀—北川、灌县—安县断裂的逆冲速率分别约0.11，0.40，0.28 mm/a，与实际观测结果基本一致，中央的映秀—北川断裂逆冲性质最为强烈，对龙门山隆升起主要贡献，龙门山左侧的巴颜喀拉块体和断裂带区域呈不断抬升趋势，右侧的四川盆地保持相对稳定。

(2)数值模拟结果显示3条主要断裂附近6.5～20.4 km深度范围内岩体大小主应力比为1.88～3.17，且表现出深度越深主应力比越小的规律，总体上与实际测试结果相符，数值模拟获得的主应力方向和倾角与利用该地区震源机制解反演得到的构造应力场特征接近。

(3)龙门山断裂带中段区域构造应力呈现低态稳定—应力累积—临界稳定的阶段式演化过程，且深部要晚于浅部达到临界稳定状态，其临界稳定阶段的应力值仍表现出降低—升高—降低的波动性，这种阶段式应力演化过程与该区域所处的应力累积—地震发生—应力再累积的地震活动特征相呼应。

(4)通过地表调查和震前GPS观测揭示的断裂活动速率分析评价活动断裂区的地震危险程度存在一定局限，需要进一步结合其深部动力学过程和区域应力状态进行综合研究。

本文侧重在研究思路与方法上，通过垂直于龙门山断裂中段走向的二维剖面模型，探讨构造挤压作用下多断裂面逆冲运动过程中的变形及应力演化规律。数值模拟尚存在一些待改进之处，例如地震的触发和破裂过程未加以考虑，文中二维模型虽然是三维情形的合理简化，但该地区地质构造十分复杂，可考虑更多因素建立三维模型。在三维接触或流变模型的基础上尝试引入合理的损伤及愈合机制，研究震前构造应力累积—震时应力释放与调整—震后应力恢复再累积全过程的应力状态变化，对于科学认识孕震发震机理是很有意义的。

感谢评审专家提出的宝贵意见和建议，感谢中国地震局第一监测中心提供的GPS观测数据。