基于离散元法的龙门山构造带隆升变形主控因素研究

王迎 李江海,† 马昌明 宋珏琛

基于离散元法的龙门山构造带隆升变形主控因素研究

王迎1,2李江海1,2,†马昌明1,2宋珏琛1,2

1.造山带与地壳演化教育部重点实验室, 北京大学地球与空间科学学院, 北京 100871; 2.北京大学石油与天然气研究中心, 北京 100871; †通信作者, E-mail: jhli@pku.edu.cn

为探明龙门山构造带隆升变形的主要控制因素, 基于龙门山构造带东西两侧下地壳物质层属性差异巨大的特征, 进行 3 组 PFC2D 离散元数值模拟对比实验, 将深度扩大至下地壳, 记录颗粒运动状态, 实现定量化分析。实验得到的变形结果及模型颗粒运动矢量图显示, 在下地壳物质属性无明显差异的条件下, 板块碰撞挤压应力及地壳厚度的差异不会在龙门山构造带形成巨大的地形高差。当下地壳黏度系数存在明显差异时, 软弱下地壳物质层颗粒相对运动速率为 1.5～2.94m/s, 平均运动速率为 1.62m/s, 大约是坚硬下地壳层颗粒平均运动速率的 54 倍。模型中部(龙门山构造带)出现隆升变形, 纵向影响范围为 94.74%, 隆升幅度为19.85%。软弱下地壳上覆的中地壳和上地壳颗粒具有较大的向上速度分量, 上地壳物质层上涌趋势明显。巴颜喀拉块体和四川盆地地壳存在 20km 的厚度差异, 使得龙门山构造带隆升幅度由 14.79% 增至 19.85%。综合分析 3 组离散元模拟实验结果, 得出巴颜喀拉地块下地壳物质层与四川盆地下地块物质层的黏度差异是龙门山构造带垂向隆升变形最关键控制因素的结论, 在下地壳黏度结构存在明显差异的前提下, 巴颜喀拉块体和四川盆地的地壳厚度差异对龙门山构造带纵向上逆冲隆升幅度有明显的促进作用。

龙门山构造带; 离散元数值模拟; 下地壳异质性; 地壳厚度; 主控因素

龙门山逆冲构造带是青藏高原东部巴颜喀拉地块与四川盆地的地质、地貌和地球物理边界, 涉及印支造山期中国大陆主体拼合以及喜马拉雅造山期青藏高原隆升两大构造事件[1‒3]。在经历破坏力极大的 2008 年汶川地震后, 作为最活跃的构造边界和青藏高原周缘生长的关键区域, 龙门山逆冲构造带成为全球构造地质学家研究逆冲构造隆升机制的热点地区。关于龙门山逆冲构造带隆升变形机制的讨论集中在地壳缩短(crustal shortening)[4‒5]和下地壳管道流(poiseuille flow)[6‒8]两大机制。地壳缩短机制认为, 板块碰撞挤压应力引起上地壳在水平方向上发生缩短, 由此产生的断层及其相关褶皱构造和逆冲推覆构造是龙门山地区地壳加厚隆升和地形高差巨大的主要原因。下地壳管道流机制认为, 青藏高原深部存在黏度低的易发生韧性流变变形的物质层, 其东南向流动受到四川盆地阻挡, 在龙门山地区形成陡峭的地形。前者将研究重点集中在块体间的相互作用方面, 缺少对块体内部运动和应力状态的分析; 后者则是在理想状态下建立流变学模型, 忽略地壳缩短量以及青藏高原隆升所产生的高地形重力载荷的影响。

前人通过低温热年代学、大地电磁探测和地震动力学数值模拟等手段, 得出地壳缩短[9]、地壳黏度结构[10‒12]和地壳厚度差异[12‒13]是晚新生代龙门山构造隆升变形重要影响因素的结论。本文基于前人对该地区的地质和地球物理研究成果, 利用离散元模拟技术, 对颗粒的运动状态进行记录和定量分析, 获取更多、更准确的动力学信息, 通过单因素变量控制方法, 探讨挤压背景下的下地壳黏度差异和地壳厚度对龙门山构造带隆升变形的影响, 分析龙门山构造带的隆升变形机理。

1 区域地质背景

