干热岩电学数值模拟方法与微观特性分析

谢关宝

1 页岩油气富集机理与有效开发国家重点实验室， 北京 102206 2 中石化石油工程技术研究院有限公司， 北京 102206

0 引言

干热岩能源是一种高效低碳清洁能源，目前全球发现的干热岩储量是所有油气总储量的数十倍(冉恒谦等，2010；杨方等，2012；陆川和王贵玲，2015；许天福等，2012，2016；付亚荣等，2018)，其内部不存在流体或仅有少量流体(或致密不透水)的高温岩体，储层温度一般在180 ℃以上，裂缝或断层较为发育.作为一种重要的非常规资源(或能源)，研究其高温条件下岩石物理性性质，可指导干热岩热储资源的高效开发；但在一些特殊情况下，如裂缝岩石、高温超高温条件，常规岩石物理实验难以开展，岩石物理响应及特征研究面临极大挑战.

近年来，数字岩心技术得到迅猛发展，通过构建数字岩心模型，分析岩石物理性质，已成为常规岩石物理实验的有益补充.考虑岩石粒度分布，构建多孔隙数字岩心模型，可以精细表征岩石内部微观结构，模拟实际岩石的导电特性(Bakke and Øren，1997； Øren and Bakke，2003).基于X-CT扫描构建碳酸盐岩模型，可分析碳酸盐岩力学特性(Arns et al.，2002，2004，2005a,b).对于页岩及其他致密岩性，可以采用Nano-CT(Sisk et al.，2010；Song et al.，2018)、FIB-SEM(Sok et al.，2009；Curtis，2010；聂昕等，2016；Salvati and Korsunsky，2017)、MAPS(李潮流等，2016，王民等，2018)等多种技术手段，构建纳米数字岩心，分析页岩孔隙结构特征.干热岩热储层位的岩石物理性质，在高温高压环境下会发生变化，尤其是高温引起的电学性质变化，已成为识别干热岩资源的重要特征之一(Spichak et al.，2007；Chen，2015).许多学者研究了电阻率与温度之间的关系(Winsauer et al.，1952；Ucok et al.，1979；Kariya and Shankland，1983；孟银生等，2010；余琪祥等，2011；童小龙等，2020).苏联科学家(Cheremensky，1982)，总结了岩石电阻率与温度之间的经验公式，但忽略了岩石热开裂引起的电阻率异常变化.基于数字岩心分析技术，改变岩石微观参数，可分析岩石的微观物理性质(Adler et al.，1990；Blunt et al.，2002；Yue et al.，2011；Faisal et al.，2017；Tan et al.，2021).基于枫丹白露砂岩孔隙网络模型，不考虑水湿岩石水膜传导作用，在低含水饱和度情况下，电阻率指数的模拟结果远大于实验结果(Knackstedt et al.，2007)；如考虑水膜传导作用，岩石电阻率模拟结果与实验结果准确吻合(Liu et al.，2009).基于格子玻尔兹曼法，可模拟气水分离过程，分析天然气在地层水中的润湿性和溶解度，分析其对储层岩石电性的影响(Jiang et al.，2011).基于有限元法，采用消息传递机制(MPI)并行算法，可模拟超过2000×2000×2000个像素的岩石电性参数(Amabeoku et al.，2013).目前，基于数字岩心的岩石物理特性数值模拟研究面临一个重大挑战，即图像分辨率与岩心尺寸(即图像分辨率与储层非均质性)之间的矛盾(孔强夫等，2015).

我国青海省共和盆地干热岩发育，地热资源丰富.本研究以该地区热储层岩心为研究对象，采用X-CT、QEMSCAN、Maps等多种实验方法，构建了多尺度、多组分、高分辨率数字岩心；采用分形布朗运动方法(fBm)，构建了裂缝性数字岩心；考虑干热岩储层高温的特点，在数值模拟中引入了温度对岩石骨架的影响，并利用有限元法计算了岩石电阻率，分析了干热岩电学特性微观响应机理，为干热岩热储层位评价和地热资源开发提供了微观理论基础.

1 多尺度数字岩心构建

研究所用岩心取自共和盆地GR1井，其热储层位主要岩性为黑云母片麻岩、花岗岩和花岗闪长岩，孔隙度、渗透率都极低，裂缝和断层较为发育(陈梓慧等，2015).研究钻取了三颗岩心，如图1所示，分别编号为No.1、No.2、No.3，对应取心深度为2980 m、3380 m和3600 m，岩性均为花岗岩.采用X射线衍射(XRD)及覆压孔渗获实验取了岩心主要矿物成分与孔隙度等参数，如表1所示.实验结果表明，研究区干热岩主要造岩矿物为石英、斜长石和钾长石，含少量黏土矿物；孔隙度非常小(最小为0.44%)，岩心较为致密，其中φ3>φ2>φ1.为了构建高分辨率数字岩心，采用QEMSCAN进行定量矿物分析，采用Maps进行微观孔隙分析，采用X-CT扫描构建三维数字模型.

