数值模拟盐体几何起伏导致的应力扰动

王明文，罗纲，孙云强，常成

（1.武汉大学测绘学院，武汉 430079；2.中国科学院大学地球与行星科学学院，北京 100049）

0 引言

在油气盆地中，盐体与油气的运移、聚集、储存、勘探和开发都密切相关[1-3]。全世界约80%的含油气盆地为盐盆地[4]，在这些盆地中都发现了大量石油和天然气。膏盐层一般并不是油气钻探的主要目的层段，但是由于盐层具有极低的渗透率[5]，是油气藏良好的盖层，因此巨厚盐层之下往往蕴藏着大型油气田，例如美国的墨西哥湾盆地、西非盐盆地和巴西深海盆地[6]。但是，在盐盆地中，相对于周围的沉积物而言，盐体的力学性质比较特殊，具体表现为力学强度低、易于蠕变及塑性变形[5,7-8]，因此，在力平衡作用下，盐体周围的应力会因为其黏弹性应力松弛而发生扰动，且孔隙流体压力也会发生变化。所以在盐膏层钻井，特别是在深层盐膏层钻井时会因地层应力、孔隙压力的变化而造成一些重大工程事故，如井眼失稳、卡钻和挤毁套管等[9]。

近年来，已有不少学者使用地质力学方法对盐盆地进行建模（演化模型和静态模型）并模拟盐体周围的应力扰动[10-11]，发现盐体的几何形状对应力扰动具有重要影响[12-16]；但是并没有给出区域背景应力场样式（挤压或拉张背景应力）、盐体的几何形状（突起或洼陷）与应力扰动的关系。本文构建了一个挤压构造背景应力背景下的静态盐构造模型，模拟计算盐盆地周围的应力变化，分析不同几何形状的盐体所导致的应力扰动模式，同时测试盐体起伏大小对应力扰动大小的影响，并以塔里木盆地库车坳陷克拉苏西剖面处的盐体为例进行了实际模拟计算。本研究可为挤压型盐盆地系统应力扰动研究及相关钻井设计提供科学依据。

1 研究区概况

库车坳陷位于塔里木盆地北部，处于塔北隆起与南天山造山带之间，西起温宿，东至库尔勒[17]，东西长约500 km，南北宽20～100 km[18]。库车坳陷在第三纪沉积了巨厚的膏盐岩[19]，受中—晚期喜马拉雅运动的强烈挤压与上覆地层载荷的差异压实作用，形成了多种类型的盐构造[20]，本区域的盐构造由 3个部分组成，从下往上分为盐下沉积层、膏盐层、盐上覆层[21]。

研究区位于库车坳陷克拉苏构造带西端，其南北向水平长度为76 km，纵向深度为15 km；拜城凹陷南北两侧发育盐底辟构造，北侧的吐孜玛扎盐丘北翼发育有断层。南秋里塔格背斜和北秋里塔格背斜是两个紧闭的箱状褶皱，它们位于盐丘之上，盐丘的厚度为5～7 km；拜城凹陷上覆层厚度为6～7 km（见图1）。

研究区盐构造为典型的挤压型盐构造[22-26]，根据各时代地层岩石磁组构、水压致裂、井径崩落等多参数实测和解释，库车坳陷现今地应力场的最大主压应力方向为近南北向[27]，从库车坳陷三叠系克拉玛依组砂岩中发育逆冲断层可以判断库车坳陷构造应力场的中间主应力方向为近东西向，最小主应力为垂向[28]。

2 模型设置

图1 塔里木盆地北缘库车坳陷主要地质构造及库车坳陷克拉苏构造带以西地震剖面[17]

在自然界中单层、双层或多层盐构造普遍存在，如江汉盆地江陵凹陷八岭山—花园多层系盐构造、滇西兰坪盆地盐构造和墨西哥湾盆地中北部的盐构造[29-31]。本文构建了双层盐构造模型，同时也模拟测试了单层的盐构造模型，发现其结果和双层盐构造模型关于盐体几何形态对应力扰动的结论一致。笔者使用商业软件建立二维平面应变有限元模型[32]（见图2），该模型假设盐体为马克斯韦线性黏弹性体[33]，沉积物为弹性固体，即无孔隙压力；假设盐体不渗透[13]，也没有孔隙压力。模型的材料参数如表1所示。

