基于热力耦合分析的火山热膨胀型裂缝定量表征

Quantitative Characterization of Fractures under Volcanic Thermal Expansion Based on Thermal-Mechanical Coupling Analysis

摘要

渤南地区是渤海海域近年来油气勘探开发的热点，新生界火山活动为油气成藏提供优势条件的同时，也给油气田开发带来了难题.内部非均质、高温剧烈的火山通道影响了周缘地层的裂缝发育，引发了严重的裂缝性井漏，造成巨大的经济损失，亟需对火山周缘膨胀型裂缝进行有效表征.根据测井、岩心、地震等资料，识别刻画火山岩岩相特征，分析周缘地层裂缝发育情况；通过岩心矿物分析、薄片观察等手段，确定火山活动温度场边界条件；基于能量守恒定律和地质力学，引入热力耦合原理，采用有限元数值模拟方法对火山活动引起的岩石热膨胀裂缝进行定量预测.基于热力耦合分析火山通道对周缘裂缝发育的影响，确定了火山对井漏风险的影响范围，对渤海油田钻井安全、火山相关油藏预测具有重要参考价值.

关键词

火山 / 热力耦合 / 裂缝定量预测 / 热膨胀 / 渤中34-9油田 / 石油地质学

Key words

volcano / thermal coupling / quantitative fracture prediction / thermal expansion / Bozhong 34-9 oil field / petroleum geology

0 引言

随着油气勘探技术向着更广、更深和更复杂的领域发展，火山发育区逐渐成为国内外油气勘探的重要领域（王月蕾等，2014；Shi et al.，2017；崔宝文等，2018）.火山是一种形态不规则、内部非均质、能量屏蔽性强的特殊地质体，其活动剧烈、内部高温，火山活动对油气成藏的建设或破坏改造具有重要意义（何家雄等，2010；林向洋等，2011；戴世立等，2012；Xie et al.， 2021）.火山活动带来的热液活动催化了烃源岩的成熟和生烃作用；产生的气孔、溶蚀孔和构造缝可以作为良好的储集空间和渗滤通道，这类裂缝型储层与砂岩孔隙型储层的不同之处是在致密、性脆、质纯的火成岩中发育各种构造裂缝（肖敦清等，2003；王维等，2017；刘晓健等，2020；尹志军等，2020）.火山岩与围岩物性的明显差异提供了良好的盖层条件和侧向封堵条件.而火山岩岩相复杂、横向结构变化快，剧烈的高温活动也对周围岩石造成破坏变形；不仅增加了火山岩岩相结构的研究难度，其超高温产生的周缘裂缝极易导致井漏发生，大大增加了油田钻井作业难度.

渤南地区黄河口凹陷南部斜坡带广泛分布了新生代火山岩，尤其是孔店组、沙河街组、东营组和馆陶组为岩浆活动的主要时期.斜坡带一般地层展布稳定，面积大，油气成藏条件差；而渤中34-9油田新生代火山活动频发，岩浆沿断裂大量喷涌，发育的大范围溢流相火山岩作为优质盖层，穿层的火山通道相对斜坡带油气的横向节流作用，其中古近系的油层多分布在东三段和沙一、二段，埋藏深度在 2 435.6~3 268.3 m，形成了大中型油田（王维等，2017；杨海风等，2018，2021；胡志伟等，2019；刘晓健等，2020；尹志军等，2020）.然而火山对油气圈闭的建设性改造作用为渤中地区提供优越成藏条件的同时，也给油田的开采钻探带来了诸多难题.火山岩岩相空间分布复杂、影响范围大，尤其火山通道相岩体强度大、贯穿性强，岩浆上涌时超高温烘烤围岩，在周围形成温度场，使岩石发生不同程度的膨胀变形和破裂，在钻井时易引发裂缝性井漏，带来了不同程度的经济损失.

以渤中34-9油田为例，裂缝性井漏发生在火山通道相周围，本文以火山通道相围岩裂缝的定量预测为研究目标，探讨火山活动对构造应力场及裂缝发育的影响.针对火山相关裂缝，前人利用地震属性（Sun et al.，2014；高宇慧等，2021；任涛等，2021）、成像测井（张兆辉等，2018）、应力场（李玉玺，2019）等方法定量描述火山岩裂缝，对火山活动造成的围岩裂缝预测研究不深.随着油气勘探开发领域的进一步扩展，井漏安全问题、火山岩相关油藏预测亟需火山通道相围岩裂缝定量预测研究.本文结合渤中34-9油田地震、钻井、测井等资料，针对引发裂缝性井漏的火山通道相进行刻画，基于能量守恒定律和地质力学，引入热力耦合原理，采用有限元数值模拟方法对火山活动引起的岩石热膨胀裂缝进行定量预测，以期为渤海油田钻井安全、火山相关油藏预测提供参考依据和借鉴.

1 区域地质概况

渤海盆地是由不同级序的旋扭应力场控制形成的不同级序帚状构造组成，渤中凹陷紧邻济阳坳陷，一起组成了以鲁西隆起为砥柱的冀鲁帚状构造体系内旋区，是一个西南收敛、东北发散的大型帚状构造体系.渤海海域南部处于帚状构造末端，油气富集，储量丰富.渤中34-9油田位于渤海海域南部，黄河口凹陷与莱北低凸起过渡的斜坡带上，西邻郯庐断裂带西支，发育了以EW向和NE-SW向为主的复杂断裂体系，整体上形成了南超北断的半地堑式箕状构造特征.钻探含油目的层为古近系东营组‒沙河街组，主要发育辫状河三角洲前缘沉积亚相，砂岩含量高.根据对测井资料的解释，沙河街组储层具有中孔、中渗的物性特征，东营组为中孔、中‒高渗的特征，砂泥岩与火成岩形成一套储盖组合，且储层分布稳定、连通性好，研究区近几年发现油气储量近8 000×10⁴ m³（杨海风等，2018）.自油田发现以来钻井皆在含油层系钻遇厚度不一的火山岩，岩性非均质性强，裂缝发育复杂，钻井过程中半数井发生裂缝性漏失，为钻井施工和油气开发带来了巨大的挑战.