龙门山构造带位于特提斯‒喜马拉雅构造域与滨西太平洋构造域的交接转换部位[14‒15], 东侧以安县‒灌县断裂与四川盆地分界, 西侧以茂县‒汶川断裂与巴颜喀拉地块分界, 南部为青藏高原东南缘的川滇构造带, 北部为秦岭造山带南缘的米仓山构造带[16](图 1(a))。龙门山构造带是典型的陆内挤压造山带, 晚三叠世之前处于被动大陆边缘拉张背景下, 晚三叠世末期转变为构造挤压环境。受印支期造山运动以及喜马拉雅造山运动的影响, 巴颜喀拉块体变形剧烈, 龙门山构造带发生强烈的冲断隆升[17]。龙门山褶皱冲断带平均海拔高度为 2000～4000km (图 1(b)), 其主峰九顶山海拔高度达到 4989km。沿NE-SW 走向, 龙门山构造带北段地形梯度小于中段和南段[18]。在研究区范围内, 龙门山构造断裂带位于莫霍面深度变化最为强烈的地带, 自西北向东南方向, 巴颜喀拉地块和四川盆地之间存在巨大的莫霍面深度差(图 1(c))。

龙门山构造带的总体长度约为 600km, 自西向东发育 3 条显露于地表的主干冲断裂: 青川‒茂汶断裂、北川‒映秀断裂和安县‒灌县断裂[19‒20]。青川‒茂汶断裂为龙门山构造带与松潘‒甘孜褶皱带的边界断裂[21]; 北川‒映秀断裂为汶川大地震的发震断裂, 其地表破裂带长度为 240～300km[22]; 安县‒灌县断裂为其南东边界。晚新生代以来, 它们都有明显的活动迹象, 且具有发生强震的可能性。它们在垂向上呈叠瓦状展布, 向四川盆地方向逆冲推覆, 在地下 20km 深处汇聚成一条大型剪切带[23‒24]。