图2 有限元网格和边界条件

表1 有限元模型材料参数

本研究中所使用的模型模拟深度为20 km，模型水平宽度为 45 km；模型中有上下两层盐体，对应着 3种盐体的几何形态，分别为：左侧都是突起的盐体，中间为凹凸组合的盐体，右侧都是洼陷的盐体。上层盐体的平均深度约为6.5 km，下层盐体的平均深度约为15.0 km。模型的左右两侧边界和底部边界都是滑移边界：底部边界只有水平向位移，没有垂向位移；左右边界只有垂向位移，没有水平向位移。上边界为自由边界。为了更好地展示盐体附近的应力扰动，对盐体附近的网格进行了加密。

模型的静态力学平衡方程为：

其中，∂σij为应力张量（i,j=1, 2, 3），ρgi为体力项。本研究中以挤压应力为正值。

模型每一步的应变增加表示为：

其中，{dεe}和{dεv}分别表示为黏性和弹性的应变增量[15,37]。黏性和弹性的关系由马克斯韦尔体线性黏弹性关系确定，具体表述为：

其中，tσ是t时刻的应力增量，{dσ}、dt分别是应力张量增量和时间增量，Q为黏度相关的矩阵，D为弹性材料矩阵。

对于模型的背景应力场，本文将模型的垂直应力（σv）设置为最小主应力σ3，水平应力（σh）和法向应力（σoop）分别对应最大主应力σ1和中间主应力2σ。初始的垂向应力等于上覆岩层的重力随深度的积分，在挤压型构造背景应力场下，水平应力设置为垂向应力的1.3倍，法向应力设置为垂向应力的1.1倍，即：

3 模拟结果

3.1 参考模型

本文假设了垂向、水平、法向的初始应力随着深度的增加而增加，所以初始的von Mises（剪应力的三维表示，见公式（8））应力也是随着深度逐渐增加而增加的（见图3a）。由于黏弹性盐体不能承受偏应力[12]，盐体内的剪应力随着时间的发展逐渐松弛减小，并在盐体周围的沉积层产生应力扰动（见图3b）。模型最终的应力状态（见图3c）是初始状态的应力（见图3a）与由于盐体的黏弹性应力松弛产生的应力扰动（见图3b）之和。本研究中挤压应力为正值。

图3 参考模型不同状态下的von Mises应力大小变化

模型中，在盐体突起处，水平应力、垂向应力和法向应力都逐渐增加（见图4a—4c）。水平应力和法向应力都是在盐体突起的核部显著增加，然后随着远离突起核部，应力扰动开始减小（见图4a、图4c）。与之不同的是，垂向应力在盐体突起的核部有明显的增加（见图4b虚线框内），但是在距离突起核部不远的两侧，应力扰动突然从挤压应力变为拉张应力（见图4b椭圆框）。水平应力、法向应力和垂向应力之所以如此变化，是因为突起的盐体附近的沉积物在水平方向从两侧分别向突起的核部移动，而在垂向上则向上移动（见图5a），由此产生了水平方向和垂向的挤压应变（见图5b、图5c）。

图4 参考模型盐体周围沉积物中的应力扰动与主应力方向

图5 盐体黏弹性应力松弛在沉积物产生的位移变化和应变扰动

与盐体突起处的应力扰动相反，在盐体洼陷处，水平应力、垂向应力和法向应力都会减小（见图4a—4c）。水平应力和法向应力都是在盐体洼陷的核部减小明显，随着远离洼陷核部，应力扰动开始减小（见图4a、图4c）；与之不同的是，垂向应力在盐体洼陷的核部有微弱的拉张应力，在垂向上远离核部一段距离后变为较大的拉张应力，出现一个拉张应力区域，但是在距离此拉张应力区域不远的两侧，应力扰动突然从拉张应力变为挤压应力（见图4b）。盐体洼陷处应力如此变化的原因在于，盐体附近的沉积物从洼陷核部向两侧移动，垂向上往下移动，由此位移产生了水平方向和垂直方向的拉张应变（见图5b、图5c）。

