博格达山位于新疆北部，整体呈东西向展布，将南部的吐哈盆地和北部的准噶尔盆地分隔开来，在区域沉积物源供给和构造演化中占据重要位置（王国灿等，2020）（图 1a）.博格达山北缘冲断带也是整个东北天山褶皱冲断带的前锋（Zhou et al.，2020），在新生代印亚碰撞的远程效应影响下变形强烈（Tapponnier and Molnar，1979；Avouac et al.，1993；Yin et al.，1998），是研究东北天山新生代变形扩展的有利区域.博格达山北缘先后发现了甘河、三台等油气田以及古牧地背斜等多个含油气构造（康竹林，1997），是有利的油气勘探区域（文志刚等，2005；邹会明，2017），新生代构造活动对其沉积层序和油气运聚有重大影响（孙自明和沈杰，2014）.因此，探索博格达山北缘新生代的变形机制不仅对研究北天山新生代扩展方式有重要意义，还可以为该区域油气勘探提供依据.

图  1  博格达山西段地质构造简图

F1.三工河断裂; F2.妖魔山断裂: F3.阜康断裂; A1.苦坝沟背斜; A2.七道湾背斜; A3.古牧地背斜；a.研究区的大地构造位置；b博格达山西段构造纲要简图，修改自Chen et al.（2015）和马德龙等（2017）

Fig.  1.  Simple geological map of west Bogda Mountain

下载: 全尺寸图片 幻灯片

众多学者基于野外地质调查、地质和地震资料解析、生长地层和热年代学分析等方法探究博格达山北缘的构造特征及构造演化过程.研究结果显示，该地区在时间上具有多旋回构造叠加的特点，发育有明显的先存构造（Tang et al.，2014；Chen et al.，2015）；在空间上呈现自山前基底卷入构造向盆地区薄皮逆冲构造转换的特征（吴建华等，2002；陈科等，2012；孙自明和王毅，2014）.前人的物理模拟和数值模拟实验结果表明，存在先存构造的位置会优先发生构造活动（Sassi et al.，1993；刘恒麟等, 2022）；滑脱层的分布特征、内聚力强度和厚度都对褶皱冲断带变形的样式和传播范围有控制作用（吴航等，2019；辛文等，2020；Li et al.，2021）；存在双滑脱层的情况下，上部滑脱层有利于构造变形的解耦（Buiter，2012；Graveleau et al.，2012）.博格达山北缘中侏罗统西山窑组滑脱层和下侏罗统八道湾组滑脱层的滑脱性能存在明显差异，可能在研究区构造演化中起到重要作用（商琳等，2013；马德龙等，2017）.然而双滑脱层的内聚力强度及其差异在褶皱冲断带变形中的作用尚未可知，先存构造和双滑脱层在研究区新生代构造变形中的具体作用仍待进一步研究.

本次研究针对博格达山北缘冲断带西段米泉地区（图 1b），采用离散元数值模拟方法，结合构造解析和应力应变分析，通过控制变量法设计5组数值模拟实验来探究先存构造和双滑脱层内聚力强度对博格达山北缘新生代构造变形的影响，为该地区构造变形特征研究和油气勘探提供依据.

1.   地质背景

