库车前陆冲断带西部古近系盐构造三维离散元数值模拟

2020-04-01李江海章雨王洪浩王殿举

石油勘探与开发 2020年1期

李江海，章雨，王洪浩，王殿举

（1.造山带与地壳演化教育部重点实验室北京大学地球与空间科学学院，北京 100871；2.北京大学石油与天然气研究中心，北京 100871）

0 引言

塔里木盆地库车前陆冲断带油气资源丰富，尤其是克深—大北含油气区区带面积达629.4 km2，天然气资源量达10 797.62×108m3[1]，具有良好的勘探前景。该区含丰富盐岩层，盐岩流动变形强烈，盐上层多发育与滑脱相关的褶皱，盐下层变形以冲断为主，发育断层相关褶皱及基底断块[2-5]。前人相关研究主要在以下几个方面：①依据薄膜理论，借助弹塑性模型，认为库车前陆冲断带西段盐构造发育的初期阶段，控制盐体运动与盐构造发育的主要因素是由拜城凹陷巨厚沉积物所导致的沉积差异负载作用，而后期影响库车坳陷盐构造发育的最主要因素则是挤压应力[6]；②基于地震解释和平衡剖面研究，提出了吐孜玛扎盐墙和却勒盐推覆体的演化模式[7]；③基于物理模拟实验认为库车盐构造变形主要受区域挤压和同构造沉积作用的影响[8]；④综合二维物理模拟实验和数值模拟研究，总结了盐岩边界、区域构造应力变化、差异负载（沉积负载和局部构造负载）3个因素对库车坳陷盐构造的影响作用[9]；⑤基于二维离散元模拟研究，认为应变速率大小和作用时间、盐岩展布形态、先存盐底辟、基底先存断裂以及基底古隆起等因素控制了库车前陆冲断带盐构造的形成[10]。然而，以上研究均是基于二维数据的研究结果，仅能对单一剖面进行模拟，不能从三维空间反映出盐岩在垂直于主应力方向上的侧向流动及构造样式，造成了盐构造形成机制认识不透彻。

相较于二维离散元模型颗粒间的赫兹模型（Hertz-Mindlin model），三维离散元数值模拟提供的线弹性模型（linear springs model），能够对颗粒和墙面单元的法向刚度和切向刚度两个参数进行设定，地下盐岩近似于牛顿流体，当颗粒与颗粒或者颗粒与墙体之间发生接触时，它们之间的作用力可通过力-位移关系计算得到，在分析大位移问题上具有更大的优越性[11]。同时，允许颗粒之间进行平行粘结，类似于有一定体积的胶结材料作用在距接触点一定范围内，既能传递力又能传递力矩，可以更好的表现不同岩石的性质，更适用研究颗粒集合体的破裂以及流动问题[12]。本文结合前人理论认识，以克深和大北两个工区为例，采用三维离散元数值模拟方法，通过对盐构造主应力方向和垂直挤压方向的构造分析，探讨库车前陆冲断带西部盐构造形成的控制因素及其形成机理。

1 地质背景

库车前陆冲断带位于塔里木盆地北部，夹持于南天山造山带与塔北隆起之间，古近纪至中新世，库车地区处于盐湖、干旱湖泊-盐湖的环境下，主要沉积了古近系库姆格列木群膏盐岩和新近系吉迪克组膏盐岩，构成了地区最主要的滑脱层，以库车为界，前陆冲断带西部主要发育库姆格列木群（E1—2km）盐岩，分布范围较广，厚度较大，最厚处超过6 000 m；东部主要发育吉迪克组（N1j）盐岩，分布面积较小，厚度也相对较小，最厚处约1 800 m[4-5,13-16]（见图1）。强烈的挤压作用和盐岩滑脱层的存在导致库车前陆冲断带构造变形在整体上具有“南北分带、东西分段、上下分层”的特征，即从西到东以乌什凹陷、拜城凹陷和阳霞凹陷为主体分为 3段；自北而南划分为北部单斜带、克拉苏—依奇克里克构造带、乌什—拜城—阳霞凹陷带、秋里塔格构造带、南部斜坡带共 5个构造变形带[3,17-19]，其中，收缩型强变形带（克拉苏—依奇克里克构造带、秋里塔格构造带）分隔弱变形带（北部单斜带、中部凹陷带、南部斜坡带）。克深—大北三维地震工区位于库车前陆冲断带西部，库姆格列木群盐岩的存在，导致本区具有分层变形的特征[6]。