剪应力的变化主要集中在盐体几何形状不规则（突起和洼陷）的位置（见图4d）。在盐体突起左侧出现挤压剪应力，右侧出现拉张剪应力；与之相反的是，盐体洼陷左侧出现拉张剪应力，右侧出现挤压剪应力；但是随着沉积物远离几何形状异常盐体，剪应力很快减小并趋近于零。剪应力的扰动（见图4d）可以由剪切变形的变化来解释（见图5d），例如，盐体突起左侧物质流动表现为逆时针方向（见图5a，左侧虚线框），右侧为顺时针方向（见图5a，右侧虚线框）。

通过将主应力与正应力进行比较发现（见图4），3个主应力与水平应力、法向应力和垂向应力的扰动模式上（挤压或者拉张）有很大的相似性，但大小有区别，表现为：水平应力对应最大主应力，垂向应力对应最小主应力，法向应力对应中间主应力。这是因为在模型设置时确定了模型初始状态的水平应力为最大主应力、法向应力为中间主应力，垂向应力为最小主应力。在盐体的应力松弛之后，剪应力变化最大的区域主要出现在盐体起伏较大的区域附近，因此主应力方向在盐体突起和洼陷附近有所变化，发生了旋转，而在水平盐层或盐体几何形状变化较小的区域主应力方向基本没有变化（见图4i）。

von Mises应力在盐体突起处增加（见图3b左侧虚线框），在盐体洼陷处减小了（见图3b右侧虚线框）；同样地，平均应力也是在盐体突起处增加，在洼陷处减小（见图4e）。已有研究表明，近盐处的安全钻井位置通常是在最小主应力大、平均应力大和剪应力小的区域[34-40]。因此难以确定盐体的突起和洼陷处是否为合适的钻井位置，因为突起处最小主应力大，同时von Mises应力也大；洼陷处的von Mises应力大，但最小主应力却小。

本文通过穿过盐体的垂直剖面图来进一步分析应力扰动情况。剖面1穿过左边的盐体（盐体都是凸状）（见图6a），在盐体突起附近，相较于远场应力（初始应力），水平应力、法向应力和垂向应力都明显增加；且水平应力和法向应力在接近盐体位置增加最大，垂向应力则在盐体突起处开始增加，在远离盐体一段距离后应力值达到最大（见图6a）。在盐体内，法向应力和水平应力都发生了减小且大小等于垂向应力，产生了各向同性的应力状态（见图6a）。同时，可以看到在剖面位置处剪应力特别小，基本为零（见图6a）。结果表明剖面上的盐体应力松弛之后的主应力方向相较于远场主应力方向无显著变化（见图6a、图6b）。在图6c中，盐体内von Mises应力和平均应力都相较于远场应力有所减小，但是von Mises应力减小更加明显，这是因为盐体内呈现一个各向同性的应力状态，不存在偏应力。沉积物中von Mises应力和平均应力在靠近盐体处达到最大值，随着远离盐体，应力逐渐减小；底面边界附近的von Mises应力甚至小于远场应力。

图6 穿过突起盐体的应力剖面图

剖面 2穿过右边的盐体（盐体都是下陷的）（见图7a），在盐体洼陷附近，相较于远场应力，水平应力、法向应力和垂向应力都有减小，水平应力减小明显（3层沉积物都可看到），垂向应力和法向应力变化不明显，只在中、下层沉积物中可见；水平应力的变化表现为在盐体附近应力最小，随着远离盐体其应力值逐渐增加。在此整个剖面的剪应力相较于其远场应力没有明显变化，也说明了在盐体应力松弛之后的主应力方向相较于远场主应力方向基本无变化（见图7a、图7b）。与剖面1处应力扰动不同的是，剖面2处的沉积物中 von Mises应力和平均应力在靠近盐体处达到最小值，随着远离盐体，应力逐渐增大（见图7c）。盐体内的各个应力变化与剖面1中类似，不再赘述。盐体水平时（无几何形状起伏），在盐体应力松弛过程中对周围沉积物的应力扰动很小甚至没有。为了进一步探究盐体的起伏程度与应力扰动之间的关系，本文通过改变盐体起伏程度的大小，设计两个对比模型。

图7 穿过洼陷盐体的应力剖面图

3.2 对比模型

3.2.1 对比模型1

为了与参考模型进行比较，在对比模型1，左边突起盐体的起伏更平缓（见图8a虚线框），其他的模型设置均与参考模型相同。