渤中凹陷经历了多期构造运动，是渤海湾盆地地壳隆升、岩浆上涌活动规模最大的区域.古新世‒始新世时期，深部岩浆上涌，地层隆起抬升，在热底辟作用下岩浆侧向流动，其产生的拖曳力使地壳拉伸减薄（滕吉文等，1997），在渤中凹陷南侧毗邻区发生张裂、断陷，发育了一系列近EW向基底断裂，形成了黄河口凹陷，应力场表现为近SN向的伸展作用；始新世‒渐新世时期，由于太平洋板块俯冲作用，在郯庐断裂带的影响下，渤中凹陷及其周边应力场表现为右旋拉张走滑，基底断裂继续活动，发育了一系列近EW向和NE-SW向走滑断裂，地幔岩浆沿走滑断裂侵入地壳、喷出地表，形成了诸多大型火山通道相和范围较广的溢流相火山岩；新近纪‒第四纪，岩浆活动明显减弱，在青藏高原向东滑移及太平洋板块俯冲的共同作用下，现今应力场表现为近EW向挤压弱走滑作用.

火山活动产生的温度一般在1 000 ℃左右，在超高温下岩石破裂现象是火山岩裂缝研究一直以来忽略的问题，但对于花岗岩、砂岩、油页岩等岩石热开裂是公认的热点问题（于永军，2015；薛东杰等，2016；史伟新，2018）；且在火山野外露头中也存在密集的膨胀型张性缝，前人根据野外火山露头模拟了火山熔岩的热作用，对其应力变化和裂缝造成的次级灾害进行了评估（李镇村等，2016）.在针对渤海第三系裂缝性井漏研究中发现，裂缝发育特征与火山通道及火山岩分布具有较强相关性，火山通道周围裂缝增多，井漏情况更为复杂.火山活动在始新世‒渐新世时期活动范围大，火山岩分布广，对裂缝发育及裂缝性井漏影响最大，是本次研究的重点.

2 火山热力作用的原理及方法

2.1　火山岩岩相特征

利用火山岩体非均质性和能量屏蔽性强的特征，通过常规地震剖面反射特征和方差体地震属性，将火山岩体从围岩中识别出来（朱红涛等，2014）.研究区火山主要以裂隙式喷发和中心式喷发为主，中心式火山喷发通过近似垂直的通道与岩浆库相连而喷出地表，在平面成孤立的点式分布；裂隙式火山岩浆沿断层上升喷出地表，以玄武质岩浆为主，岩浆粘度较低，不具爆炸性，在平面上沿断层呈串珠状分布（图1）.

研究区火山岩岩相主要表现为火山通道相、火山溢流相和爆发相.火山通道相内部岩石成层性差，振幅较弱，反射杂乱，剖面上表现为倒锥状，裂隙式火山通道相上部略宽，呈蘑菇状，平面上方差体属性表现为似圆形的高值区；火山溢流相在地震剖面上成层性好，振幅较强，反射连续，一般与火山通道相临近或相接.

通过钻井岩心和薄片观察鉴定，研究区火山岩储层多为中基性火成岩，以裂隙式或宁静式火山为主，裂隙式火山溢流相分布范围广、厚度及成分较为稳定，成层性好，且沿断裂发育一系列蘑菇状火山通道（图2），为油气储藏提供了良好的盖层和侧向封堵条件.故油气开采钻井多设计在裂隙式火山发育区，是热膨胀裂缝研究的重点.

火山喷发期次的岩性划分原则为以溢流相玄武岩作为喷发期次的开始，以火山沉积岩或正常沉积岩作为喷发期次的结束.根据岩性和测井曲线（图2），始新世‒渐新世时期火山主要发育在东营组，东二下段沉积时期可划分为2个旋回，东二上段沉积时期2个旋回，东一段沉积时期2个旋回.由地震剖面解释可知，BZ34-9-E井附近有两个火山通道，根据深浅双侧向测井数据计算，东一段‒东二下段地层的裂缝密度和裂缝开度有明显的突增特征，其突增趋势与火山喷发旋回相吻合；且对比火山发生喷发期次与裂缝密度曲线，火山活动都伴随着裂缝突然密集发育的情况，这证明了火山岩浆活动的确会对裂缝的发育造成影响（王启明等，2020；Ye et al.，2022）.

2.2　火山周缘温度场的确定

要模拟火山周围的热力作用，首先要确定火山活动时的温度，从而设置研究区模型的温度边界条件，对火山产生的温度进行模拟.多口钻井岩心显示，钻遇火山通道相储层岩石类型主要为斜长石、辉石等矿物组成的玄武岩（图3）.斜长石的结晶温度为1 050~1 228 ℃，辉石的结晶温度为1 170~ 1 215 ℃，熔点在1 300 ℃以上（靳是琴等，1987）.斜长石、辉石属于基性端元组分，是岩石中较早结晶的产物，因此斜长石、辉石结晶温度高于岩石的成岩温度，而与岩浆源区温度接近或略低，本次模拟火山底部岩浆温度边界条件取1 200 ℃.

火山喷发时，岩浆沿火山通道上升，岩浆与围岩直接接触产生热传导作用，由于二者的温度差，热能直接以原子震动的方式由高温处向四周低温处传导.热传导遵循傅里叶定律 q^n = - k \frac{\mathrm{d}T}{\mathrm{d}x}，通过热传导在通道周缘围岩中形成一个温度场，周围岩石受热膨胀，岩石自由热膨胀形变产生的附加应力改变了周缘应力场，围岩裂缝十分发育，致使井漏风险大大提高.

选取研究区火山通道相较为集中的区域为靶区，以火山机构特征在地震剖面上较为明显且裂缝发育较多的东二段为目的层，通过对方差体属性的识别，结合测井数据元素相关性分析，对比刻画火山岩范围，提取东二段不同深度的火山通道切片形状.采用有限元建模方法，基于多个火山通道切片进行火山通道相模型构建（图4），并通过布尔运算等方法嵌入地层模型中，构建研究区有限元模型，模型划分为22 294个节点和108 820个单元（图5）.本文主要对热膨胀作用产生破裂现象的机理和影响范围进行研究，且火山通道相影响范围有限，研究区地层发育较为均质，故有限元分析中不同岩石均设置各向同性的导热系数和热膨胀系数.