图1 库车前陆冲断带构造分区及盐岩层厚度分布图

2 克深—大北三维工区地震解释

基于三维地震数据解释可以看出，克深和大北工区地下盐岩的空间展布存在较大差异（见图2）。在东部的克深1-2工区，盐岩主要集中在北部由盐下逆冲断层形成的断层三角带中，南部则发育盐撤向斜，盐岩在向斜核部最薄，在部分区域甚至出现盐焊接构造；而在西部的大北工区，盐岩主要集中在南部的盐背斜（大宛齐背斜）中，在背斜核部盐岩最厚可超过6 000 m。在大北工区北部，盐岩受吐孜玛扎断层影响，发育盐刺穿构造。在部分区域，盐岩随逆冲断层冲破地表，形成盐墙（吐孜玛扎盐墙）构造。在位于克深 1-2工区和大北工区之间的克深 5工区，地下盐岩空间展布具有过渡的性质，在北部和南部均发育有盐岩厚度中心。

图2 克深—大北三维地震工区联合深度剖面

以过克深 4井的骨干剖面（深度剖面）为例（见图3），克深工区变形主要集中在北部的挤压端，盐上层、盐岩层及盐下层变形样式存在较大区别。盐上层变形以宽缓的大规模向斜为主，并显示了生长地层的特征，在靠近挤压端的部分区域，发育有滑脱褶皱和冲断层；盐岩层厚度从北向南变化较大，盐岩多集中在靠近北部的断层三角带中，而南部盐岩层厚度较薄；盐下层发育叠瓦状逆冲断层，越靠近挤压端，盐下地层的抬升幅度越大，断层的断距也越大，远离挤压端的部分区域，发育有双冲构造。另外，从东到西，克深工区在平行主应力方向剖面上的变形样式大致相似。

以过大北 3井的骨干剖面（深度剖面）为例（见图4），大北工区盐岩主要集中在南部的盐背斜中，在剖面北部盐岩发生底辟形成盐墙（吐孜玛扎盐墙）构造，在盐背斜与盐墙接触的部位，盐岩减薄。盐下层变形样式与克深工区类似，发育叠瓦状逆冲断层。盐上层在远离挤压端的位置发育背斜构造，在盐背斜和盐墙接触的部位发育小规模向斜构造，向斜的一翼被盐墙截切。从东到西，大北工区在平行主应力方向剖面上的变形样式存在一定区别，主要表现在工区中部平行主应力方向的剖面上南部大宛齐背斜的隆升幅度较大，盐岩层厚度较大，而向两侧大宛齐背斜隆升幅度逐渐减小，盐岩层厚度也逐渐减小。

图3 库车前陆冲断带西部克深工区过克深4井骨干地震剖面解释图（剖面位置见图1，据文献[10]）

图4 库车前陆冲断带西部大北工区过大北3井骨干地震剖面解释图（剖面位置见图1）

3 离散元数值模拟方法和初始模型设计

3.1 离散元方法和实验参数

离散元方法是一种基于离散颗粒间接触准则的数值模拟方法。通过时间-位移有限差分方法，计算颗粒在牛顿定律下移动，能有效模拟弹塑性变形过程[20]。离散元方法基于应力与位移准则进行微粒位移变形计算，并允许微粒间的滑动与破裂，适应于大应变量构造变形样式分析模拟实验，被广泛应用于模拟浅层地壳结构、断裂系及剪切带变形过程[21-23]。本文采用离散元三维模拟的线弹性模型，其颗粒和墙面单元具有法向刚度和切向刚度两个基本参数，当颗粒与颗粒或者颗粒与墙体之间发生接触时，它们之间的作用力可通过力-位移关系计算得到。