对比参考模型中的应力扰动结果（见图4、图8），发现此模型中的盐体应力松弛产生的应力扰动模式与参考模型一致：盐体突起处水平应力、法向应力、垂向应力、主应力、平均应力和von Mises应力均增加，洼陷处这几个应力均变小；突起左侧和洼陷右侧出现挤压型剪应力，突起右侧和洼陷左侧出现拉张型剪应力。所不同的是，由于左侧突起的盐体起伏更平缓，使得对沉积物中的应力扰动比参考模型中的小（见图4a—4c、图8a—8c）。

图8 模型1中盐体周围沉积物中的应力扰动

穿过盐体的剖面图可以更清楚地表现出盐体起伏程度对应力扰动的影响（见图9）。结果显示，盐体内的应力变化趋势与参考模型中的应力变化基本一致（见图6、图9）。在盐体突起附近的沉积物中，水平应力、法向应力和垂向应力的扰动均小于参考模型中对应的应力扰动。以水平应力为例，盐体与沉积物交界面处，参考模型的应力增加（σh=53.1 MPa，见图6a A点）要大于本模型（σh=44.5 MPa，见图9a B点）。此外，本模型剖面处的剪应力也无显著变化，基本为0（见图9a）。同样地，主应力、平均应力、von Mises应力的扰动与参考模型相比，也是因为盐体突起更平缓而有所减小（见图6b—6c、图9b—9c）。

3.2.2 对比模型2

图9 穿过突起盐体的应力剖面图

相比于参考模型，此模型的特点是，只有模型右边洼陷盐体的起伏发生了改变（起伏更陡峭），其他的模型设置部分均与参考模型一致。

将本模型的应力扰动结果与参考模型相比较（见图4、图10）可以看出，由于本模型中盐体洼陷得更加陡峭，水平应力、垂向应力和法向应力都减小得更显著（见图10a—10c）。由图11可以看出，本模型的水平应力、法向应力和垂向应力变化曲线与参考模型基本一致，所不同的是在中间沉积层中，水平应力的变化比较明显（见图11a），且在靠近盐体的地方，水平应力下降更多（见图7a、图11a）；剪应力同远场没有变化，接近于0。同参考模型相比，因为洼陷盐体处的起伏变陡，在中间沉积层中，最大主应力的变化比较明显（见图11b），而且在靠近盐体的地方，最大主应力下降得更多（见图11b）；类似地，本模型中的平均应力和 von Mises应力在盐体洼陷附近的下降也要比参考模型多。

3.3 克拉苏构造带西剖面盐构造模型

本文分析了几种不同几何样式的概念性模型的应力扰动，接下来分析实际剖面的应力扰动结果。模型材料参数和边界条件等均与参考模型相同。

类似于参考模型中的von Mises应力分析结果（见图3），由于盐体内的应力松弛，盐体内的剪应力随着时间的发展逐渐松弛减小，并在盐体周围的沉积层产生应力扰动（见图12b）。模型最终状态的应力（见图12c）是初始状态的应力（见图12a）与由于盐体的黏弹性应力松弛产生的应力扰动之和。

模型模拟结果显示，在此盐盆地中（相当于大区域的盐体洼陷），盐体上表面水平应力变小（见图13a），但是垂向应力在盆地中央区域（盐体平缓）有略微增加，往盆地两侧（左右两侧存在尺度较小的洼陷）开始出现垂向拉张应力（见图13b），说明在大尺度的盐体构造中，应力扰动要同时受到整个盐构造与局部小尺度盐构造的双重影响。盐体下表面呈凸状，此处出现了水平挤压应力，但是随着远离突起核部，应力逐渐减小；垂向应力在此处也具有类似的情况；法向应力的扰动与也与水平应力类似，只是应力扰动的值更小（见图13a、图13c）。在图13d中，剪应力的变化主要出现在盐体几何形状异常的附近。

von Mises应力在盐体突起处增加，如在盐体下表面突起处（见图13i中间虚线方框）、北秋里塔格背斜处（见图13i左边虚线框）以及北吐孜玛扎盐丘附近（见图13i右边虚线框）；平均应力在盐体下表面突起处和北秋里塔格背斜处也变大，但是在北吐孜玛扎盐丘附近变小）（见图13e），因此难以确定北秋里塔格背斜处是否是安全的钻井位置（von Mises应力和平均应力都变大），但是确定北吐孜玛扎盐丘附近不是合适的钻井位置（von Mises应力变大和平均应力减小）。