(a)研究区地形及断裂分布; (b)龙门山断裂带地形剖面; (c)龙门山断裂带莫霍面深度。地形数据来源于地理空间数据云(http://www.gscloud. cn), 地震数据来源于中国地震台网(https://news.ceic.ac.cn), 莫霍面深度数据来源于Crust1.0 (https://igppweb.ucsd.edu/～gabi/crust1.html)

2 离散元数值模拟

2.1 离散元方法及实验参数

Cundall 等[25]于 1979 年提出离散元(discrete ele-ment merhod)方法的理论基础。该数值模拟方法基于颗粒单元间的接触准则, 利用时间步长和有限差分方法求解系统内每个颗粒单元的牛顿运动学方程[25‒26]。求解时, 将地质目标体离散为质量有限的圆盘状的刚性颗粒, 颗粒间为柔性接触, 接触范围极小(点接触), 接触位置可发生重叠, 且存在特殊的链接强度。通过调整颗粒强度、颗粒之间的黏结强度以及内摩擦系数等微观参数, 设定和模拟不同的岩石类型和岩性组合。基于动态应力松弛显式差分法求解运动方程和动力方程, 采用搜索算法搜索颗粒体接触条件, 并计算接触受力状态, 将颗粒单元间的不平衡力重新作用到节点上, 迭代至整个系统不平衡力足够小或者节点位移趋于平衡为止, 达到有效模拟自然重力沉积环境下非连续介质大应变量构造变形的目的[27‒28]。近年来, 离散元数值模拟方法广泛应用于构造地质研究, 如伸展变形[29]、板块俯冲变形[30]和褶皱冲断带构造变形[31]等。本文利用 PFC (particle flow code) 2D 软件进行离散元数值模拟。

2.1.1几何参数

随着观测数据的增加和地球物理学方法(接收函数、面波层析成像、体波层析成像和深地震测深等)的进步, 研究区地壳结构模型愈发精细。研究表明, 龙门山地处莫霍面向西倾斜的陡变带, 其西侧巴颜喀拉块体的上地壳、中地壳和下地壳厚度分别为 20～22, 26～30 和 14～18km, 地壳总厚度达到60～70km; 其东侧四川盆地的上地壳、中地壳和下地壳厚度分别为 10～12, 10～15 和 15～20km, 地壳总厚度仅为 30～40km[32‒34]。依据研究区实际地壳结构, 本文离散元数值模型与实际地质体的尺寸按 1:1000 的比例设置, 长 400m, 宽 100m。为模拟莫霍面的西倾形态, 将模型底部设置为非水平边界, 在轴的 100～275m 区间为倾斜边界, 倾角大约为6°。墙体法向刚度和切向刚度均设置为 1010N/m。在 400m´100m 区域内生成随机粒子, 颗粒半径在0.8～1.0m 之间, 不同粒径的颗粒数量服从高斯分布。加载重力加速度为 9.81m/s2, 允许颗粒聚集, 自然压实。达到平衡状态时, 删除顶部多余颗粒, 仅保留 400m´60m范围内颗粒。

2.1.2微观物理学参数

本文研究区下地壳实际黏度系数高达 1019～1021Pa·s[35‒37], 如果在离散元模型中如此设定, 单个模型的模拟运行时间会长达数年。本文参照前人的模型设定方法以及宏观参数与微观参数关系测定结果[38‒39], 并进行大量参数标定, 在满足模拟需求和保证模拟有效性的前提下, 为了提高计算效率, 将二维离散模型中下地壳物质层颗粒的黏度系数缩小1010倍, 设定为 109～1011Pa·s。

合理的参数是保证数值模拟实验结果有效性的前提。在通过二维颗粒流分析程序(PFC2D)进行离散元数值模拟实验时, 多通过双轴力学模拟实验的方法来确定颗粒的微观参数[40‒41]。从前人的研究结果中获得的研究区岩石力学参数如下: 巴颜喀拉块体以四川盆地上地壳物质层的内聚力分别为 30 和50MPa[42], 内摩擦角分别为 10°和 30°[43], 杨氏模量均为 8.75×1010Pa[44], 泊松比均为 0.25[45]。巴颜喀拉地块以及四川盆地中下地壳物质层的杨氏模量均为1.0×1011Pa[44], 泊松比均为 0.25[45]。以上述岩石力学参数为基础, 进行双轴力学数值模拟实验调试, 施加轴向恒定压缩速率, 平台伺服施加指定围压, 利用 Servo 系统和 Get_gain 函数获得岩石试样破坏全应力‒应变曲线。通过围压‒最大轴主应力对应的应力摩尔圆及包络线, 对比实际的材料参数, 确定模型微观输入参数, 获得各层 PFC2D 微观参数。本文采用双轴数值试验设定实验参数, 并参考文献[30,40,46‒ 49]。模型各层微观物理学参数见表 1。

2.1.3边界条件

地表 GPS 速度观测数据显示, 印度板块以 40mm/a 的速度沿喜马拉雅山脉向欧亚板块推进, 青藏高原受到强烈挤压发生连续变形, 中部的巴颜喀拉块体整体呈东南向运动, 速率约为 21mm/a[50‒51]。依据动力学相似原则设定实验模型变形速率, 则Ramberg (m)系数(无量纲)的计算公式[52]如下:

式中,,,,和分别表示密度、重力加速度、几何长度、黏度以及变形速率, 下角标m和n分别代表实验模型和自然界地质原型。在理想条件下, 地质原型和实验模型的m值应相同或属于同一量级。将本文地质原型及表 2 中参数代入式(1), 得到实验模型左侧的挤压速率为 6×10−3mm/s。

2.2 实验设计

2.2.1黏度结构对地表变形的影响

在龙门山构造带深部, 黏度呈现强烈的横向变化, 东西两侧下地壳存在 1～2 个量级的黏度差异。巴颜喀拉块体下地壳黏度范围为 1.48´1017～1.12´1021Pa·s, 平均值为 1.0´1019Pa·s; 四川盆地下地壳黏度范围为 4.09×1019～7.08×1020Pa·s, 平均值为 1.0×1020Pa·s[35‒37]。为考察下地壳黏度结构对龙门山构造带隆升变形的影响, 本文设计两组对比模型: 模型 1和模型 2。

模型 1 将 400m´60m 范围内的随机粒子划分为 4 个颗粒集, 自上而下分别为上地壳、中地壳和下地壳, 各层颗粒沉积密度均沿轴正方向减小。下地壳分为左右两段, 左段软弱下地壳的颗粒黏度系数设为 109Pa·s, 右段坚硬下地壳的颗粒黏度系数设为 1010Pa·s, 分别对应巴颜喀拉地块下地壳和四川盆地下地壳(图 2(a))。上地壳和中地壳内设标志层, 便于后期识别构造变形。为了保证力和力矩的有效传递, 下地壳层颗粒接触类型为线性模型, 上地壳层和中地壳层颗粒接触类型为线性平行黏结模型。模型右侧固定, 左侧以 6´10−3mm/s的速度挤压。记录挤压推进 2.5%, 5%和 7.5%时的变形情况(图 3)。

表1 离散元模拟微观物理学参数

表2 实验材料参数及模型相似比

说明: 模型各层密度与自然界地质体实际密度一致; 黏度数值参考文献[35,53‒54]。

模型 2 仅在纵向上做颗粒集的分层划分, 将400m´60m 范围内的随机粒子分为 3 个颗粒集(图2(b))。下地壳物质层黏度无差别, 为 1010Pa·s。其他设置(地壳厚度、推进速率和颗粒接触类型)均与模型 1 相同。记录挤压推进 2.5%, 5%和 7.5%时的变形情况(图 4)。

2.2.2地壳厚度差异对地表变形的影响

在印度洋板块的北向挤压下, 青藏高原具有明显的 NNE 向运动趋势, 受西伯利亚板块和塔里木地块的阻挡, 发生快速隆起, 其中部的巴颜喀拉块体成为巨厚陆壳体[55‒56]。

为探讨巴颜喀拉块体与四川盆地地壳厚度差异对龙门山构造带隆升变形的影响, 本文设计模型 3, 与模型 1 进行对比。

模型 3 保留 400m´40m范围内的颗粒, 分为 4 个颗粒集。自上而下划分为上地壳、中地壳和下地壳 3 层。下地壳划分为软弱下地壳和坚硬下地壳, 黏度系数均与模型 1 一致, 各物质层颗粒沉积密度沿轴正方向保持不变。模型 3中地壳厚度均一, 即莫霍面不存在起伏, 因此模型3 的底板设置为水平边界(图 2(c))。其他设置(黏度系数、推进速率和颗粒接触类型)均与模型 1 相同。记录挤压推进 2.5%, 5.0%和 7.5%时的变形情况(图 5)。

3 离散元数值模拟结果

提取二维离散元模拟变形结果(图 3～5)及颗粒运动矢量图(图 6), 后者记录了模型颗粒在末位置处的累积速度。模型 2 和 3 在近挤压端出现速度场奇点异常现象(图 6(b)和(c)中红色箭头所示), 奇点数量占比低于 0.05%, 对模型整体变形趋势及速度场的影响可忽略不计。

如图 3 所示, 模型 1 中, 随着左端的持续挤压, 当应力超过一定范围时, 颗粒在互斥作用下发生剪切错动和分离, 宏观上表现为地层垂向的逆冲抬升和水平方向的缩短变形。在挤压量达到 2.5%的早期挤压变形阶段, 左段黏度小的软弱下地壳受左侧挤压作用发生横向流动, 带动其上部物质层同步右移, 受到右段黏度大的不易流动的坚硬下地壳阻挡, 出现堆积聚集现象, 在中部(下地壳物质层属性变化处)出现褶曲, 高角度逆冲断裂 F1和 F2相继形成。F2与随后形成的反向调节断层 F3构成“V”字形的断层组合, 模型左部形成正断层 F4。在挤压量达到5.0%的中期挤压变形阶段, 模型左段断层数量增多, 靠近挤压端形成冲断断层 F7, 模型左段和右段差异变形态势逐渐增强。挤压量达到 7.5%时, 断层数量未增加, F1, F2和F6的断距小幅度增加。

图2 离散元模拟实验的初始模型设计

图3 离散元模拟实验模型1变形过程

图4 离散元模拟实验模型2变形过程

图5 离散元模拟实验模型3变形过程

图6 离散元模拟实验颗粒运动矢量对比

如图 4 所示, 模型 2 中右段发育褶曲叠瓦状逆冲断裂, 而在物质层厚度相对较大的左段和中部未发生明显的变形。模型颗粒运动速率整体上较小, 为 0～0.005m/s。沿轴正方向, 模型颗粒的向上速度分量逐渐增大, 远挤压端颗粒显示较强的上涌趋势。

如图 5 所示, 模型 3 中的变形集中在黏度系数较小的软弱下地壳区段, 近挤压端出现反冲断层, 模型中部(下地壳物质层属性变化处)形成变形强烈的隆起构造。软弱下地壳层颗粒具有较大的水平速度分量, 相对运动速率为 1.0～1.75m/s, 其上覆中地壳和上地壳物质层颗粒具有较大的垂向速度分量, 相对运动速率为 1.25～2.25m/s。坚硬下地壳上覆中地壳和上地壳物质层颗粒运动速率趋近零。

4 讨论

模型 1 左右两段下地壳的黏度和厚度存在明显差异, 左侧推进挤压作用导致模型中部(龙门山构造带)发生隆升变形, 纵向影响范围为 94.74%, 隆升幅度为 19.85%。以下地壳物质层属性变化处为界, 其左右两侧的模型颗粒运动方向和相对运动速率存在较大的差异, 右侧颗粒仅在有限空间范围内做低速运动, 左侧颗粒的相对运动速率整体上较大, 且在水平方向上具有明显的分段性。软弱下地壳物质层颗粒具有较大的水平运动速度分量, 上地壳和中地壳物质层具有较大的纵向运动速度分量(图 3和图 6(a))。

上述结果与前人通过分析该地区的重力异常数据、三维成像反演结果及 GPS 水平位移速度等资料得出的认识匹配度良好。在欧亚板块框架下, 巴颜喀拉地块的构造运动特征较为复杂, 主要处于NW-SE向拉张和NE-SW向挤压变形状态, 块体内部的运动速率及其差异性明显, 自西向东运动速率逐渐减小[57]。中国地壳运动观测网络数据[58]显示, 巴颜喀拉块体中部运动速率约为 20mm/a, 东部约为 10mm/a。巴颜喀拉地块是青藏高原深部软弱物质发生东南向“逃逸”及块体相互作用、变形的主要承载体[57‒69]。四川盆地深部存在高速度、低电导率的下地壳物质层, 对青藏高原软弱下地壳的东南向“逃逸”起拦截作用[4,60]。

模型 2 中地壳厚度存在明显差异, 左右两段下地壳的黏度系数相同。在模型左侧的挤压下, 变形出现在物质层厚度较薄的右端, 全空间区域颗粒相对运动速率均小于 1m/s。综合分析模型 1 和 2 的变形结果(图 3 和 4)和颗粒运动速度场(图 5)可知, 巴颜喀拉地块下地壳物质层与四川盆地下地壳物质层黏度结构的差异是龙门山构造带纵向隆升变形的先决条件。在下地壳黏度结构没有明显差异的前提下, 来自巴颜喀拉块体的南东向挤压应力以及龙门山构造带东西两侧地壳厚度的差异, 对龙门山构造带纵向地貌形态的塑造作用不明显。

模型 3 中, 左右两段下地壳的黏度存在明显的差异, 地壳厚度相同。在左侧的挤压作用下, 近挤压端发育反冲断层, 模型中部(龙门山构造带)发育构造三角带, 纵向影响范围为 93.17%, 隆升幅度为14.79%。

综合分析模型 1 和模型 3 的变形结果(图 3 和 5)及颗粒运动速度场(图 6)可知, 巴颜喀拉块体与四川盆地地壳厚度的差异对龙门山构造带垂向隆升变形起到一定的促进作用。在下地壳黏度结构存在明显差异的前提下, 龙门山东西两侧地壳厚度相同时, 隆升幅度仅为 14.79%; 龙门山东西两侧地壳存在20km 的厚度差时, 隆升幅度可以达到 19.85%。与模型 3 相比, 模型 1 左侧地壳、中地壳以及下地壳物质层颗粒相对运动速率在水平方向具有明显的分段性, 纵向上具有较高的一致性。由此推测, 地壳厚度差异产生的压力梯度能够促进软弱下地壳的水平运动, 进而影响上、中地壳的运动状态。在下地壳物质层属性变化处, 模型 1 和 3 的颗粒相对运动速率及变形趋势均急剧降低, 说明龙门山东西两侧下地壳黏度的差异是其构造带纵向隆升变形的先决条件。

5 结论

本文基于 PFC2D 离散元数值模拟对比实验, 针对地壳黏度结构和地壳厚度, 研究龙门山构造带隆升变形的影响作用, 得到以下主要结论。

1)下地壳异质性是龙门山构造带纵向隆升变形的最为重要的控制因素。在巴颜喀拉地块与四川盆地下地壳物质层黏度系数没有明显差异的情况下, 地壳厚度的差异不会形成龙门山的纵向高地貌形态。