数值模拟研究中，首先设定模型不同的边界条件，经过多次调试，建立基础几何模型。模型与实际地质体尺寸按1︰100比例设置，在长250 m（x轴）、宽250 m（y轴）、高150 m（z轴）的空间范围内，按照颗粒半径为1.5～2.0 m随机生成颗粒，后设置重力因素，让颗粒在重力状态下自然压实并达到平衡，然后逐步删除顶部的颗粒，仅保留250 m×250 m×60 m范围内的颗粒。后对颗粒进行属性赋值，球心z轴坐标范围为0～20 m的颗粒设置为盐下层，范围在20～35 m的颗粒设置为盐岩层，范围在35～60 m的颗粒设置为盐上层。选取每层的中间层位设置标志层，便于对后期变形构造进行识别。

合理的参数是数值模拟实验的基础。在通过三维颗粒流分析程序（PFC3D）进行离散元数值模拟实验时，多通过三轴力学模拟实验的方法来确定颗粒的微观参数。三轴力学模拟实验模型是由一个圆柱形墙面围绕生成的颗粒体以及上下两面压缩板组成的。其中，上下两面压缩板在试样压缩过程中，上墙面对试样进行压缩剪切，下墙面保持不动；圆柱形墙面在压缩剪切过程中通过伺服系统维持恒定围压的工作。应力应变数据的收集与整理是根据跟踪墙面所受力的情况以及墙面相对位移决定的。通过大量的三轴力学模拟实验，并不断调整模型中颗粒的参数，根据比较模拟结果与试验室内的结果，两者吻合度很高时，数模模型中的各参数值就是实际材料的各微观参数值。但在具体操作过程中，由于缺少库车前陆冲断带西部盐上层、盐岩层及盐下层完整的三轴力学实验数据，基于塔里木物探中心地质力学室岩石试验及前人研究可获得的库车前陆冲断带西部盐上层、盐岩层和盐下层部分岩石力学参数如下：盐上层和盐下层岩石密度分别为2 400 kg/m3和2 600 kg/m3，岩石的弹性模量、泊松比和内聚力均为60 MPa、0.25和10 MPa，内摩擦角分别为35°和40°。盐岩层黏度为1×1019Pa，体积模量和剪切模量分别为1×104MPa和1×103MPa[9,24-26]。本文在实验参数的选取中采用三轴力学模拟实验调试并参考前人相关模拟参数的方法[27-29]，具体采用的颗粒相关参数可见表1。

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

本文认为地下盐岩近似于牛顿流体，采用的模型为线性模型，而盐上地层和盐下地层则均采用常用的平行粘结模型[30-32]，颗粒之间通过平行粘结进行接触，既能传递力，也能传递力矩。在盐岩颗粒与盐上地层及盐下地层接触的边界，颗粒之间采用的模型也是平行黏结模型，具体参数可见表1。此外，在模拟过程中墙体的法向刚度和切向刚度均设置为 1×1010N/m，墙体与颗粒之间的摩擦系数设置为0.5，模型右侧刚性墙体以1 m/s的速率向x轴负方向进行挤压。

3.2 初始模型设计

在以上工作的基础上，实际工区的初始模型建立，参考了前人物理、数值模拟实验以及平衡剖面研究总结的演化模式[7,33]及盐岩变形的影响因素，包括：盐岩分布[6]、基底构造[15]、先存断裂、先存底辟[10]以及南天山差异推覆作用影响[34-35]。

3.2.1 克深工区初始模型

模型在挤压端设置了 2条先存断裂，分别为过点（170，0，20）与y轴平行，与xOy平面夹角呈 45°的平面以及过点（205，0，25）与y轴平行，与xOy平面夹角呈60°的平面。先存断裂导致靠近挤压端的基底发生抬升，抬升呈阶梯状，被抬升的盐下层厚度分别为25 m、30 m，而远离挤压端盐下层的厚度为20 m。基底形态的改变同样影响了盐岩层的厚度，在靠近挤压端盐岩层厚度最薄，为5 m，而远离挤压端盐岩层的厚度为15 m，盐岩层顶部与水平面平行（见图5）。