首先，对样品进行二维大面积扫描电镜成像(Maps，分辨率为250 nm，图2a)和QEMSCAN矿物成分分析(图2b)，确定岩心孔隙尺寸分布特征及矿物分布情况.根据岩心孔径分布特征，选择合适的分辨率进行X-CT扫描. 然后，根据矿物分布特征来确定纳米级别的扫描区域(图2b中的黄圈)，钻取子样品进行纳米CT扫描，1号样品子样品记为GRY1，2号样品子样品记为GRY2，3号样品子样品记为GRY3-1.最后，由Maps图像提取干热岩样品孔隙半径特征，发现其孔隙尺寸很小，主要分布范围为0.25～4.0 μm，后续采用的X-CT扫描分辨率为500 nm，其三维重构图如图2c所示，但此分辨率下孔隙识别度很低，孔隙度仅0.12%.因此，将X-CT扫描分辨率提高至65 nm，子样品名称为GRY3-2，其三维重构图如图2d所示，此时孔隙较为容易识别，孔隙度为1.39%.

图1 干热岩柱塞样(a) 1号样(25 mm); (b) 2号样(25 mm); (c) 3号样(25 mm).Fig.1 HDR plug cores(a) No.1 25 mm plug sample; (b) No.2 25 mm plug sample; (c) No.3 25 mm plug sample.

表1 干热岩XRD和孔隙度实验分析Table 1 XRD and porosity measurement experiment of HDR

Maps扫描区域与QEMSCAN扫描区域一致，将两者图像重叠对比，可以确定Maps图像中每个像素点对应的灰度、矿物类型和微孔隙发育程度.根据Maps中矿物灰度值大致确定矿物密度排序，颜色越黑代表矿物密度越低，颜色越亮白代表矿物密度越高.通过3个Maps图像分析可得到矿物的密度，密度大小排序为：ρ石英≤ρ斜长石<ρ钾长石<ρ绿泥石+方解石<ρ黄铁矿.这个密度排序为后续X-CT重建三维灰度图像的多阈值分割提供了参考依据.

数字岩心组分划分通常采用单阈值进行分割，0表示孔隙格架，1表示骨架，但单一阈值分割方式会丧失某些矿物成分信息.因此，在X-CT扫描图像的基础上，结合QEMSCAN矿物识别技术和XRD矿物组分分析，依据前述Maps获取的矿物密度排序，对干热岩数字岩心进行多阈值分割.数字岩心矿物组分分割结果显示如图2e、图2f所示.可以看出，干热岩样品矿物以石英、斜长石和钾长石为主，含少量绿泥石、蒙脱石等黏土矿物，含极少量铁矿物，且岩心孔隙度均很小.从表2可以看到，从数字岩心得到的各矿物组分与XRD矿物分析结果接近，因此，可将数字岩心各矿物组分的电学综合响应定义为骨架电阻率，根据理论方法和岩石物理实验，构建骨架电阻率与温度关系.通过多次扫描实验构建的数字核心分别命名为GRS1、GRS2、GRS3-1、GRS3-2(65 nm).

2 电学特性微观数值模拟方法

电学特性的数值模拟通常采用有限元方法(FEM)进行.考虑温度对岩心电阻率的影响，将温度对岩心电学特性的影响，转换为其对岩心骨架的影响；将不同温度下的骨架电导率作为有限元输入，模拟不同温度下整体数字岩心电流，得到整个数字岩心电阻率.

表2 干热岩矿物成分对比分析Table 2 Companion of the mineral composition of HDR from XRD and constructed digital cores

图2 多种扫描实验及多组分高分辨三维数字岩心(a) Maps扫描结果； (b) QEMSCAN扫描结果； (c) 500 nmX-CT扫描结果； (d) 65 nmX-CT扫描结果； (e) 三维多组分模型； (f) 三维孔隙模型.Fig.2 Multiple scanning experiments and multi-component high-resolution 3D digital cores(a) High-resolution scanning by Maps; (b) mineral analysis by QEMSCAN; (c) Reconstruction of 3D grayscale model with 500 nm X-CT; (d) reconstruction of 3D grayscale model with 65 nm X-CT; (e) 3D multi component digital core, including potash (light purple), mixture of quartz and plagioclase (white), chlorite (red), and ferriferous minerals (yellow) ; (f) 3D pore space.