博格达山是天山的东北支，平均海拔4 000 m以上，东西向延伸超过250 km，总体上为向北凸出的弧形复背斜（陈科等，2012）.作为中亚造山带的一部分，博格达山记录了3次主要的构造事件：古生代与古亚洲洋闭合有关的碰撞和拼合、中生代欧亚大陆陆内造山以及新生代印亚碰撞的远程效应（Tapponnier and Molnar，1979；Şengör et al.，1993；Jahn et al.，2004）.关于博格达山的构造演化还存在一些争议，但可以确定的是该地区至少在早侏罗世就已经呈现正地貌特征（Greene et al.，2005；Chen et al.，2015；张妍等，2015），并在中生代持续变形（沈传波等，2006），现今地形是晚中新世以来在古造山带基础上强烈变形发育的继承性构造（王宗秀和李涛，2004；Zhou et al.，2020）.博格达山广泛出露中-上石炭统至新近系，该地区综合地层柱状图如图 2所示.中-上石炭统为一套海相中-基性侵入岩和火山岩，主要出露在博格达山复背斜核部（Tang et al.，2014）.二叠纪沉积环境由石炭纪的浅海体系演变为滨岸体系，中-下二叠统主要由含有腕足类和珊瑚化石的长石砂岩、粉砂岩、凝灰岩和石灰岩透镜体组成，上二叠统主要为非海相碎屑沉积物（Wartes et al.，2002）.三叠系呈现向上逐渐变细的沉积旋回，下三叠统主要由红色砾岩和一些泥质砂岩组成，中-上三叠统为灰绿色泥岩、泥质粉砂岩和砂质泥岩夹砾岩，整体反映了水深增加和湖泊扩张的环境（李忠和彭守涛，2013）.侏罗系以砂、泥岩沉积为主，层序保存完整，可划分为下侏罗统八道湾组和三工河组，中侏罗统西山窑组和头屯河组，以及上侏罗统齐古组和喀拉扎组，其中八道湾组和西山窑组是研究区主要的构造滑脱层（商琳等，2013）.下部八道湾组（J₁b）厚约670 m，夹薄煤层，上部西山窑组（J₂x）厚约600 m，富含煤层（郭召杰等，2011；陈科等，2012），岩性特征表明两套滑脱层的内聚力存在显著差异.不整合覆盖在侏罗系之上的白垩系只在研究区北侧古牧地地区出露，主要由灰绿色泥岩、砂岩组成.上覆新生界广泛分布在博格达山的北麓，由陆相泥岩、砂岩和砾岩组成.中-上二叠统、上三叠统、和中下侏罗统是该地区主要的烃源岩层系（文志刚等，2005）.

图  2  博格达山北缘综合地层柱状图

修改自马德龙等（2017）；Zhou et al.(2017)

Fig.  2.  Comprehensive stratigraphic column of northern Bogda Mountain

下载: 全尺寸图片 幻灯片

研究区位于山前褶皱冲断带，构造变形较为复杂，地震反射资料品质差，导致前人对研究区深层结构的认识存在较大的差异.对冲断带山前部分的争议主要集中在断裂是否切穿基底，吴建华等（2002）认为七道湾背斜及其深部构造属于基底卷入构造；孙自明和王毅（2014）认为博格达山山前是一系列沿二叠系滑脱层发育的逆冲推覆构造.冲断带盆地部分发育深浅两套冲断系统，浅部沿中侏罗统西山窑组滑脱层发育突破型断层传播褶皱（陈莹莹等，2022），而深部构造样式存在不同认识：（1）深部发育双重堆垛构造，具有统一的中侏罗统西山窑组顶板，底部具有基底卷入成分（吴建华等，2002）；（2）深部发育构造三角带（孙自明和王毅，2014）；（3）深部构造是以中侏罗统西山窑组和二叠系上芨芨槽群芦草沟组为顶、底板的双重构造（商琳等，2013）；（4）根据贯穿古牧地背斜的长山1井的油气来源判断双重构造未向下贯穿至二叠系，深部构造为以侏罗系西山窑组和八道湾组为顶、底板的双重构造（邹会明，2017）.本文综合前人的研究进展，根据米泉1井、长山1井和牧6井的钻测井数据、地表露头以及油气来源资料，在对研究区最新的二维和三维地震资料解释的基础上，建立了横穿整个博格达山北缘冲断带的地质剖面AA'（图 3）.

图  3  研究区典型剖面图

剖面位置见图 1

Fig.  3.  Typical geological profile of the study area

下载: 全尺寸图片 幻灯片

剖面自造山带向盆地方向呈现基底卷入构造向薄皮滑脱构造转换的特征，在博格达山山前发育三工河、妖魔山、米泉3条高角度基底断裂.野外观测到的不整合特征及热年代学工作都揭示博格达山及其北缘经历了多期构造隆升（Wartes et al.，2002；沈传波等，2006；朱文斌等，2006；汪新伟等，2007；陈科等，2012；Tang et al.，2014），因此山前的3条高角度基底断裂都经历了多期活动，对新生代构造变形来说属于先存断裂.盆地区冲断带发生上下变形解耦现象，浅部沿中侏罗统西山窑组煤系地层滑脱发育断层传播褶皱，断面由南向北逐渐变陡直至突破地表形成古牧地断裂，根据断裂上盘的生长地层判断，断裂的主要活动时间在中新统之后（商琳等，2013）；深部是以中侏罗统西山窑组和下侏罗统八道湾组煤系地层为顶、底板滑脱层的双重构造，马德龙等（2017）认为深部双重构造的变形时间主要为晚侏罗世-早白垩世和晚白垩世，这和米泉地区最新三维地震剖面揭示的双重构造之上发育白垩系和侏罗系、古近系和白垩系两个不整合面是吻合的.