此外，还在上述模型的基础上添加了颗粒的同沉积作用。本模型中同沉积颗粒共设置了 2批，在初始模型中设置了第1批同沉积颗粒，第1批同沉积颗粒共5 000个，大小为1.5～2.0 m，颗粒在x轴范围为0～125 m、y轴范围为0～250 m、z轴范围为100～125 m内随机生成（见图5），颗粒的自身属性与盐上层颗粒属性保持一致，在重力作用下，逐渐下落至稳定，用以模拟同沉积作用在挤压过程中对于模型的影响。模型在运行100 000步时，设置第2批同沉积颗粒，第2批同沉积颗粒的生成范围、属性均与第 1批同沉积颗粒相同。

图5 库车前陆冲断带西部克深工区三维离散元模型图

3.2.2 大北工区初始模型

大北工区南部大宛齐背斜处发育古盐湖的中心[36]，因此本模型在远离挤压端的位置增加了盐岩层的厚度，且盐岩层的厚度向周缘呈阶梯状递减，用来模拟古盐湖沉积。其中在x轴范围为30～110 m、y轴范围为85～165 m，盐岩厚度最厚，为27 m；在x轴范围为20～120 m、y轴范围为75～175 m，盐岩层厚度为23 m；在x轴范围为10～130 m、y轴范围为65～185 m，盐岩层厚度为19 m；在其余范围内，盐岩层厚度为15 m（见图6）。此外，大北工区北部吐孜玛扎盐墙由先存盐底辟发育而来[7]，因此本模型在靠近挤压端的位置设置了先存线状盐底辟。先存线状盐底辟位于盐岩层之上，其延伸方向与y轴方向平行，底辟最右端距挤压墙体65 m，宽度为20 m，高度为20 m，底辟顶部在xOz平面上具有一定的弧度。

然后在模型靠近挤压端的位置设置了 5条先存断裂，分别为过点（125，0，20）、（150，0，20）、（175，0，20）、（200，0，20）、（225，0，20）与y轴平行，与xOy平面夹角呈45°的平面。本模型与克深工区模型不同，未设置靠近挤压端的基底隆起，先存断裂未刺穿到盐岩层内部。此外，本模型同样添加了颗粒的同沉积作用，同沉积颗粒共设置了 2批，在初始模型中设置了第1批同沉积颗粒，第1批同沉积颗粒共2 000个，颗粒大小为1.5～2.0 m，颗粒在x轴范围为100～150 m、y轴范围为0～250 m、z轴范围为100～120 m随机生成（见图6）。同样，模型在运行100 000步时，设置第2批同沉积颗粒，第2批同沉积颗粒的生成范围、属性均与第1批同沉积颗粒相同。

图6 库车前陆冲断带西部大北工区三维离散元模型图

4 离散元数值模拟变形结果分析

4.1 克深工区模型模拟结果

从主应力方向的剖面（过点（0，125，0），法向量为（0，1，0））中可以看出，模型的变形主要集中在挤压端，且变形是前展式的，靠近挤压端的盐下层先存断裂最先活化，基底早期隆起发生进一步抬升，随着挤压作用的持续进行，模型的缩短量进一步增加，先存断层倾角逐渐增大，断距逐渐增加，刺穿到盐岩层内部。且向远离挤压端的方向，盐下地层新的逆冲断层逐渐发育，形成叠瓦状逆冲构造（见图7）。靠近挤压端的盐上地层在挤压作用下形成大规模背斜，而同沉积作用导致远离挤压端的盐上地层厚度增加，并形成大规模的向斜构造。随着挤压作用的持续进行，向斜和背斜的结合部位逐渐变得陡立，在模型运算到200 000步时，其翼部发生倒转，形成倒转背斜和倒转向斜构造（见图7）。

考虑到后期的剥蚀去顶作用，对模型中部平行主应力方向的剖面进行顶部截切（见图8b）。基于对比可以看出，模型运算140 000步后过模型中部平行主应力方向剥蚀去顶后的剖面与克深工区过克深 4井骨干地震剖面（见图3）在形态上具有较大的相似性。