2.1 电阻率与温度关系

干热岩大多属于酸性侵入岩，骨架导电能力很差，可看作绝缘体.当温度T升高时，岩心电阻率ρcore会降低，岩石电阻率和温度有如下关系(Cheremensky，1982)：

(1)

式中，α是经验系数；ρ0初始电阻率(Ωm)；T0为初始温度(℃)；ρcore是岩心电阻率(Ωm)；T为温度(℃).

选取1号、2号、3号岩心子样本测量其不同温度下电阻率，样品分别命名为No.1-1、No.1-2、No.1-3、No.2、No.3-1、No.3-2；测量样品单轴加压(12 MPa)，利用电阻丝对样品加热，在测量温度点稳定15 min；测量温度从室温开始，至500 ℃，每50 ℃测量一次，每个测量点测量5次，取平均值作为该测量点数值，结果如图3所示.电阻率的对数与温度大致成线性相关特征，与式(1)的描述一致.拟合可以得式(1)中系数α(本地区为0.02)，ρ0在150 ℃时为7.0×108Ωm.

图3 不同温度干热岩电阻率Fig.3 The measured HDR resistivity results at different temperatures

干热岩体可视作一种多孔介质复合材料，包括骨架和孔隙两部分.假定孔隙不导电，则岩心电阻率ρrock与骨架电阻率ρskeleton有如下关系式：

(2)

式中，φ为孔隙度(小数)；ρskeleton为骨架电阻率(Ωm)；ρrock为岩石电阻率(Ωm).

通过岩心测量的岩心电阻率，可求解骨架电阻率，将其代入有限元算法中，即可模拟不同温度下干热岩岩心电学参数.

2.2 有限元算法

有限元基本思想是变分原理，即将求解每个像素上的电压分布问题，转化为求解系统线性能量极值问题，从而确定数字岩心有效电导率.为使能量En取极小值，需满足En对变量um(结点电压)的偏导数均为零，即：

(3)

在三维数字岩心中，由于其自身就是由离散的像素构成，选定每一个像素为一个单元，单元的形状为立方体，每个像素通过8个顶点与相邻单元连接，像素的能量由各个节点上的电压确定，则像素内的能量为

(4)

式中，σpq电导率张量；e(ex，ey，ez)是像素内点(x，y，z)处的局部场.

在有限元方法中每个像素的能量可以表示为

(5)

式中，矩阵Drs是三维模型中的8×8矩阵，也称为“刚度矩阵”，ur是像素上第r(1，2,…，8)个节点处的电压；us是像素上第s(1，2,…，8) 个节点处的电压.

将所有像素能量相加，可以计算整个三维数字岩心的总能量，将系统总能量取最小值，计算所有节点电压.如果已知每个像素单元的电导率张量，可以求得每个单元内的电流.然后，对三维数字岩心中所有单元电流取体积平均，得到三维数字岩心的等效电流.最后，根据电流与电压的关系，得到岩样等效电阻率.

3 电学微观响应特征分析

设置模拟温度10～500 ℃，每50 ℃模拟一次，数字岩心代表体积元设置为400体素，模拟结果以Z-Axis为准，计算不同温度下的数字岩心电阻率，模拟结果如图4所示.随温度升高，干热岩电阻率呈减小变化趋势，与式(1)的描述是一致的，且GRS1、GRS2数值模拟结果与实验结果也相吻合；但GRS3-1、GRS3-2的数值模拟结果与实验结果不一致，在其物理实验中间段，变化范围较大且呈现非线性特征，与式(1)相差较大.

是什么原因导致3号岩心电阻率变化趋势与1号、2号不同呢？与1号、2号样品相比，3号样品埋藏深度更深，储层温度更高，热应力超过岩心抗压强度临界值，导致热开裂或岩心破裂，如图5所示.此外，在3号样品的Maps和X-CT图像中可以清楚地看到一些微裂缝.实验室加热和原生地下高温条件的共同作用导致3号样品的裂缝更加发育，甚至破裂.

图4 不同温度下干热岩电阻率数值模拟结果(a) GRS1； (b) GRS2； (c) GRS3-1； (d) GRS3-2.Fig.4 Resistivity numerical simulation results of HDR at different temperatures(a) GRS1； (b) GRS2； (c) GRS3-1； (d) GRS3-2.

图5 3号岩样裂缝显示(a) Maps图像； (b) 65 nm X-CT图像； (c) 岩心开裂照片.Fig.5 Fractures in maps image; X-CT image and sample of No. 3 HDR(a) Maps image; (b) 65 nm X-CT image; (c) Rock sample after high temperature heating.