2)在下地壳黏度结构存在明显差异的前提下, 巴颜喀拉块体与四川盆地的地壳厚度的差异对龙门山构造带的纵向逆冲隆升幅度有明显的促进作用。龙门山东西两侧地壳厚度存在 20km 的厚度差时, 隆升幅度可以达到 19.85%; 龙门山东西两侧地壳厚度相同时, 隆升幅度仅为 14.79%。

[1]Clark M K, Bush J W M, Royden L H. Dynamic topography produced by lower crustal flow. against rheological strength heterogeneities bordering the Ti-betan Plateau. Geophysical Journal International, 2005, 162(2): 575‒590

[2]Harrowfield M J, Wilson C J L. Indosinian deforema-tion of the Songpan-Garze Fold Belt. northeast Tibe- tan Plateau. Journal of Structural Geology, 2005, 27: 101‒117

[3]Wilson C J L, Harrowfield M J, Reid A J. Brittle modification of Triassic architecture in eastern. Tibet: implications for the construction of the Cenozoic plateau. Journal of Asian Earth Sciences, 2006, 27(3): 341‒357

[4]Hubbard J, Shaw J H. Uplift of the Longmen Shan and Tibetan plateau, and the 2008 Wenchuan (=7.9) earthquake. Translated World Ssmology, 2009, 458: 194‒197