对比远离挤压端（过点（100，0，0），法向量（1，0，0））和靠近挤压端（过点（150，0，0），法向量（1，0，0））垂直主应力方向的剖面，模拟结果可以看出，在远离挤压端垂直主应力方向的剖面上，各层变形较弱，盐内幕标志层基本未发生变形，而在靠近挤压端垂直主应力方向的剖面上，由于剖面横切过盐岩的厚度中心，盐岩层厚度较大，且盐岩层顶面呈现出不规则的样式（见图9）。此外，在靠近挤压端垂直主应力方向的剖面上，盐内幕标志层可见复杂的褶皱样式，部分区域存在标志层错断的现象，表明靠近挤压端的位置，盐岩存在较强的侧向流动。

图10计算了模型140 000步颗粒运动矢量。从图中可以看出，随着挤压作用的持续进行，模型中颗粒靠近挤压端运动速率较大，远离挤压端运动速率非常小，运动速率最大的区域为盐上层靠近挤压端的位置，盐岩层至盐下层颗粒运动速率逐渐减小，说明构造变形基本集中在靠近挤压端，从盐下层至盐上层变形逐渐增强，本区岩石容易破裂形成逆冲断裂。此外，盐岩层颗粒具有相对较大的向上的速度分量，导致挤压区岩层逐渐上拱形成背斜，且垂直主应力方向上颗粒有向中部运动的速度分量，且越靠近中部、运动速率越大，说明背斜内部盐岩层较厚并具有复杂的内幕构造。

图7 克深工区典型盐构造演化剖面图（主应力方向）

图8 克深工区典型盐构造演化剖面图（模型中部位置）

图9 克深工区典型盐构造演化剖面图

图10 克深工区典型盐构造演化颗粒运动矢量图

4.2 大北工区模型模拟结果

从大北工区模拟结果主应力方向剖面中可以看出（见图11），模型的变形是前展式的，变形更多的集中在挤压端。盐下层变形相对简单，先存断裂发生再活化，断层的断距逐渐增加，且靠近挤压端的断层倾角逐渐增大。在模型运算150 000步及200 000步过模型中部平行主应力方向的剖面上，可见新的基底断裂发育。盐岩层的变形以流动变形为主，持续的挤压作用导致靠近挤压端的先存盐底辟在剖面上的宽度逐渐变窄，但底辟高度逐渐增加，形成盐墙构造，且盐墙逐渐向远离挤压端发生倾斜。在远离挤压端，前期设置的盐岩厚度中心在挤压作用下形成盐背斜。而同沉积作用形成的差异负载导致模型中部的盐岩层发生减薄，盐岩更多的向远离挤压端的盐背斜聚集。盐上层变形以褶皱变形为主，在挤压作用下，靠近挤压端和远离挤压端的盐上层发育背斜，而模型中部在同沉积作用下发生沉降，形成向斜构造。同时，盐底辟两侧盐上层的标志层可见错断现象。

基于模型中部平行主应力方向剥蚀去顶后的剖面与大北工区过大北 3井骨干地震剖面，在形态上具有较大的相似性（见图12）。除盐下层逆冲断层的具体数量及盐墙的形态与实际地震剖面存在一定差别外，模拟结果中盐上层、盐岩层及盐下层的变形样式与实际地震剖面非常类似。

图11 大北工区典型盐构造演化剖面图（主应力方向，地层划分及挤压方向同图6）

图12 大北工区典型盐构造演化剖面图（模型中部位置）

图13 大北工区典型盐构造演化剖面图

远离挤压端垂直主应力方向选取过点（70，0，0）、法向量（1，0，0）的剖面，剖面横切过盐背斜；靠近挤压端垂直主应力方向选取过点（155，0，0）、法向量（1，0，0）的剖面，剖面横切过盐墙。从图13中可以看出，在远离挤压端垂直主应力方向的剖面上，各层变形较弱，盐内幕标志层仅在局部形成小规模褶皱，表明在远离挤压端盐岩的侧向流动较弱。而在靠近挤压端垂直主应力方向的剖面上，由于剖面横切过盐底辟，剖面上盐岩展布较为复杂，部分盐体在剖面上与盐源层分离。且由于边界墙体的阻碍作用，在剖面的两侧，可见盐上层标志层被盐底辟错断的现象。盐内幕标志层相比远离挤压端剖面变形更加复杂，褶皱变形的幅度更加强烈，在部分区域可见标志层错断的现象，表明在靠近挤压端盐岩的流动性更强。