4 讨论

本文设计并模拟了 3个洼陷与突起组合的双层盐体模型，发现盐体几何形状与盐体周围沉积物中的应力扰动有很大关联，这与之前的很多研究结论一致[15-16,41]。模型模拟结果显示，在一个挤压背景应力场下，突起的盐体和洼陷的盐体所产生的应力扰动是相反的（挤压和拉张），而且突起和洼陷的曲率变化会对应力扰动有影响。模型的结果也适用于实际的盐构造模型（塔里木盆地库车坳陷盐构造）。

钻井稳定性分析要综合考虑最小主应力、平均应力和von Mises应力的扰动[15,42-45]。根据Mohr-Coulomb或者Drucker-Prager准则，平均应力的增加意味着材料可以承受更大的剪应力，材料越不容易产生变形或者破裂。von Mises应力的增加则意味着材料更加接近屈服。所以理想的钻井轨迹的设计应该是在平均应力大且von Mises应力小的地方。本文结果显示洼陷起伏较小盐体处是比较理想的钻井位置，因为和远场应力相比，von Mises应力减小得很多，而平均应力减小得不大（见图3b、图4e）。

图10 模型2中盐体周围沉积物中的应力扰动

图11 穿过洼陷盐体的应力剖面图

图12 von Mises应力扰动结果

本文的二维平面应变模型也存在一些局限性，比如，没有考虑第三维的盐体几何形状与运动；局部的温度异常会产生温度应力[46]；盐体周围的断裂带也会导致应力和孔隙压力的扰动。地质环境下，孔隙流体压力的变化会引起有效应力的变化和钻井设计，如果考虑这个因素，钻井位置可能会变化，本文没有包含孔隙流体压力的模拟[15]。因此，在分析实际钻井设计时，需要注意本研究模型模拟结果的适用性。随着地质与地球物理研究的发展，对上述因素有了更多的资料后，可以进一步模拟研究这些因素对盐体周围应力扰动的影响和作用，构建一个更加接近现实的模型，模拟并得到一个更加现实的应力场和孔隙流体压力场，更好地指导钻井设计。

图13 库车坳陷盐体周围沉积物中的应力扰动

5 结论

首先使用地质力学建模方法设计了 3个突起、洼陷相互组合双层盐体模型和一个塔里木盆地库车坳陷盐构造模型，并模拟了盐体周围的应力扰动情况。其次研究了挤压构造应力环境下盐体的应力松弛对其周围沉积物所产生的应力扰动，得到如下结论：

盐体的应力扰动与盐体的几何形状有重要的关系。在挤压背景应力场下，突起的盐体附近会产生挤压型应力扰动；洼陷的盐体附近会产生拉张型应力扰动；水平状盐体（无几何形状异常）周围无应力扰动。

盐体突起和洼陷的起伏程度变化也会影响应力扰动。突起或洼陷的起伏程度越大，应力扰动越大；反之，突起或洼陷的起伏程度越小，应力扰动越小。

库车坳陷盐体构造中，应力扰动同时受到整体盐构造几何形态与局部小尺度盐构造几何形态的双重影响。如果不考虑孔隙流体压力时，最佳的钻井位置是在平均应力大与von Mises应力小的位置，即本文中盐体下陷处；如果考虑孔隙流体压力，最优的钻井位置可能会变化。

符号注释：

D——弹性材料矩阵；{dε}——应变增量，无因次；{dεe}——黏性应变增量，无因次；{dεv}——弹性的应变增量，无因次；E——弹性模量，GPa；Q——黏度相关的矩阵；t——时间，s；σt——应力增量，MPa；σh——水平应力，MPa；σv——垂向应力，MPa；σVM——von Mises应力，MPa；σoop——法向应力，MPa；σ1——最大主应力，MPa；σ2——中间主应力；σ3——最小主应力，MPa；η——黏度，Pa·s；v——泊松比，无因次；σij——应力分量，Pa；gi——重力加速度，m/s2；ρ——密度，kg/m3。