[5]Tapponnier P, Xu Z Q, Roger F, et al. Oblique step-wise rise and growth of the Tibet Plateau. Science, 2001, 294: 1671‒1677

[6]Clark M K, Royden L H. Topographic ooze: building the eastern margin of Tibet by lower crustal flow. Geology, 2000, 28(8): 703‒706

[7]Cook K L, Royden L H. The role of crustal strength variations in shaping collisional orogens, with app-lication to Tibet. Journal of Geophysical Research: Solid Earth, 2008, 113(B8): B08407

[8]颜丹平, 孙铭, 巩凌霄, 等. 青藏高原东缘龙门山前陆逆冲带复合结构与生长. 地质力学学报, 2020, 26(5): 615‒633

[9]谭锡斌, 徐锡伟, 李元希, 等. 龙门山中段中央断裂和前山断裂的晚新生代垂向活动性差异及其构造意义. 地球物理学报, 2015, 58(1): 143‒152

[10]朱涛, 詹艳, 孙翔宇, 等. 四川龙门山断裂带高精度地壳/岩石圈黏度结构及其动力学意义. 地球物理学报, 2020, 63(1): 196‒209

[11]Liu C, Zhu B J, Yang X L, et al. Crustal rheology control on earthquake activity across the eastern mar-gin of the Tibetan Plateau: insights from numerical modelling. Journal of Asian Earth Sciences, 2015, 100: 20‒30

