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

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

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

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

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

速度加载至67万年,垂直位移分布如图5所示,汶川—茂县断裂、映秀—北川断裂、灌县—安县断裂出露地表处上下两盘垂直方向上的相对位移演化曲线如图6所示,F₁,F₂,F₃的垂直方向相对滑移量分别为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条断裂的逆冲速率关系为v_F2>v_F3>v_F1,汶川震后调查表明(何宏林等,2008; 张培震,2008),沿汶川—茂县没有发生破裂,而在映秀—北川—平通—南坝一带产生了严重的地表破裂和错动,沿灌县—安县也形成了较长的地表破裂,可以看到,本文模拟结果支持了上述地震地质研究结论。

已有研究认为龙门山地区自西向东形成系列的逆冲断层和推覆构造,逆冲断层将青藏高原东部地壳向东逃逸转换为隆升(Xu et al,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

研究人员在龙门山断裂带区域开展了包括应力解除、水压致裂等方法在内的大量原地应力测试工作(安其美等,2004; 陈群策等,2012; 秦向辉等,2013)。受汶川地震同震应力效应影响,汶川地震后龙门山断裂带实测地应力沿断裂走向具有明显分段特征,本文选取断裂带中段剖面AA'附近的地应力实测点(测试钻孔位置见图1),整理其中岩心相对完整的深部测试段获得的主应力数据列于表2。为探讨区域构造运动对地应力场的影响,在有限元模型F₁,F₂,F₃这3条断裂上下两盘的不同深度位置上设置若干组主应力监测点,提取各监测点主应力量值列于表3,其中,中间主应力为面外主应力。为检验模拟结果的可靠性,将监测点的特征值主应力比λ=σ₁/σ₃ 与实测点的特征值k_Hv=σ_H/σ_v 进行对比,λ和k_Hv均表征最大水平主应力作用强度。由表2,3可以看到,模型监测点λ的范围为1.88～3.17,且基本上呈现深度越深比值越小的特点; 实测点k_Hv的范围为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年汶川M_S8.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条断裂监测点处最大主应力倾角的平均值,F₁上下盘分别为3.2°和7.3°,F₂上下盘分别为11.1°和8.8°,F₃上下盘分别为12.3°和2.7°,均为小角度近水平,可以看到数值模拟得到的主应力方向和倾角与震源机制解反演结果较接近。应力形因子R=(σ₁-σ₂)/(σ₁-σ₃)反映了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万年的低态稳定阶段构造运动后,各监测点主应力量值开始显著增加,但进一步达到临界状态所需时间不同,深部要晚于浅部进入这个临界相对稳定阶段。F₃断层6.5 km、F₁断层6.7 km、F₂断层8.1 km处监测点达到临界稳定状态的时间约为45～50万年,F₁断层15.2 km处监测点达到临界稳定状态的时间约为60万年,F₁断层20.4 km处监测点则看不到明显的稳定阶段,60万年后主应力量值仍表现出降低—升高—降低的波动性,说明当活动断裂带系统的构造应力达到一个临界值后,应力将处于不断调整与重新分配的过程,应力的累积和释放会形成一个相对稳定的自平衡状态,且深度越深,这种自发式的应力调整活动越频繁、越剧烈。临界稳定阶段出现的应力波动被认为是震前的自发式调整,是发震的必要前提条件,若该阶段岩体在无法承受更高水平的应力累积时快速破裂,应力迅速降低,这种长期累积的应变能突然释放即预示着地震的发生。已有研究(闻学泽,2018; 孙云强,罗纲,2018)认为断裂构造带在经过长时间的应力累积期后,将进入能量释放的地震丛集期,而在发震后应力又会以一定速度继续累积,2008年汶川、2013年芦山以及2017年九寨沟地震事件也初步证明龙门山地区处于应力累积—地震发生—应力再累积的地震周期。由此看出,数值模拟得出的低态稳定—应力累积—临界稳定的阶段式构造应力演化过程与该地区现实中的地震周期性活动特征相呼应。

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

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

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

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

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

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

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

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