综上所述，侏罗系八道湾组和西山窑组煤系滑脱层内聚力差异显著，可能对研究区新生代变形有着重要影响；山前3条基底卷入断裂以及古牧地地区深部双重构造在新生代之前就开始发育，对于新生代变形来说属于先存构造，对研究区新生代构造变形同样有控制作用.以横穿整个博格达山北缘冲断带的地质剖面AA'为地质模型，本文在设置先存构造的基础上，通过改变不同内聚力强度双滑脱层组合，利用数值模拟实验探讨研究区新生代构造变形特征及控制因素.

2.   实验方法与模型设计

2.1   离散元方法原理

离散元（Discrete Element Method，简称DEM）方法的基本思想是将材料视为由若干个离散单元组成的初始弹性颗粒系统，通过给定合适的微观参数来模拟不同材料的运动学和力学行为（Cundall and Strack，1979），可以有效地模拟弹塑性形变，表征非连续介质破裂和大尺度变形过程，已经被国内外学者广泛应用于褶皱-冲断带构造变形过程及机制的研究中（Hardy and Finch，2005；Morgan，2015；辛文等，2020；Li et al.，2021）.

离散元数值模拟的计算主要是采用时间步进的有限差分法，首先检索颗粒之间的接触关系，再根据接触力学模型计算颗粒在某时刻t的受力情况，进而采用跳蛙法求解牛顿力学方程，更新颗粒在t+∆t时的位置，完成一个时间步的计算，之后不断循环上述步骤，直至计算结束.本文颗粒间接触力计算采用Hertz-Mindlin接触模型，使用离散元数值模拟软件VBOX完成实验（Li et al.，2018，2021；李长圣，2019），基于实验结果讨论先存断层和滑脱层对博格达山北缘冲断带构造变形的控制作用.

2.2   实验模型设计

参考剖面AA'地质单元尺度、剖面长度、沉积序列与地层厚度，设计了长50 km、高6.5 km的初始模型.根据构造解析结果，在距挤压端0 km、8 km、16 km处分别设置3条倾角45°的基底卷入型先存断层（F1、F2、F3），并根据断层活动性分别设置2.5 km、2.5 km、0.5 km高的基底抬升.在古牧地地区侏罗系内设置4条倾角30°的先存断层（F4、F5、F6、F7），用来模拟中生代先存构造.研究区存在上、下两套滑脱层，为探究两套滑脱层的内聚力强度对冲断带构造变形的影响，设计5组对比实验（图 4）.实验1是不设置滑脱层的空白对照实验，实验2设置两套强内聚力滑脱层，实验3将上部设置为强内聚力滑脱层，下部设置为弱内聚力滑脱层（简称“上强下弱”模型），实验4将上部设置为弱内聚力滑脱层，下部设置为强内聚力滑脱层（简称“上弱下强”模型），实验5设置两套弱内聚力滑脱层.为控制单一变量，实验中滑脱层的厚度均设置为600 m，下滑脱层之下能干层厚1 km，滑脱层之间能干层厚0.8 km、上滑脱层之上能干层厚3.5 km，造山带根部地层被基底卷入断层错断抬升，所以滑脱层从F3下盘才开始设置.实验模型中灰色代表地层，红色代表断层，黄色代表弱内聚力滑脱层，白色代表强内聚力滑脱层，需要指出的是，内聚力强度越大的滑脱层滑脱性越弱.

图  4  实验初始模型设计

图a~e分别代表实验1~5

Fig.  4.  Pictures showing initial model design

下载: 全尺寸图片 幻灯片