[12]England P, Houseman G. Finite strain calculations of continental defornation. 2. comparison with the India-Asia collision zone. Journal of Geophysical Research: Solid Earth and Planets, 1986, 91(B3): 3664‒3676

[13]李廷栋. 青藏高原隆升的过程和机制. 地球学报, 1995(1): 1‒9

[14]Enkin R J, Yang Z Y, Chen Y, et al. Paleomagnetic constraints on the geodynamic: history of the Major Blocks of China from the permian to the present. Journal of Geophysical Research: Solid Earth, 1992, 97(B10): 13953‒13989

[15]Yin A, Harrison T M. Geologic evolution of the Himalayan-Tibetan Orogen. Annual Review of Earth and Planetary Sciences, 2000, 28: 211‒280

[16]李智武, 宋天慧, 王自剑, 等. 川西‒龙门山盆山系统走向差异演化的变形、隆升和沉积记录及关键构造变革期讨论. 成都理工大学学报(自然科学版), 2021, 48(3): 257‒282

[17]Jia D, Wei G, Chen Z, et al. Longmen Shan fold-thrust belt and its relation to the western Sichuan Basin in central China: new insights from hydrocar-bon exploration. AAPG Bulletin, 2006, 90(9): 1425‒ 1447

[18]邓宾, 何宇, 黄家强, 等. 龙门山褶皱冲断带扩展生长过程——基于低温热年代学模型证据. 地质学报, 2019, 93(7): 1588‒1600

[19]金文正, 汤良杰, 杨克明, 等. 川西龙门山褶皱 冲断带分带性变形特征. 地质学报, 2007, 81(8): 1072‒1080

[20]Jin W, Tang L, Yang K, et al. Segmentation of the Longmen Mountains thrust belt, Western Sichuan Foreland Basin, SW China. Tectonophysics, 2010, 485: 107‒121

[21]郑勇, 李海兵, 王焕, 等. 印支期龙门山断裂带的逆冲‒推覆构造和沉积响应. 地质论评, 2018, 64(1): 45‒61

[22]李海兵, 王宗秀, 付小方, 等. 2008 年 5 月 12 日汶川地震(Ms 8.0)地表破裂带的分布特征. 中国地质, 2008, 35(5): 803‒813

[23]何祥丽, 李海兵, 王焕, 等. 蠕滑断裂带岩石组成和构造特征分析: 以龙门山灌县‒安县断裂带为例. 岩石学报, 2020, 36(10): 299‒314

[24]Burchfiel B C, Royden L H, Van der Hilst R D, et al. A geological and geophysical context for the Wen-chuan earthquake of 12 May 2008, Sichuan, People’s Republic of China. GSA Today, 2008, 18(7): 4‒11

[25]Cundall P A, Strack O D L. A discrete numerical mode for granular assemblies. Géotechnique, 1979, 29(1): 47‒65

[26]Dean S L, Morgan J K, Fournier T. Geometries of frontal fold and thrust belts: insights from discrete element simulations. Journal of Structural Geology, 2013, 53(8): 43‒53

[27]Morgan J K, McGovern P J. Discrete element simu-lations of gravitational volcanic deformation: 1. de-formation structures and geometries. Journal of Geo-physical Research: Solid Earth, 2005, 110(B5): 1671‒ 1677