图14为模型运算140 000步颗粒运动矢量。从图中可以看出，随着挤压作用的持续进行，模型中颗粒靠近挤压端运动速率较大，远离挤压端运动速率较小，但速率减慢相对平缓，远离挤压端仍具有一定的运动速率，说明在先存线状盐底辟及初始盐岩沉积中心的影响下，地层受挤压变形的范围更大。先存底辟、盐岩沉积中心及其盐上层颗粒均有相对较大的向上的速度分量，导致靠近挤压端底辟逐渐上拱形成盐墙、远离挤压端盐岩及上覆层逐渐上隆。此外，垂直主应力方向上两侧颗粒有向中部运动的速度分量，中部颗粒则表现为向上和沿挤压力方向向前的运动趋势，导致盐岩逐渐向中部聚集并在中部上拱，盐内幕构造更加复杂。

图14 大北工区典型盐构造演化颗粒运动矢量图

5 讨论

三维离散元数值模拟过程较为复杂，计算量与二维离散元数值模拟相比呈指数增长，对特定研究区的定量模拟难度较大[20-23]。库车前陆冲断带西部克深、大北工区盐构造样式复杂且研究区范围较大，在模拟中主要采用定性模拟的方法，通过设定模型不同的边界条件，来探讨本区盐构造形成的主控因素。与二维模拟结果相比，三维模拟可以更好地探讨不同因素对于盐构造及盐内幕构造的影响，能够对主应力方向剖面及垂直主应力方向的剖面模拟结果进行综合分析，两个方向上盐构造均有明显的分段性，且模拟结果与实际地质现象较为符合。

研究区是由于印度板块和欧亚板块碰撞，导致天山新生代再活化，并在南天山山前形成前陆冲断带[37]。克深工区和大北工区处于相同的盆山构造体系，盐构造在主应力方向上变形整体上是前展式的，主要集中在山前，局部存在垂向迁移。克深工区山前的挤压应力导致基底先存断裂发生再活化，使早期隆升的基底进一步发生抬升[15]，随着挤压应力的向前传播，新的逆冲断层在盐下不断发育，形成叠瓦状逆冲构造。盐岩层的形态受到山前挤压作用和同沉积作用的双重影响，远离山前的同沉积作用会对盐岩层形成差异负载，导致远离山前的盐岩具有盐撤的趋势，而山前的挤压作用导致盐岩层有向远离山前流动的趋势，在同沉积作用和挤压作用的双重作用下，远离山前的盐岩层厚度基本保持稳定，而在山前，盐岩同时具有沿断裂隆升部位向上运动的速度分量，并在盐下基底断层形成的盐三角带中聚集，形成本区盐岩的厚度中心。大北工区来自山前的挤压应力同样导致基底先存断裂发生再活化，逐渐发育新的逆冲断层并形成叠瓦状逆冲构造。此外，渐新统苏维依组—第四系砾岩为山前河流-洪积扇相沉积，前积作用导致盐上存在差异负载，进而形成先存底辟[7]。挤压应力导致先存线状盐底辟进一步发育，盐岩在沿主应力方向运动时相对于周围地层有较大的垂向运动的速度分量，形成吐孜玛扎盐墙。据文献[36]，大宛齐背斜处发育有古盐湖的沉积中心，原地积聚有较厚的盐层，盐岩后期同沉积作用形成的差异负载导致工区中部的盐岩层发生减薄，盐岩继续在大宛齐盐背斜处聚集并上隆，形成工区的厚度中心。

垂直主应力方向的剖面显示，在挤压端位置，两个工区两侧盐岩有向中部运动的速度分量，中部则表现为向上和沿挤压力方向向前的运动趋势，且越靠近中部运动速率越大，导致盐岩逐渐向中部聚集并在中部上拱，盐内幕构造更加复杂，形成的盐构造以复杂褶皱样式为主，且存在与源盐分离的构造组合，按此趋势，在部分区域（尤其在工区中部），盐岩会随逆冲断层冲破地表；而在远离挤压端，盐岩仍有向中部运动的速度分量，但侧向流动能力相对降低，各层变形减弱，在非厚层原始盐岩沉积区，盐内幕标志层仅在局部形成小规模褶皱。