本文中实验参数的设置参照前人的双轴实验结果（Morgan，2015），通过调整颗粒的细观参数及颗粒间的粘结参数来构建不同的地质单元（表 1）.实验中能干层的参数设置对应的岩石粘聚力Co为19.0 MPa，内摩擦角ϕ为19.3°（李长圣，2019）.模型左侧刚性挡板以2 m/s的速度由南东向北西挤压，时间步长0.05 s，总挤压量为15 km.

表  1  离散元数值模拟实验参数

Table  Supplementary Table   Experimental parameters of discrete element numerical simulation

地质单元	颗粒的细观参数					颗粒间的粘结参数				色标
	颗粒半径 （m）	剪切模量 （109 Pa）	泊松比	摩擦系数	密度 （103 kg/m3）	杨氏模量 (108 Pa)	剪切模量 (108 Pa)	抗拉强度 (107 Pa)	剪切强度 (107 Pa)
能干层	60/80	2.9	0.2	0.3	2.5	2.0	2.0	2.0	4.0
弱滑脱层	60/80	2.9	0.2	0.0	2.2	/	/	/	/
强滑脱层	60/80	2.9	0.2	0.1	2.2	/	/	/	/
断层	60/80	2.9	0.2	0.0	2.5	/	/	/	/
基底	60/80	2.9	0.2	0.4	2.5	2.0	2.0	2.0	4.0

下载: 导出CSV 

| 显示表格

3.   实验结果

3.1   实验1：无滑脱层模型

实验1是不设置滑脱层的空白对照试验，模拟褶皱冲断带对先存构造的响应特征，记录运算0 step、50 000 steps、100 000 steps和150 000 steps，分别对应挤压缩短量0 km、5 km、10 km和15 km，模拟结果如图 5所示.当挤压量达5 km时，先存断层F1和F2率先活化使得断层上盘地层持续变形形成背斜a1和a2，同时发育反冲断层f8将F1错断（图 5b）.当挤压量达10 km时，先存断层F3活化并冲出地表，F3上盘形成多条次级反冲断层（图 5c）.当挤压量增加至15 km，变形仍集中于冲断带根部，先存断层F4至F7基本未活化，前陆方向也未发生明显变形（图 5d）.

图  5  实验1模拟结果与构造解释

Fig.  5.  Simulation results and structure explanation of Experiment 1

下载: 全尺寸图片 幻灯片

3.2   实验2：双强内聚力滑脱层模型

实验2在实验1的基础上设置两套强内聚力滑脱层，记录运算0 step、50 000 steps、100 000 steps和150 000 steps，分别对应挤压缩短量0 km、5 km、10 km和15 km，模拟结果如图 6所示.当挤压量达5 km时，先存断层F1和F2率先活化使得断层上盘地层持续变形形成背斜a1和a2，同时发育反冲断层f8将F1错断（图 6b）.当挤压量达10 km时，沿两套滑脱层分别发生变形解耦，下滑脱层控制地层整体的向前滑动，而上滑脱层控制上能干层的滑动和褶皱变形，最终先存断层F4向上扩展并冲出地表，同时在F4上盘形成反冲断层f9并构成箱型褶皱a3.先存断层F3也被卷入变形，但未向上扩展（图 6c）.当挤压量增加至15 km，变形主要沿上滑脱层向前扩展，并在前端形成反冲断层f10，先存断层F5、F6和F7轻微活化（图 6d）.

图  6  实验2模拟结果与构造解释

Fig.  6.  Simulation results and structure explanation of Experiment 2

下载: 全尺寸图片 幻灯片

3.3   实验3：“上强下弱”滑脱层模型

实验3在实验1的基础上设置强内聚力的上部滑脱层和弱内聚力的下部滑脱层，记录运算0 step、50 000 steps、100 000 steps和150 000 steps，分别对应挤压缩短量0 km、5 km、10 km和15 km，模拟结果如图 7所示.当挤压量达5 km时，先存断层F1和F2率先活化并在断层上盘分别发育反冲断层f8和f9，形成背斜a1和a2（图 7b）.当挤压量达10 km时，变形迅速沿下滑脱层传播并沿先存断层F4突破地表，新发育的反冲断层f10和F4之间形成背斜a3.先存断层F3在下滑脱层处被错断，未向上扩展（图 7c）.当挤压量继续增加至15 km，先存断层F5至F7依次活化，变形沿下滑脱层持续向前传递并从F6处突破地表，形成背斜a4（图 7d）.整个变形过程中上滑脱层未发挥作用.