[28]Naylor M, Sinclair H D, Willett S, et al. A discrete element model for orogenesis and accretionary wed-ge growth. Journal of Geophysical Research: Solid Earth, 2005, 110(B12): 1‒16

[29]周易, 于福生, 刘志娜. 分层伸展叠加变形二维离散元模拟. 大地构造与成矿学, 2019, 43(2): 213‒225

[30]王一丹, 于福生, 刘志娜, 等. 板块俯冲变形过程二维离散元模拟——对东海陆架盆地成因启示. 海洋地质与第四纪地质, 2019, 39(5): 166‒176

[31]李长圣. 基于离散元的褶皱冲断带构造变形定量分析与模拟[D]. 南京: 南京大学, 2019

[32]He R, Shang X, Yu C, et al. A unified map of Moho depth andpsratio of continental China by receiver function analysis. Geophysical Journal International, 2014, 199(3): 1910‒1918

[33]Wang W, Wu J, Fang L, et al. Crustal thickness and Poisson’s ratio in southwest China based on data from dense seismic arrays. Journal of Geophysical Research: Solid Earth, 2017, 122(9): 7219‒7235

[34]高天扬, 丁志峰, 王兴臣, 等. 利用接收函数、面波频散和ZH振幅比联合反演青藏高原东南缘地壳结构及其动力学意义. 地球物理学报, 2021, 64(6): 1885‒1906

[35]朱涛, 詹艳, 孙翔宇, 等. 四川龙门山断裂带高精度地壳/岩石圈黏度结构及其动力学意义. 地球物理学报, 2020, 63(1): 196‒209

[36]孙玉军, 董树文, 范桃园, 等. 中国大陆及邻区岩石圈三维流变结构. 地球物理学报, 2013, 56(9): 2936‒2946

[37]Diao F, Wang R, Wang Y, et al. Fault behaviorand lower crustal rheology inferred from the first seven years of postseismic GPS data after the 2008 Wen-chuan earthquake. Earth and Planetary Science Letters, 2018, 495: 202‒212

[38]Egholm D L, Sandiford M, Clausen O R, et al. A new strategy for discrete element numerical models: 2. sandbox applications. Journal of Geophysical Re-search: Solid Earth, 2007, 112: B05204

[39]Liu Z, Koyi H A. The impact of a weak horizon on kinematics and internal deformation of a failure mass using discrete element method. Tectonophysics, 2013, 586: 95‒111

[40]李维波, 李江海, 王洪浩, 等. 库车前陆冲断带克拉苏构造带变形影响因素分析——基于离散元数值模拟研究. 大地构造与成矿学, 2017, 41(6): 1001‒ 1010

[41]刘畅, 陈晓雪, 张文, 等. PFC数值模拟中平行粘结细观参数标定过程研究. 价值工程, 2017, 36(26): 204‒207

[42]Chéry J, Zoback M D, Hasani R. An integrated me-chanical model of the San Andreas fault in central and northern California. Journal of Geophysical Research: Solid Earth, 2001, 106(B10): 22051‒22066

[43]He J K, Lu S J, Wang W M. Thre-dimensional mecha-nical modeling of the GPS velocity field around the northeastern Tibetan plateau and surounding regions. Tectonophysics, 2013, 584: 257‒266

[44]Zhu S B, Zhang P Z. FEM simulation of interseismic and coseismic deformation asociated with the 2008 Wenchuan Earthquake. Tectonophysics, 2013, 584: 64‒80

[45]尹力, 罗纲. 有限元数值模拟龙门山断裂带地震循环的地壳变形演化. 地球物理学报, 2018, 61(4): 1238‒1257

[46]Finch E, Hardy S, Gawthorpe R. Discrete element modelling of contractional fault-propagation folding above rigid basement fault blocks. Basin Research, 2004, 25(4): 515‒528

[47]Ai J, Langston P A, Yu H S. Discrete element model-ling of material non-coaxiality in simple shear flows. International Journal for Numerical and Analytical Methods in Geomechanics, 2014, 38: 615‒635

[48]Zhang Y, Teng J, Wang Q, et al. Density structure and isostatic state of the crust in the Longmenshan and adjacent areas. Tectonophysics, 2014, 619/620(6): 51‒57

[49]Zhu S, Zhang P. FEM simulation of interseismic and coseismic deformation associated with the 2008 Wen-chuan Earthquake. Tectonophysics, 2013, 584: 64‒80

[50]Liang S, Gan W, Shen C, et al. Three-dimensional velocity field of present-day crustal motion of the Tibetan Plateau derived from GPS measurements. Journal of Geophysical Research: Solid Earth, 2013, 118(10): 5722‒5732