通过火山岩矿物分析，进行岩浆源区温度边界条件设置.火山通道底部约束温度设为1 200 ℃，理论上距离火山无限远的地层温度为当时地层普通埋藏温度，设置足够远的围岩边界为平均埋藏温度，考虑当时埋深，设置边界温度为40 ℃.基于ANSYS热分析模块，对靶区内火山造成的温度场进行模拟.以几个火山为中心，温度呈放射状向四周减弱，沿断裂形成的大型火山影响的范围更大，小型火山影响范围小，其与应力场影响的强弱有关.

2.3　火山热应力作用模拟

火山通道内的岩浆产生的高温使周围岩石温度升高，产生了一个放射状温度场，在温度场的作用下，岩石发生不同程度的膨胀变形，由热膨胀作用产生的附加应力形成了一个热应力场.

热应力场的应力强度高值区主要集中在火山通道周围地区，是岩石变形破裂的优势区域，也是最大主应力的高值区；各火山之间地区由于温度较高，也是应力集中的区域.模型剖面显示，在温度场的作用下，火山周围发生明显形变，地层上拱，应力强度由火山向两侧逐渐减弱.由最大主应力方向可知，应力方向呈放射状向周围递减.

3 基于热力耦合作用的火山围岩裂缝模拟

3.1　热力耦合作用的本构关系描述

火山活动作用是涉及多个物理场互相耦合的复杂问题，本文研究的火山是在构造应力场作用下，沿活动的走滑断裂上涌而形成的；且在拉张应力环境下，岩浆对地层上拱影响较小，等效于构造应力场作用中的主要研究重点在于应力场和热应力场的耦合.在温度场和岩石变形场耦合的问题中，热膨胀作用使应力场改变而造成周围岩体破裂，是本次研究的中心问题，故只研究热膨胀对应力场的单向耦合作用，不考虑应力场的改变反过来对温度场的影响.

岩石在一定载荷下发生形变到破坏，是一个损伤积累并伴随能量部分耗散的过程，现于热力耦合作用下的岩石破裂过程引入损伤变量，用来反映岩石裂缝发育过程中的应力‒应变关系.基于广义胡克定律，岩石在发生破裂损伤时，其内部的应力σ和应变ε服从以下关系：