图  7  实验3模拟结果与构造解释

Fig.  7.  Simulation results and structure explanation of Experiment 3

下载: 全尺寸图片 幻灯片

3.4   实验4：“上弱下强”滑脱层模型

实验4在实验1的基础上设置弱内聚力的上部滑脱层和强内聚力的下部滑脱层，记录运算0 step、50 000 steps、100 000 steps和150 000 steps，分别对应挤压缩短量0 km、5 km、10 km和15 km，模拟结果如图 8所示.当挤压量达5 km时，先存断层F1和F2率先活化使得断层上盘地层持续抬升，同时发育反冲断层f8将F1错断（图 8b）.当挤压量达10 km时，先存断层F3重新活化并突破地表，F3上盘形成反冲断层f9，前陆方向沿两套滑脱层分别发生变形解耦，下滑脱层控制地层整体的向前轻微滑动，而上滑脱层控制上能干层强烈滑动和褶皱变形，沿上滑脱层发育逆冲断层f10（图 8c）.当挤压量继续增加至15 km，上滑脱层继续发挥主要作用，上能干层沿上滑脱层滑动并从断层f12冲出地表，在其上盘形成反冲断层f13和背斜a4.先存断层F5至F7微弱活化（图 8d）.

图  8  实验4模拟结果与构造解释

Fig.  8.  Simulation results and structure explanation of Experiment 4

下载: 全尺寸图片 幻灯片

3.5   实验5：双弱内聚力滑脱层模型

实验5在实验1的基础上设置两套弱内聚力滑脱层，记录运算0 step、50 000 steps、100 000 steps和150 000 steps，分别对应挤压缩短量0 km、5 km、10 km和15 km，模拟结果如图 9所示.当挤压量达5 km时，先存断层F1和F2率先活化使得断层上盘地层持续抬升，同时发育反冲断层f8将F1错断（图 9b）.当挤压量达10 km时，沿两套滑脱层分别发生变形解耦，下滑脱层控制地层整体的向前滑动，而上滑脱层控制上能干层的滑动和褶皱变形，最终先存断层F4向上扩展并冲出地表，同时在F4上盘形成反冲断层f9并构成箱型褶皱a3，先存断裂F5也开始重新活化（图 9c）.当挤压量继续增加至15 km，先存断层F6开始活化，上能干层继续沿上滑脱层向前扩展，形成断层f10和背斜a4（图 9d）.

图  9  实验5模拟结果与构造解释

Fig.  9.  Simulation results and structure explanation of Experiment 5

下载: 全尺寸图片 幻灯片

4.   结果分析与讨论

实验1~5的模拟结果和体积应变对比如图 10所示，体积应变图中红色代表右行剪切，蓝色代表左行剪切，颜色越深表示变形越强烈.对比各组实验结果和实际剖面，探讨先存断层和双滑脱层对褶皱冲断带构造变形的影响及对博格达山北缘冲断带新生代演化过程的启示.

图  10  模拟结果及体积应变对比

a1、a2对应实验1；b1、b2对应实验2；c1、c2对应实验3；d1、d2对应实验4；e1、e2对应实验5

Fig.  10.  Comparison of simulation results and volume strain

下载: 全尺寸图片 幻灯片

4.1   先存断层的影响

5组数值模拟实验结果显示，在挤压过程中，应变主要是以先存断层活化的形式向前扩展，先存断层吸收了大部分挤压缩短量，而新生断层大部分是在先存断层上盘形成的反冲断层或者是先存断层的延伸和扩展，先存断层对冲断带的变形样式起到重要控制作用.实验1不存在滑脱层，先存断层F3很快向上扩展并突破地表，而在实验2~5中，存在滑脱层的情况下，应力优先沿滑脱层传播，F3的活动性明显减弱（图 10），先存断层的活动性在一定程度上受到与之相交的滑脱层的影响.

4.2   双滑脱层的影响

对比5组实验的变形传播范围发现，实验1挤压应变集中在冲断带根部，基本没有向前陆方向传播（图 10a2），而实验2~4应力均能通过滑脱层快速向前陆方向传播（图 10b2、10c2、10d2、10e2），滑脱层是应力的快捷传播通道.