图6为3号样品的2个子样品不同温度下的电阻率数值模拟结果.由图可见，随着温度升高，其电阻率在300 ℃后发生较大差异变化.图6中虚线指示的是两个子样品的平均电阻率，当温度达到300 ℃时，电阻率变化趋势发生了变化，故推测此时岩样产生了裂缝.因此，加热过程中产生裂缝的温度边界约为300 ℃.

图6 3号样品不同温度下电阻率： 热开裂温度边界约为300 ℃Fig.6 Resistivity measurement results of two subsamples from the No.3 sample at different temperatures, showing the temperature boundary of thermal cracking is about 300 ℃.

为进一步分析热裂解产生裂缝对电阻率的影响，研究采用分形布朗运动(fBm)方法在3号数字岩心中人为添加裂缝(Voss，1988； Kim， 2007)，形成带裂缝数字岩心模型，命名为GRS3-3.因研究区裂缝多为低倾角，因此模拟中添加的裂缝均为水平裂缝.为了模拟加热和热开裂实际情况，模拟中随着温度的升高，生成的裂缝数量也相应增加.具体模拟过程及参数如下：当300 ℃时，添加第一条裂缝，其长度200体素，开度为2体素；当400 ℃时，添加第二条裂缝，450 ℃时添加第三条裂缝；随着温度升高，先前生成裂缝的尺度会增大(假设长度每升高50 ℃增长100体素)；当温度升至450 ℃(或第一条裂缝贯穿岩心)时，则不再增加裂缝条数.整个添加裂缝的过程如图7a—图7c所示，最后构建的裂缝性数字岩心模型如图7d所示.

对裂缝性数字岩心进行电学特性微观数值模拟，结果如图8所示.由图可见，在300 ℃前，电阻率仍呈线性下降趋势；加入裂缝后，由于裂缝内空气不导电，电阻率减少趋势变缓，并保持稳定(在岩心矿物相变前，如矿物脱水).此时，含裂缝数字岩心GRS3-3模拟的电阻率，与岩石实验测得的平均电阻率变化趋势一致，这说明裂缝性数字岩心构建与数值模拟方法是正确的.在获取岩样或进行高温岩石物理实验受限的情况下，利用上述数字岩心方法，能够更方便地分析干热岩高温电学性质.

图7 不同温度增加裂缝模型示意图(a) 350 ℃时裂缝模型； (b) 400 ℃时裂缝模型； (c) 450 ℃时裂缝模型 ； (d) 450 ℃时多裂缝模型.Fig.7 Schematic diagram of adding fractures at different temperatures(a) Adding the first fracture at 350 ℃；(b) Adding the second fracture at 400 ℃；(c) Adding the third fracture at 450 ℃； (d) The combined digital core with multipal fractures at 450 ℃.

图8 3号岩样电响应数值模拟结果Fig.8 Resistivity simulation results of two sub-cores from No.3 sample after adding fractures

4 结论

为分析不同温度下干热岩电学特性与响应机理，构建了多尺度、多组分、高分辨率数字岩心三维模型，分析了多尺度孔隙结构及和矿物成分特征，并研究了不同温度下干热岩电阻率变化规律.具体认识如下：

(1)采用X-CT、Maps、QEMSCAN和XRD多种实验分析方法，构建了三维多尺度、多组分、高分辨率数字岩心模型，分析了干热岩储层孔隙特征和矿物成分特征，为研究干热岩电学特性提供了有效载体.

(2)干热岩电阻率随温度升高呈线性下降趋势.对于深埋层伴有热开裂的岩石，300 ℃前电阻率仍呈线性下降，300 ℃后电阻率变化幅度减小并趋于稳定.

(3)裂缝性数字岩心模拟与高温岩石物理实验结果一致，表明岩心加热过程中可加热过程中可能出现热开裂.

(4)数字岩心技术可快速研究干热岩高温石物理特性，可有效解决现有岩石物理实验无法定量分析的问题，为研究干热岩热储层高温岩石特性提供了一种新方法.

干热岩高温电学特性在300 ℃出现拐点现象，文中利用裂缝性数字岩心进行了验证，物模结果与数模结果有较好的一致性；但温度引起的干热岩“裂纹到缝网再到破裂”这个过程是比较复杂的，还需综合考虑干热岩体的缝网拓展机理、岩心矿物相变、矿物会脱水、原位地层应力条件等因素，才能够对干热岩储层岩石物理性质有更深刻、更准确的认识.