[51]张培震, 邓起东, 张国民, 等. 中国大陆的强震活动与活动地块. 中国科学: D辑, 2003, 33(4): 12‒20

[52]Ramberg H. Gravity, deformation and the Earth’s crust in theory, experiments and geological applica-tion. 2nd edition. Kolkata: Academic Press, 1981

[53]Diao F, Wang R, Wang Y, et al. Fault behavior and lower crustal rheology inferred from the first seven years of postseismic GPS data after the 2008 Wen-chuan earthquake. Earth and Planetary Science Let-ters, 2018, 495: 202‒212

[54]Huang M H, Bürgmann R, Freed A M. Probing the lithospheric rheology across the eastern margin of the Tibetan Plateau. Earth and Planetary Science Letters, 2014, 396: 88‒96

[55]Hao M, Li Y, Zhuang W. Crustal movement and strain distribution in East Asia revealed by GPS observations. Scientific Reports, 2019, 9(1): 16797

[56]瞿伟, 高源, 陈海禄, 等. 利用 GPS 高精度监测数据开展青藏高原现今地壳运动与形变特征研究进展. 地球科学与环境学报, 2021, 43(1): 182‒204

[57]李承涛, 苏小宁, 孟国杰. 巴颜喀拉块体东北缘GPS应变率空间分布特征及其与 2017 年九寨沟Ms 7.0地震的关系. 地震, 2018, 38(2): 37‒50

[58]徐锡伟, 闻学泽, 陈桂华, 等. 巴颜喀拉地块东部龙日坝断裂带的发现及其大地构造意义. 中国科学:地球科学, 2008, 38(5): 529‒543

[59]陈长云, 任金卫, 孟国杰, 等. 巴颜喀拉块体东部活动块体的划分、形变特征及构造意义. 地球物理学报, 2013, 56(12): 4125‒4141

[60]Wang Z, Wang X, Huang R, et al. Structural hetero-geneities in Southeast Tibet: implications for regio-nal flow in the lower crust and upper mantle. Inter-national Journal of Geophysics, 2012, 315: 1‒12

Main Controlling Factors of Uplift Deformation of Longmenshan Structural Belt: Insight from Discrete Element Method

WANG Ying1,2, LI Jianghai1,2,†, MA Changming1,2, SONG Juechen1,2

1. Key Laboratory Orogenic Belts and Crustal Evolution (MOE), School of Earth and Space Sciences, Peking University, Beijing 100871; 2. Institute of Oil and Gas, Peking University, Beijing 100871; † Corresponding author, E-mail: jhli@pku.edu.cn

In order to explore the main controlling factors of uplift deformation of Longmenshan structural belt, based on the differences in the properties of lower crust material layer between the east and west sides of Longmenshan structural belt, three groups of PFC2D discrete element numerical simulation are carried out to realize quantitative analysis.The experimental deformation results and the model particle motion vector map show that under the condition of no obvious difference in the material properties of the lower crust, the existence of plate collision and compression stress and crustal thickness difference will not form a huge topographic elevation difference in the Longmenshan structural belt. When there are obvious differences in the viscosity coefficient of the lower crust, the relative value of the particle movement rate of the weak lower crust material layer is 1.5‒2.94 m/s, and the average movement rate is 1.62 m/s, which is about 54 times of the average movement rate of the particles of the hard lower crust layer. Uplift deformation occurs in the middle of the model (Longmenshan structural belt), with a vertical influence range of 94.74% and a uplift amplitude of 19.85%. The particles of the middle crust and upper crust overlying the weak lower crust have a large upward velocity component, and the upward trend of the material layer of the upper crust is obvious. There is a 20km thickness difference between Bayankala block and the crust of Sichuan Basin, which increases the uplift amplitude of Longmenshan structural belt from 14.79% to 19.85%. Based on the comprehensive analysis of three discrete element simulation experiments, it is concluded that the viscosity difference between the material layer of the lower crust of Bayan Kara block and the material layer of the underground block of Sichuan Basin is the most key control factor for the vertical uplift deformation of Longmenshan structural belt. On the premise that there are obvious differences in the viscosity structure of the lower crust, the crustal thickness differences between the Bayan Kara block and the Sichuan Basin significantly promote the vertical thrust uplift amplitude of the Longmenshan structural belt.

Longmenshan structural belt; discrete element numerical simulation; heterogeneity of the lower crust; crustal thickness; main controlling factors

10.13209/j.0479-8023.2022.056

2021–07–15;

2021–10–15

国家科技重大专项(2016ZX05033002, 2016ZX05033001)资助