实验2双强内聚力滑脱层模型，应力同时沿两套滑脱层传播，其中上滑脱层的应力传播速度明显更快，传播距离也更远，控制着上能干层的远距离滑移（图 10b1、10b2）.实验5双弱内聚力滑脱层模型，相对于实验2应力传播更远，滑脱层活动性更强，但其构造变形模式和实验2相似（图 10e1、10e2）.对比结果显示，当两套相同强度滑脱层存在的情况下，两套滑脱层都发挥作用，但应力优先沿上滑脱层传播.

实验3“上强下弱”滑脱层模型，下滑脱层成为变形扩展的唯一通道，上滑脱层基本没有起到拆离作用（图 10c1、10c2）.实验4“下强上弱”滑脱层模型，上滑脱层是主要的应变传播通道，而下滑脱层有利于先存断层F4的活化，起到次要作用.因此，当存在两个不同内聚力滑脱层的情况下，应力会优先沿弱内聚力滑脱层传递，当下滑脱层内聚力较弱时，上滑脱层将基本不发挥作用.

4.3   对博格达山北缘冲断带构造变形的启示

博格达山北缘存在下侏罗统八道湾组和中侏罗统西山窑组两套滑脱层（商琳等，2013），其中八道湾组含煤条带和煤层，西山窑组富含煤层和煤纹，符合“上弱下强”模型.将实际地质剖面与实验结果对比发现，只有实验4的结果与实际地质剖面最为吻合，图 10a1、10a2、10a3、10a4模拟结果中与实际地质剖面中的博格达冲断隆起、苦坝沟背斜、七道湾背斜、古牧地背斜相对应，先存断层F1和F2的强烈活化变形也对应新生代博格达山的剧烈隆升（Chen et al.，2015）.古牧地地区下部双重构造主要在中生代活动，而上部背斜是沿西山窑组煤系滑脱层长距离滑脱并在古牧地冲出地表形成的断层传播褶皱（汪新伟等，2007；马德龙等，2017），模拟实验4结果与该认识相吻合（图 10d1、10d2）.以上结论说明先存断裂和双滑脱层的内聚力强度差异确实对褶皱冲断带构造变形有重要影响.从实验结果分析可知，下侏罗统八道湾组强内聚力滑脱层有助于先存双重构造的复活，而中侏罗统西山窑组弱内聚力滑脱层控制着浅部古牧地背斜的形成，两套滑脱层共同控制盆地区冲断带深浅构造的变形解耦，并且“上弱下强”的双滑脱层模式是控制研究区新生代构造变形的关键因素.

4.4   实验模型的局限性

本文实验采用的离散元数值模拟方法能很好地模拟弹塑性形变，表征剪切变形特征，但实际褶皱冲断带构造变形过程受到更多复杂因素影响，如盆山边界形态、基底性质差异等.本文采用的二维模拟技术没有考虑冲断带沿走向的变化，实验设计简化了地层岩性和产状，也没有考虑同构造沉积和剥蚀的影响.此外，由于先存断层的产状和分布范围难以准确厘定，故实验中先存断层的倾角和位置取了估计值.虽然存在上述近似和简化，离散元数值模拟实验仍然优化了物理模拟存在的比例化和流变学的问题，可以体现出褶皱冲断带的几何学和运动学特征及变形机制.

5.   结论

通过离散元数值模拟的方法，本文模拟了先存构造和不同内聚力强度双滑脱层组合对褶皱冲断带构造变形的影响，结合博格达山北缘冲断带构造特征，讨论其新生代变形机制，主要结论如下：

（1）滑脱层是应力传播的有利通道，双滑脱层的内聚力强度是褶皱冲断带构造变形的重要控制因素.

（2）存在两套同等内聚力强度滑脱层的情况下，应力优先沿上滑脱层传播.

（3）当两套滑脱层内聚力强度相差较大时，变形优先沿弱内聚力滑脱层传播，尤其是当下滑脱层内聚力强度较低时，上滑脱层可能不发挥作用.

（4）通过与地质模型对比，认为先存构造和“上弱下强”的双滑脱层模式是控制博格达山北缘新生代构造变形的关键因素.