\left\{\begin{array}{c}{\epsilon }_{\mathrm{1}}=\frac{\mathrm{1}}{\left[\mathrm{1}-D\right]E}\left[{\sigma }_{\mathrm{1}}-\mu \left({\sigma }_{\mathrm{2}}+{\sigma }_{\mathrm{3}}\right)\right]\\ {\epsilon }_{\mathrm{2}}=\frac{\mathrm{1}}{[\mathrm{1}-D]E}\left[{\sigma }_{\mathrm{2}}-\mu \left({\sigma }_{\mathrm{1}}+{\sigma }_{\mathrm{3}}\right)\right]\\ {\epsilon }_{\mathrm{3}}=\frac{\mathrm{1}}{[\mathrm{1}-D]E}[{\sigma }_{\mathrm{3}}-\mu \left({\sigma }_{\mathrm{1}}+{\sigma }_{\mathrm{2}}\right)]\end{array}\right.，（1）

式中，D为损伤因子；E为弹性模量；μ为泊松比.

岩石在荷载因素作用下产生的不可恢复主应变为岩石在破坏过程中的唯一耗能机制，则将岩石在破坏开始时某点单元体的x，y，z三个方向耗能率叠加：

\psi = \sum {\sigma }_i {\dot{\epsilon }}_i^N = - \frac{\dot{D}}{{\left[\mathrm{1}-D\right]}^{\mathrm{2}}E} [ {\sigma }_{\mathrm{1}}^{\mathrm{2}} + {\sigma }_{\mathrm{2}}^{\mathrm{2}} + {\sigma }_{\mathrm{3}}^{\mathrm{2}} - \mathrm{2} \mu \left({\sigma }_{\mathrm{1}}{\sigma }_{\mathrm{2}}+{\sigma }_{\mathrm{2}}{\sigma }_{\mathrm{3}}+{\sigma }_{\mathrm{1}}{\sigma }_{\mathrm{3}}\right) ]，（2）

式中： {\sigma }_i ( i = \mathrm{1,2} , \mathrm{3} )为破坏刚开始时该点的名义主应力， {\dot{\epsilon }}_i^N ( i = \mathrm{1,2} , \mathrm{3} )为由损伤引起的不可恢复主应变率，为式（1）对D的导数.

在平面应力状态下， {\sigma }_{\mathrm{2}} = {\sigma }_{\mathrm{3}}的约束条件为：

F \left({\sigma }_{\mathrm{1}},{\sigma }_{\mathrm{3}}\right) = \frac{\mathrm{1}}{\mathrm{4}G} [ \frac{{\sigma }_{\mathrm{1}}-{\sigma }_{\mathrm{3}}}{\mathrm{c}\mathrm{o}\mathrm{s}\phi } - ( {\sigma }_{\mathrm{1}} + {\sigma }_{\mathrm{3}} {)\mathrm{t}\mathrm{a}\mathrm{n}\phi ]}^{\mathrm{2}} + \alpha \frac{({\sigma }_{\mathrm{1}}+\mathrm{2}{\sigma }_{\mathrm{3}}{)}^{\mathrm{2}}}{\mathrm{18}K} + \beta，（3）

式中： \alpha为岩石初始屈服的剪切应变能， \beta为岩石初始屈服的体积应变能，G为剪切模量，K为体积模量.

基于最小耗能原理，引用拉格朗日方程为：

\frac{\partial [\psi +\lambda F\left({\sigma }_{\mathrm{1}}-{\sigma }_{\mathrm{3}}\right)]}{\partial {\sigma }_i} = \mathrm{0}， i=1,3，（4）

故损伤因子演化方程为：

D = \mathrm{1} - \mathrm{e} \mathrm{x} \mathrm{p} [ \frac{{\lambda }_{\mathrm{1}}\left(A_{\mathrm{2}}{B}_{\mathrm{1}}-A_{\mathrm{1}}{B}_{\mathrm{2}}\right){\sigma }_{\mathrm{3}}{\epsilon }_{\mathrm{1}}}{A_{\mathrm{2}}{\epsilon }_{\mathrm{1}}-\mathrm{2}{A}_{\mathrm{1}}{\epsilon }_{\mathrm{3}}} + c_{\mathrm{0}} ]，（5）

式中： A_{\mathrm{1}} = \frac{{(\mathrm{1}-\mathrm{s}\mathrm{i}\mathrm{n}\phi )}^{\mathrm{2}}}{\mathrm{2}G{(\mathrm{c}\mathrm{o}\mathrm{s}\phi )}^{\mathrm{2}}} + \frac{\alpha }{\mathrm{9}K}， A_{\mathrm{2}} = \frac{\mathrm{2}\alpha }{\mathrm{9}K} - \frac{\mathrm{1}}{\mathrm{2}G}， B_{\mathrm{1}} = \frac{\mathrm{1}}{\mathrm{2}G} - \frac{\mathrm{2}\alpha }{\mathrm{9}K}， B_{\mathrm{2}} = \frac{-{(\mathrm{1}+\mathrm{s}\mathrm{i}\mathrm{n}\phi )}^{\mathrm{2}}}{\mathrm{2}G{(\mathrm{c}\mathrm{o}\mathrm{s}\phi )}^{\mathrm{2}}} - \frac{\mathrm{4}\alpha }{\mathrm{9}K}.

在温度作用下岩石发生热膨胀变形，变形产生的附加应力即为热应力，当构造应力场与温度场产生的热应力存在单向耦合作用时，热力耦合作用的岩石变形存在受热自由膨胀的应变量.该应变量是在特定的温度变化值下，通过岩石热膨胀系数 \alpha计算得到，对于各向同性岩石，3个方向的热膨胀应变皆可表示为：

{\epsilon }_{\mathrm{1}}^T = {\epsilon }_{\mathrm{2}}^T = {\epsilon }_{\mathrm{3}}^T = \alpha ∆ T，（6）

式中，α为岩石3个方向上相等的热膨胀系数； ∆ T为温度变化值.则热力耦合作用下，岩石材料的总变形量为：

\left\{\begin{array}{c}{\epsilon }_{\mathrm{1}}^H={\epsilon }_{\mathrm{1}}+{\epsilon }_{\mathrm{1}}^T=\frac{\mathrm{1}}{\left[\mathrm{1}-D\right]E}\left[{\sigma }_{\mathrm{1}}-\mu \left({\sigma }_{\mathrm{2}}+{\sigma }_{\mathrm{3}}\right)\right]+\alpha ∆T\\ {\epsilon }_{\mathrm{2}}^H={\epsilon }_{\mathrm{2}}+{\epsilon }_{\mathrm{2}}^T=\frac{\mathrm{1}}{[\mathrm{1}-D]E}\left[{\sigma }_{\mathrm{2}}-\mu \left({\sigma }_{\mathrm{1}}+{\sigma }_{\mathrm{3}}\right)\right]+\alpha ∆T\\ {\epsilon }_{\mathrm{3}}^H={\epsilon }_{\mathrm{3}}+{\epsilon }_{\mathrm{3}}^T=\frac{\mathrm{1}}{[\mathrm{1}-D]E}[{\sigma }_{\mathrm{3}}-\mu \left({\sigma }_{\mathrm{1}}+{\sigma }_{\mathrm{2}}\right)]+\alpha ∆T\end{array}\right.，（7）

其中， {\sigma }_{\mathrm{1}}、 {\sigma }_{\mathrm{2}}、 {\sigma }_{\mathrm{3}}为耦合应力， {\epsilon }_{\mathrm{1}}、 {\epsilon }_{\mathrm{2}}、 {\epsilon }_{\mathrm{3}}为耦合应力引起的应变.

基于损伤变量和广义胡克定律，在热力耦合作用下岩石的损伤本构方程为：

\left\{\begin{array}{c}{\sigma }_{\mathrm{1}}=\frac{E}{\mathrm{1}-E}\alpha ∆T\left(\mathrm{1}-D\right)+\mu \left({\sigma }_{\mathrm{2}}+{\sigma }_{\mathrm{3}}\right)\\ {\sigma }_{\mathrm{2}}=\frac{E}{\mathrm{1}-E}\alpha ∆T\left(\mathrm{1}-D\right)+\mu ({\sigma }_{\mathrm{1}}+{\sigma }_{\mathrm{3}})\\ {\sigma }_{\mathrm{3}}=\frac{E}{\mathrm{1}-E}\alpha ∆T\left(\mathrm{1}-D\right)+\mu ({\sigma }_{\mathrm{1}}+{\sigma }_{\mathrm{2}})\end{array}\right..（8）

3.2　热力耦合作用下的应力场

基于地质背景和构造形迹分析，始新世‒渐新世火山活动时期，研究区受郯庐断裂带和太平洋板块运动的右旋伸展走滑作用，形成了一系列NEE向为主的张性缝和张剪缝.岩石力学参数是岩石的基本物理性质，是应力场模拟及裂缝预测的基础，通过测井参数计算优选相关岩石力学参数计算模型，计算出单井动态岩石力学参数，结合调研和岩石力学实验计算静态力学参数，对单井连续的动态岩石力学参数进行动‒静态校正，得到连续的静态力学参数.在单井力学参数约束下，通过波阻抗、瞬时频率、方差体等地震属性，构建岩石力学参数体，提取目的层相关的力学参数并求取平均，分别赋予断层、火山通道和地层不同属性的岩石力学参数；并经过反复尝试加载，对靶区应力场进行模拟，从而确定了火山活动时期的力学边界条件（图8），具体为：模型底部约束，加载东、西两个面上42 MPa拉伸应力，以及 50 MPa的顺时针剪切走滑作用应力.赋予各单元不同的力学参数并精细划分网格后，基于热力耦合下的应力‒应变关系和有限元模拟原理，获得火山热膨胀和构造应力场耦合作用下的三维应力场模拟结果.

由于始新世‒渐新世期间的右旋伸展走滑作用，总体上应力高值区沿NW-SE方向长条形展布，中部火山对周围应力场影响较大，最大主应力和最小主应力分布趋势差异较大（图9）.最大主应力受火山膨胀作用影响，在火山周围表现为挤压应力，在区域右旋走滑作用下，南北两侧呈现椭圆状分布的高值区，由于EW方向的伸展作用，东西两侧为小范围低值区；最小主应力在火山周围表现为拉张应力，南北两侧表现为低值区，东西两侧表现为高值区，有利于张性缝或张剪缝的发育.火山周围应力强度表现较大，尤其与断层相连部位，高值区呈近E-W方向展布，是裂缝发育有利区域.

基于应力场模拟结果分析，火山热膨胀作用对周围地层的应力分布影响较为明显，应力场分布明显收到火山通道和断层的分割作用.BZ34-9-E井距离火山通道较近，所处区域应力强度大，应力较为集中，属于裂缝发育区，大大增加了井漏发生的风险.

3.3　热膨胀裂缝定量预测

根据地质力学相关理论，裂缝参数与应力场存在一定的定量关系.岩石在受力变形时，内部积累的应变能释放率达到产生裂缝体表面所需能量时，岩石发生破裂（季宗镇等，2010；冯建伟等，2011）.前人基于地质力学中的裂缝表面能及应变能理论（Feng et al.，2019， 2021），优选相关岩石裂准则，建立表征单元体裂缝模型，推导出裂缝密度和开度参数计算公式.

挤压应力条件下：

\left\{\begin{array}{c}{\sigma }_{\mathrm{P}}=\frac{\mathrm{2}{C}_{\mathrm{0}}\mathrm{s}\mathrm{i}\mathrm{n}\mathrm{2}\theta +\left(\mathrm{1}+\mathrm{c}\mathrm{o}\mathrm{s}\mathrm{2}\theta \right){\sigma }_{\mathrm{3}}}{\mathrm{1}-\mathrm{c}\mathrm{o}\mathrm{s}\mathrm{2}\theta }\\ D_{\mathrm{v}\mathrm{f}}=\frac{E\left({\sigma }_{\mathrm{1}}{\epsilon }_{\mathrm{1}}+{\sigma }_{\mathrm{2}}{\epsilon }_{\mathrm{2}}+{\sigma }_{\mathrm{3}}{\epsilon }_{\mathrm{3}}\right)-k{\sigma }_{\mathrm{P}}\left(k{\sigma }_{\mathrm{P}}-\mathrm{2}\mu \left({\sigma }_{\mathrm{2}}+{\sigma }_{\mathrm{3}}\right)\right)}{\mathrm{2}E\left({\mathrm{J}}_{\mathrm{0}}+{\sigma }_{\mathrm{3}}b\right)}\\ D_{\mathrm{l}\mathrm{f}}=\frac{\mathrm{2}{D}_{\mathrm{v}\mathrm{f}}{L}_{\mathrm{1}}{L}_{\mathrm{3}}\mathrm{s}\mathrm{i}\mathrm{n}\theta \mathrm{c}\mathrm{o}\mathrm{s}\theta -L_{\mathrm{1}}\mathrm{s}\mathrm{i}\mathrm{n}\theta -L_{\mathrm{3}}\mathrm{c}\mathrm{o}\mathrm{s}\theta }{L_{\mathrm{1}}^{\mathrm{2}}\mathrm{s}\mathrm{i}\mathrm{n}\mathrm{2}\theta +L_{\mathrm{3}}^{\mathrm{2}}\mathrm{c}\mathrm{o}\mathrm{s}\mathrm{2}\theta }\\ b=\frac{\epsilon f}{D_{\mathrm{l}\mathrm{f}}}=\frac{\left|{\epsilon }_{\mathrm{3}}\right|-\left|{\epsilon }_{\mathrm{0}}\right|}{D_{\mathrm{l}\mathrm{f}}}\end{array}\right..（9）

拉张应力条件下：

\left\{\begin{array}{c}{D}_{\mathrm{v}\mathrm{f}}=\frac{E\left({\sigma }_{\mathrm{1}}{\epsilon }_{\mathrm{1}}+{\sigma }_{\mathrm{2}}{\epsilon }_{\mathrm{2}}+{\sigma }_{\mathrm{3}}{\epsilon }_{\mathrm{3}}\right)-{\sigma }_{\mathrm{t}}^{\mathrm{2}}}{\mathrm{2}E\left({\mathrm{J}}_{\mathrm{0}}+{\sigma }_{\mathrm{3}}b\right)}\\ D_{\mathrm{l}\mathrm{f}}=\frac{\mathrm{2}{D}_{\mathrm{v}\mathrm{f}}{L}_{\mathrm{1}}{L}_{\mathrm{3}}\mathrm{s}\mathrm{i}\mathrm{n}\theta \mathrm{c}\mathrm{o}\mathrm{s}\theta -L_{\mathrm{1}}\mathrm{s}\mathrm{i}\mathrm{n}\theta -L_{\mathrm{3}}\mathrm{c}\mathrm{o}\mathrm{s}\theta }{L_{\mathrm{1}}^{\mathrm{2}}\mathrm{s}\mathrm{i}\mathrm{n}\mathrm{2}\theta +L_{\mathrm{3}}^{\mathrm{2}}\mathrm{c}\mathrm{o}\mathrm{s}\mathrm{2}\theta }\\ b=\frac{\left|{\epsilon }_{\mathrm{3}}\right|-\left|{\epsilon }_{\mathrm{t}}\right|}{D_{\mathrm{l}\mathrm{f}}}\end{array}\right. , \mathrm{或} D_{\mathrm{l}\mathrm{f}} = D_{\mathrm{v}\mathrm{f}} .（10）

式中： D_{\mathrm{v}\mathrm{f}}为构造裂缝体密度， D_{\mathrm{l}\mathrm{f}}为构造裂缝线密度， {\sigma }_{\mathrm{P}}为最大主应力临界值， b为裂缝开度， C_{\mathrm{0}}为岩石的内聚力， \theta为岩石破裂角， k {\sigma }_{\mathrm{P}}对应产生微裂缝时的最大主应力（k为常数）， {\sigma }_{\mathrm{t}}为岩石抗张强度， {\epsilon }_{\mathrm{t}}为岩石弹性变形的最大张应力， {\epsilon }_{\mathrm{0}}为产生微裂缝时的最小主应变， L_{\mathrm{1}}、 L_{\mathrm{3}}为最大主应力、最小主应力方向的单元网格长度.

通过对ANSYS软件的二次开发，导入裂缝参数计算模型，基于耦合应力场模拟结果，计算裂缝密度和开度参数（图10），依据有限元模型划分单元长度100 m，裂缝参数预测结果的分辨率在 100 m左右.在伸展走滑作用下，火山通道周围地层裂缝密度的高值区呈NE-SW方向长椭圆状展布，断层周围也是裂缝发育高值区；裂缝开度高值区分布不连片，在断层和火山周围开度最大，以火山为中心向外逐渐减小.BZ34-9-E井区受火山通道影响裂缝明显较为发育，裂缝密度在2条/m左右，裂缝开度为0.1 mm左右，井漏风险大.从平面上可知，在火山热膨胀作用下，耦合应力对围岩裂缝发育明显，给钻井生产带来较大的风险.

4 讨论

在火山沿断裂发生喷发活动时，岩浆与围岩接触发生的热传递，基于岩石的热传递效果，在地层中产生以火山为温度高点的温度场.岩石在不同温度下，发生不同程度的膨胀作用，膨胀变形使岩石之间发生挤压作用，产生以火山为中心放射状的应力场（图7），但应力值相对于区域地应力场较小；在地应力场叠加耦合后，对周围应力场的影响局限在一定范围内，故膨胀作用影响裂缝发育的区域也有限.在钻井工作时，可通过对地震剖面和钻遇岩心对火山通道相位置进行判断，利用井位设计或侧钻工艺使井轨迹避开裂缝发育区域，避免钻井漏失风险.

BZ34-9-E井钻至火山岩附近时发生裂缝性漏失，由深浅双侧向测井数据解释裂缝参数可知，裂缝多发育在火山喷发活动中间，东二段火山活动周期长，裂缝在垂向地层中间区域最为发育.由于研究区东营组具有砂泥互层岩性组合特征，对应力场扰动较大，且火山通道形状并不规律，模拟建模时需对其进行平滑处理，故预测的裂缝参数结果与实际数据会产生一定误差.通过对模拟裂缝密度参数体进行重采样，绘制相关曲线与测井解释数据相对比（图11），发育趋势与模拟结果趋势基本一致.

由平面裂缝模拟结果和过BZ34-9-E井剖面裂缝发育情况来看（图12），火山热膨胀作用影响周围岩石破裂的范围为200~300 m.一般来说，通过热力耦合模拟对火山通道周围裂缝进行模拟预测的工作较为复杂，在钻井工程应用不便，且火山通道形状及规模不同，对应力场和裂缝的影响也不同.利用古温度的可测性，使用温度场对裂缝影响范围进行标记，可简单对火山通道周围的裂缝发育范围进行预测，避开井漏高风险区，对打井靶点必经过火山通道周围的井进行周缘裂缝定量表征，对井漏高风险点提前做好堵漏或侧钻的准备.对比温度场与裂缝密度发育情况，地层温度在882.33 ℃以上时，裂缝密度在2.2条/m以上，裂缝基本在火山通道100 m范围内，是井漏风险最高的区域；地层温度在800 ℃以上时，裂缝发育密度在1.9条/m以上，在距通道300 m以内的范围内，主要受火山活动控制形成裂缝密度高值区.

在钻井设计时，通过地震同相轴识别和方差体属性刻画确定火山通道的位置和形状特征，尽量避开火山通道200~300 m的距离；若钻井靶点落在火山通道附近，在钻遇地层的温度为800 ℃、距通道100~300 m范围内时，应放慢钻速，调节钻井液压力，准备堵漏剂防止井漏的发生；若钻遇地层温度为882.33 ℃、距火山通道100 m以内范围时，在准备堵漏剂的同时，必要时进行侧钻工作.

5 结论

（1）始新世‒渐新世时期，研究区火山活动剧烈，根据方差体属性和测井双侧向解释，裂缝在火山通道发育程度高，且裂缝参数高值点与火山活动期相吻合，对井漏影响最大.

（2）渤中地区火山活动时，岩浆与围岩存在的温度差造成热传递作用产生温度场，岩石在不同温度下发生膨胀作用产生热应力场，与构造应力场耦合，基于损伤变量和广义胡克定律，构建了岩石的应力‒应变关系本构方程.

（3）热膨胀作用影响裂缝区域在距火山通道200~300 m范围内.火山通道周围地层温度在800 ℃以上时，漏失风险较高，若钻遇地层温度为882.33 ℃，距火山通道100 m范围时，裂缝密度在2.2条/m以上，漏失风险最高，应考虑堵漏和侧钻工作.

参考文献

[1] Cui, B.W., Meng, Q., Bai, X.F., et al., 2018. Petroleum Exploration Progress and Suggestions for North Songliao Basin. Petroleum Geology & Oilfield Development in Daqing, 37(3): 1-9 (in Chinese with English abstract).

[2] Dai, S.L., Jiang, C.J., Zhang, Y.G., 2012. Methodologies for Predicting Effective Reservoir of Volcanic in Anda Area of Xujiaweizi Depression. Geological Science and Technology Information, 31(2):67-71 (in Chinese with English abstract).

[3] Feng, J. W., Qu, J. X., Zhang, P. F., et al., 2021. Development Characteristics and Quantitative Prediction of Multiperiod Fractures in Superdeep Thrust-Fold Belt. Lithosphere, 2021(1): 1-28. https://doi.org/10.2113/2021/8895823

[4] Feng, J. W., Shi, S., Zhou, Z. H., et al., 2019. Characterizing the Influence of Interlayers on the Development and Distribution of Fractures in Deep Tight Sandstones Using Finite Element Method. Journal of Structural Geology, 123: 81-95. https://doi.org/10.1016/j.jsg.2019.03.009

[5] Feng, J.W., Dai, J.S., Ma, Z.R., et al., 2011. The Theoretical Model between Fracture Parameters and Stress Field of Low-Permeability Sandstones. Acta Petrolei Sinica, 32(4):664-671 (in Chinese with English abstract).

[6] Gao, Y.H., Li, S.Y., Ba, Z.C., et al., 2021. Characteristics and Quantitative Prediction of Carboniferous Volcanic Reservoir Fractures in Well Block Hong-172 of Junggar Basin. Petroleum Geology & Oilfield Development in Daqing, 40(2): 30-41 (in Chinese with English abstract).

[7] He, J.X., Zhu, Y.H., Weng, R.N., et al., 2010. Characters of North-West Mud Diapirs Volcanoes in South China Sea and Relationship between Them and Accumulation and Migration of Oil and Gas. Earth Science, 35(1):75-86 (in Chinese with English abstract).

[8] Hu, Z.W., Yang, H.F., Han, Z.J., et al., 2019. Characteristics of Volcanic Rocks and Their Control Actions on Hydrocarbon Accumulation of Bozhong 34-9 Oilfield in Bohai Sea. Geological Science and Technology Information, 38(3): 208-218 (in Chinese with English abstract).

[9] Ji, Z.Z., Dai, J.S., Wang, B.F., 2010. Quantitative Relationship between Crustal Stress and Parameters of Tectonic Fracture. Acta Petrolei Sinica, 31(1): 68-72 (in Chinese with English abstract).

[10] Jin, S.Q., Li, D.C., Li, X.Z., et al., 1987. The Metamorphic Characteristics of Granulite Facies Rocks in the Laixi-Pingdu Area in Eastern Shandong Province. Acta Geological Sinica, 61(3): 240-252 (in Chinese with English abstract).

[11] Li, Y. X., 2019. Characterization and Prediction of Fractures in Carboniferous Volcanic Rocks in Shishi Oilfield, Junggar Basin (Dissertation). Southwest Petroleum University, Chengdu (in Chinese with English abstract).

[12] Li, Z.C., Shi, Y.L., Zhang, Z.Y., et al., 2016. Calculation of Thermal Stress during the Cooling of Lava Tubes. Journal of University of Chinese Academy of Sciences, 33(3): 412-420 (in Chinese with English abstract).

[13] Lin, X.Y., Su, Y.P., Zheng, J.P., et al., 2011. Characteristics and Controlling Factors of Volcanic Reservior in Kelameili Area of Jungger Basin. Geological Science and Technology Information, 30(6):28-37 (in Chinese with English abstract).

[14] Liu, X.J., Wang, Q.B., Dai, L.M., et al., 2020. Reservoir Characteristics and Formation Mechanisms of Lacustrine Carbonate and Volcanics Mixing Sediments, Laizhouwan Sag. Earth Science, 45(10): 3579-3588 (in Chinese with English abstract).

[15] Ren, T., Qin, J., Zhou, Y., et al., 2021. Prediction of Volcanic Fractures in Che 476 Well Block in Northwestern Margin of Junggar Basin by Using Pre-Stack Azimuthal Anisotropy. Petroleum Geology and Engineering, 35(1):1-7 (in Chinese with English abstract).

[16] Shi, J. A., Sun, G. Q., Zhang, S. C., et al., 2017. Reservoir Characteristics and Control Factors of Carboniferous Volcanic Gas Reservoirs in the Dixi Area of Junggar Basin, China. Journal of Natural Gas Geoscience, 2(1): 43-55. https://doi.org/10.1016/j.jnggs.2017.03.001

[17] Shi, W. X., 2018. Experimental Study on Thermal Cracking of Fluid Saturated Rocks (Dissertation). China University of Petroleum, Qingdao (in Chinese with English abstract).

[18] Sun, Q. L., Wu, S. G., Cartwright, J., et al., 2014. Neogene Igneous Intrusions in the Northern South China Sea: Evidence from High-Resolution Three Dimensional Seismic Data. Marine and Petroleum Geology, 54: 83-95. https://doi.org/10.1016/j.marpetgeo.2014.02.014

[19] Teng, J.W., Zhang, Z.J., Zhang, B.M., et al., 1997. Geophysical Fields and Background of Exceptional Structure for Deep Latent Mantle Plume in Bohai Sea. Chinese Journal of Geophysics, 40(4): 468-480 (in Chinese with English abstract).

[20] Wang, Q.M., Du, X.F., Wan, L.W., et al., 2020. Characteristics of Mixed Sedimentary Development and Main Control Factors in Lower Third Member of Shahejie Formation of Southern Slope Belt of Laizhouwan Sag, Bohai Sea. Earth Science, 45(10): 3645-3662 (in Chinese with English abstract).

[21] Wang, W., Yang, X.H., Ye, J.R., et al., 2017. Characterization and Controlling Factors of Sandstone Reservoirs in Basalt Development Areas: An Example from Comparison between BZ34-9 and KL6 Areas in South of Bohai Sea. Earth Science, 42(4): 570-586 (in Chinese with English abstract).

[22] Wang, Y.L., Chen, X.G., Hao, Z.W., et al., 2014. Fracture Prediction of Carboniferous Igneous Rock in Junggar Basin, Hashanxi Area. Geological Science and Technology Information, 33(2): 66-71 (in Chinese with English abstract).

[23] Xiao, D.Q., Wang, G.Z., Wei, A.J., et al., 2003. Study on Igneous Reservoir Forming Conditions and Distribution in Huanghua Depression. Special Oil & Gas Reservoirs, 10(1):59-61, 115 (in Chinese with English abstract).

[24] Xie, J., Hu, X. A., Li, B. C., et al., 2021. Reservoir Characteristics and Main Controlling Factors of the Mesozoic Volcanic Rocks in the D Oilfield in Southern Gentle Slope Zone of the Laizhouwan Sag. Frontiers in Earth Science, 9: 663401. https://doi.org/10.3389/feart.2021.663401

[25] Xue, D.J., Hu, B., Fang, Y., et al., 2016. Paths on Thermal-Mechanical Coupling Failure of Meso-Scale Granite under Uni-Axial Tension and Numerical Analysis. Journal of China Coal Society, 41(9): 2212-2222 (in Chinese with English abstract).

[26] Yang, H.F., Lü, D.Y., Sun, Y.H., et al., 2021. The Fault System and Its Tectonophysics Simulation in the Eastern Huanghekou Sag in Bohai Bay Basin. Earth Science, 46(7):2391-2402 (in Chinese with English abstract).

[27] Yang, H.F., Xu, C.G., Niu, C.M., et al., 2018. Joint Control of the Cenozoic Faults and Magma on Hydrocarbon Accumulation in BZ34-9 Oilfield in Huanghekou Sag, Bohai Bay Basin. Oil & Gas Geology, 39(5):1056-1064 (in Chinese with English abstract).

[28] Ye, T., Chen, A. Q., Hou, M. C., et al., 2022. Lithofacies Characteristics and Controlling on Volcanic Reservoirs in the Basement: A Case Study of the Offshore Bohai Bay Basin, Eastern China. Journal of Petroleum Science and Engineering, 209: 109860. https://doi.org/10.1016/j.petrol.2021.109860

[29] Yin, Z.J., Li, C., Lyu, W.J., et al., 2020. Lithofacies Characteristics and Distribution Prediction of Paleogene Volcanic Rocks in BZ34-9 Oilfield. China Offshore Oil and Gas, 32(5):63-72 (in Chinese with English abstract).

[30] Yu, Y. J., 2015. Study on Damage and Fracture Mechanics Behavior of Oil Shale under Thermal Mechanical Coupling (Dissertation). Taiyuan University of Technology, Shanxi (in Chinese with English abstract).

[31] Zhang, Z.H., Du, S.K., Chen, H.Y., et al., 2018. Quantitative Characterization of Volcanic Fracture Distribution Based on Electrical Imaging Logging: A Case Study of Carboniferous in Dixi Area, Junggar Basin. Acta Petrolei Sinica, 39(10):1130-1140 (in Chinese with English abstract).

[32] Zhu, H.T., Liu, Y.M., Wang, Y.L., et al., 2014.Volcanic Eruption Phases and 3-D Characterization of Volcanic Rocks in BZ34-9 Block of Huanghekou Sag, Bohai Bay Basin. Earth Science, 39(9): 1309-1316 (in Chinese with English abstract).

[33] 崔宝文, 蒙启安, 白雪峰, 等, 2018. 松辽盆地北部石油勘探进展与建议. 大庆石油地质与开发, 37(3): 1-9.

[34] 戴世立, 姜传金, 张元高, 2012. 徐家围子断陷安达地区火山岩有效储层预测方法. 地质科技情报, 31(2):67-71.

[35] 冯建伟, 戴俊生, 马占荣, 等, 2011. 低渗透砂岩裂缝参数与应力场关系理论模型. 石油学报, 32(4):664-671.

[36] 高宇慧, 李思远, 巴忠臣, 等, 2021. 准噶尔盆地红172井区石炭系火山岩储层裂缝特征及定量预测. 大庆石油地质与开发, 40(2): 30-41.

[37] 何家雄, 祝有海, 翁荣南, 等, 2010.南海北部边缘盆地泥底辟及泥火山特征及其与油气运聚关系. 地球科学, 35(1): 75-86.

[38] 胡志伟, 杨海风, 韩自军, 等, 2019. 渤海海域渤中34-9油田火山岩特征及其对油气成藏的控制作用. 地质科技情报, 38(3): 208-218.

[39] 季宗镇, 戴俊生, 汪必峰, 2010.地应力与构造裂缝参数间的定量关系. 石油学报, 31(1): 68-72.

[40] 靳是琴, 李殿超, 李宪洲, 等, 1987. 胶东莱西‒平度一带麻粒岩相岩石的变质作用特征. 地质学报, 61(3): 240-252.

[41] 李玉玺,2019.准噶尔盆地石西油田石炭系火山岩裂缝特征及预测（硕士学位论文）.成都：西南石油大学.

[42] 李镇村, 石耀霖, 张泽阳, 等, 2016. 熔岩管冷却过程中的热应力计算. 中国科学院大学学报, 33(3): 412-420.

[43] 林向洋, 苏玉平, 郑建平, 等, 2011. 准噶尔盆地克拉美丽气田复杂火山岩储层特征及控制因素. 地质科技情报, 30(6):28-37.

[44] 刘晓健, 王清斌, 代黎明, 等, 2020. 莱州湾凹陷火山岩‒湖相碳酸盐岩混积岩储层特征和成因. 地球科学, 45(10): 3579-3588.

[45] 任涛, 秦军, 周阳, 等, 2021. 利用叠前方位各向异性预测准噶尔盆地西北缘车476井区火山岩裂缝. 石油地质与工程, 35(1):1-7.

[46] 史伟新,2018.流体饱和岩石热开裂实验研究（硕士学位论文）.青岛：中国石油大学.

[47] 滕吉文, 张中杰, 张秉铭, 等, 1997. 渤海地球物理场与深部潜在地幔热柱的异常构造背景. 地球物理学报, 40(4): 468-480.

[48] 王启明, 杜晓峰, 宛良伟, 等, 2020. 渤海海域莱南斜坡带沙三下亚段混合沉积发育特征及主控因素. 地球科学, 45(10): 3645-3662.

[49] 王维, 杨香华, 叶加仁, 等, 2017. 玄武岩发育区砂岩储层特征及控制因素: 渤海南部BZ34-9区与KL6区块对比. 地球科学, 42(4): 570-586.

[50] 王月蕾, 陈学国, 郝志伟, 等, 2014. 准噶尔盆地哈山西地区石炭系火成岩裂缝预测. 地质科技情报, 33(2): 66-71.

[51] 肖敦清, 王桂芝, 韦阿娟, 等, 2003. 黄骅坳陷火成岩成藏特征研究. 特种油气藏, 10(1):59-61, 115.

[52] 薛东杰, 胡本, 方园, 等, 2016. 热力耦合路径下花岗岩细观破裂特征及数值分析. 煤炭学报, 41(9): 2212-2222.

[53] 杨海风, 吕丁友, 孙永河, 等, 2021. 渤海湾盆地黄河口凹陷东洼断裂体系发育特征及其变形过程的构造物理模拟. 地球科学, 46(7):2391-2402.

[54] 杨海风, 徐长贵, 牛成民, 等, 2018. 渤海湾盆地黄河口凹陷BZ34-9油田新生界“断裂‒岩浆”联合控藏作用. 石油与天然气地质, 39(5):1056-1064.

[55] 尹志军, 李超, 吕文杰, 等, 2020. 渤中34-9油田古近系火山岩岩相特征与分布预测. 中国海上油气, 32(5):63-72.

[56] 于永军,2015.热力耦合作用下油页岩损伤与断裂力学行为研究（硕士学位论文）.山西：太原理工大学.

[57] 张兆辉, 杜社宽, 陈华勇, 等, 2018. 基于电成像测井的火山岩裂缝分布定量表征: 以准噶尔盆地滴西地区石炭系为例. 石油学报, 39(10):1130-1140.

[58] 朱红涛, 刘依梦, 王永利, 等, 2014.渤海湾盆地黄河口凹陷 BZ34-9区带火山岩三维刻画及火山喷发期次. 地球科学, 39(9): 1309-1316.

基金资助

中海石油（中国）有限公司天津分公司科研项目(CCL2020TJX0NST1300)

中国海洋石油有限公司“十四五”重大科技项目(KJGG2022-0405)