基于电-力-声多物理场耦合数值模型的含水合物多孔介质声速和衰减特性研究*

徐 源，张欢欢，邢兰昌†，魏 伟，韩维峰，魏周拓，邓少贵

基于电-力-声多物理场耦合数值模型的含水合物多孔介质声速和衰减特性研究*

徐 源1，张欢欢1，邢兰昌1†，魏 伟2，韩维峰2，魏周拓3，邓少贵3

（1. 中国石油大学（华东）控制科学与工程学院，山东 青岛 266580；2. 中国石油勘探开发研究院新能源研究所，河北 廊坊 065007；3. 中国石油大学（华东）地球科学与技术学院，山东 青岛 266580）

在物理模拟实验中对水合物微观赋存模式和饱和度进行准确控制和评价尚存在技术困难，仅依赖实验技术研究含水合物沉积物声学特性、建立储层参数解释模型存在局限性。采用基于有限元的数字岩石物理技术，针对悬浮、接触和胶结三种典型的水合物微观赋存模式分别建立多孔介质的三维电−力−声多物理场耦合模型，考察了微观赋存模式和水合物饱和度对多孔介质声速和衰减的影响规律，对比了声速数值模拟与理论模型计算结果，建立了声波衰减参数与水合物饱和度之间的关系式。研究结果表明：（1）对于三种水合物赋存模式，由于水合物相比孔隙水具有更高的弹性模量，多孔介质的声速随着水合物饱和度的增大而增大；水合物的存在导致声波在传播过程中遇到更多不连续的声阻抗界面，声衰减随着水合物饱和度的增大而近似线性增大；（2）悬浮和接触赋存模式条件下，水合物饱和度对多孔介质的声速和衰减影响规律基本一致；对于相同的水合物饱和度，胶结模式条件下含水合物多孔介质具有更高的声速和更小的声衰减；（3）通过合理选择参数值，利用权重方程与Lee改进的Biot-Gassmann Theory（BGTL）模型估算的含悬浮和接触模式水合物多孔介质的声速较为准确；通过等效介质理论模型C计算的含胶结模式水合物多孔介质的声速更为准确。研究结果可为获取复杂条件下含水合物沉积物的声学特性提供数值建模方法，为基于声波测井数据的水合物储层精细评价提供理论支撑。

天然气水合物；数字岩石物理；微观赋存模式；水合物饱和度；声速；声衰减

0 引 言

天然气水合物是由水和天然气分子在低温高压环境下结合形成的与冰相似的化合物，其广泛分布于海洋沉积物和冻土区域中[1-4]，是一种极具开发前景且储量巨大的新型清洁能源[5-7]。地球物理测井是准确识别和评价水合物储层的重要手段，其中声波测井是评价水合物饱和度（h）的重要方法之一[8-9]。研究含水合物沉积物的声学响应特性，深入探究水合物饱和度、微观赋存模式等因素对沉积物声波传播速度及其能量衰减的影响规律，将为声波测井数据解释、水合物储层识别与评价提供理论基础。

天然气水合物在储层中具有多种微观赋存模式，不同的微观赋存模式会对沉积物的声学、电学、力学、热学等物性参数产生不同程度的影响[10-13]。为了准确评价储层中水合物的饱和度，在研究水合物饱和度与沉积物中声速和衰减之间的定量关系时，必须考虑水合物微观赋存模式对沉积物声速及衰减的影响规律。国内外有许多学者采用室内模拟实验方法研究了含水合物（模拟）沉积物的声学特性[14-20]。文献[14-17]研究了水合物的存在对不同类型沉积物（固结、松散）声速的影响。PRIEST等[18-19]研究了不同体系（过量气、定量气、溶解气）中水合物的生成对沉积物声速的影响规律。SAHOO等[21]和LIU等[22]分别研究了声波频率和水合物饱和度对沉积物声波衰减特性的影响规律。虽然研究者提出了制备具有不同水合物饱和度、不同微观分布模式的含水合物沉积物样品的实验方法，但是对于所测试的整个样品而言，实验中所得到的饱和度、微观分布模式是否准确、是否具有重复性、是否方便控制，仍然是值得进一步探讨的问题。而数值模拟方法能够准确、方便地控制水合物饱和度、微观赋存模式等影响因素，同时数值模拟方法具有重复性好、时间成本低的显著优势。

数值模拟方法在岩石物理领域具有广泛的应用，逐渐形成了数字岩石物理技术。利用数字岩石物理技术不仅能够实现岩石孔隙尺度建模，也能够研究微观因素对岩石声速、电阻率、渗透率等物理参数的影响[23]。针对岩石声学特性的数值模拟方法主要包括旋转交错网格有限差分方法、声格子方法和有限元方法。SAENGER等[24]采用旋转交错网格有限差分方法模拟出声波在各向异性、黏弹性介质中的传播过程。DEL VALLE-GARCÍA等[25]利用声格子方法研究了声波在多孔介质中的传播过程和规律。MAKARYNSKA等[26]利用有限元方法研究了不同含水饱和度岩石的声学特性。ZHU等[27]利用有限差分方法实现了油气储层介质弹性波场的数值模拟。孙建孟等[28]、刘向君等[29]采用有限元方法计算了数字岩心的弹性参数。YANG等[30]采用旋转交错网格有限差分方法模拟了地震波在海底天然气水合物沉积地层中的传播。卜庆涛等[31]基于实验条件对水合物模拟储层开展了声学特性数值模拟，发现水合物储层声速的实验与数值模拟结果对应关系良好。目前，基于数字岩石物理技术研究微观赋存模式、水合物饱和度等因素对含水合物沉积物声学特性的影响规律鲜有报道。

本文采用基于有限元的数字岩石物理技术，针对悬浮、接触和胶结三种水合物微观赋存模式分别建立多孔介质的三维电−力−声多物理场耦合模型，研究了微观赋存模式和水合物饱和度对多孔介质声速和衰减的影响规律。研究结果可为含水合物多孔介质声学特性数值模拟提供建模方法，为声波测井在水合物储层精细评价中的应用提供理论依据和模型支撑。

1 含水合物多孔介质多物理场有限元建模

1.1 几何结构与材料

采用多物理场耦合计算平台（COMSOL Multiphysics）来建立含水合物多孔介质的有限元数值模型。数值模型的几何结构包含压电晶片和多孔介质两部分。如图1所示，两块压电晶片位于多孔介质两侧，分别用于模拟发射和接收声波。两块压电晶片的几何尺寸和材料相同，晶片厚度为1 mm，边长为1.624 mm。水饱和多孔介质为正方体，边长为1.624 mm。几何结构包含64个大小相同的球形石英砂颗粒，半径为0.21 mm，石英砂颗粒之间规则相切排布。多孔介质的孔隙度为相互连通的孔隙体积与总体积的比值，本文所构建的多孔介质几何结构的孔隙体积约为1.8 mm3，孔隙度为0.42。

依据ECKER[32]提出的悬浮、接触和胶结三种典型的水合物微观赋存模式（如图2所示），构建含水合物多孔介质模型的三维几何结构。图3所示为三种水合物微观赋存模式条件下的含水合物多孔介质的三维立体结构。含悬浮模式水合物多孔介质截面中，代表水合物的9个球体均匀分布于孔隙空间中，水合物与石英砂颗粒不接触。含接触模式水合物多孔介质截面中，代表水合物的9个球体与石英砂颗粒相接触（相交）。含胶结模式水合物多孔介质截面中，代表水合物的薄层包裹在石英砂颗粒外表面，相邻石英砂颗粒表面的水合物层相互接触（相交），将石英砂颗粒胶结在一起。

图1 水饱和多孔介质数值模型几何结构

图2 沉积物中水合物的三种微观赋存模式（改自文献[32]）

图3 含不同微观赋存模式水合物的多孔介质模型三维立体结构：（a）悬浮模式；（b）接触模式；（c）胶结模式

为得到不同的水合物饱和度，对于悬浮和接触模式，保持代表水合物的球体个数不变，改变球体的直径；对于胶结模式，改变石英砂颗粒表面水合物层的厚度。受到三维模型几何结构的限制（保持水合物为单一球形颗粒），悬浮模式条件下模型所模拟的水合物饱和度存在上限，当水合物与石英砂相切时水合物饱和度为23% ～ 24%。考虑到地层中水合物饱和度过高时，水合物的赋存模式倾向于其他模式（如接触模式），因此选取含2%、4%、6%、8%、10%、12%、14%、16%、18%、20%、22%等11个水合物饱和度的多孔介质进行数值求解。悬浮与接触模式条件下，以上11种水合物饱和度对应的水合物颗粒半径的范围为68.47 ～ 152.14 μm。胶结模式条件下对应的水合物层厚度范围为1.16 ～ 13.18 μm。

在分别定义石英砂、孔隙水、水合物和压电晶片等材料的基础上，根据表1设定多孔介质中各组分的声速和密度等参数。压电晶片材料为锆钛酸铅（PZT-5H），密度为7 500 kg/m3。

表1 多孔介质各组分声学特性参数[33]

1.2 物理场控制方程与边界条件

采用静电、固体力学、压力声学多物理场耦合方法实现数值建模方案。采用声−结构边界耦合压力声学和固体力学两个物理场，基于压电效应实现静电和固体力学两个物理场的耦合。

静电物理场中控制方程为高斯方程，如式（1）所示。

式中：为电位移，C/m2；为电场强度，V/m；为电势，V；v为自由电荷密度，C/m3。

固体力学物理场控制方程为运动方程和压电材料本构方程，计算得到位移、应力和应变值，方程如式（2）所示。

式中：为压电晶片材料的密度，kg/m3；为位移场，m；V为体积力，N/m3；tt为压电晶片振动加速度，m/s2；为应力，Pa；为应变，无量纲；E与分别为压电晶片材料的弹性矩阵（Pa）和耦合矩阵（C/m2）；0为真空介电常数，8.854 × 10−12F/m；S为压电晶片材料的相对介电常数，无量纲。

压力声学物理场控制方程为时域波动方程，如式（3）所示。

式中：0为声波传播介质的密度，kg/m3；为介质的声速，m/s；t为总声压，Pa。

以上三个物理场的边界条件设定如下：①对于静电物理场，将用作发射声波的压电晶片内端面定义为电势边界，电势为矩形波调制谐波后的电压脉冲，脉冲幅值0为100 V，频率为；将用作接收声波的压电晶片内端面定义为悬浮电位边界（floating potential）；将两块压电晶片的外端面定义为接地边界。②对于固体力学物理场，将位移场初始值设为0，将两块压电晶片外端面定义为辊支承边界，防止其沿厚度方向移动。③对于压力声学物理场，为保证声波在抵达多孔介质几何外边界时沿法向产生的反射最少，将多孔介质的外边界定义为平面辐射波边界。

声−结构（本模型中的多孔介质−压电晶片）边界条件如式（4）所示。

式中：为边界上的单位法向量，无量纲；A为边界上单位面积所受的力，Pa。

1.3 网格划分与求解器

使用自由四面体网格对压电晶片几何结构进行网格划分。采用三种数量级别的网格、网格尺寸及数量如表2所示。利用稳态（静态）求解器对以上描述静电、固体力学物理场的微分方程进行求解。

表2 压电晶片部分网格尺寸及数量

使用自由四面体网格对多孔介质几何结构进行网格划分。采用多种数量级别的网格、网格尺寸及数量如表3所示。使用与瞬态求解器（广义方法）求解上述静电、固体力学、压力声学物理场中的微分方程，选择并行稀疏直接求解方法（multi-frontal massively parallel sparse direct solver, MUMPS）。在获取多孔介质声速的计算模型中，将求解时长设定为2 μs，时间步长为1 ns；在获取多孔介质声衰减的计算模型中，将求解时长设定为20 μs，时间步长为10 ns。使用组合参数扫描方法对表3中不同数量级别的网格依次进行计算。

表3 多孔介质部分网格尺寸及数量

1.4 网格选取和首波到达时刻提取方法

使用表2所示三种数量级别的网格测试了压电晶片中心频率的模型计算值对网格的依赖性，根据测试结果选取合适的网格尺寸进行后续研究。根据计算结果得到三种网格下压电晶片振动位移随激励电信号频率的变化关系，最大振动位移对应的频率即为压电晶片中心频率。结果表明，压电晶片的中心频率在三种数量级别网格下均为642.5 kHz，由于网格数量成倍增加，可认为基于这三种网格所得到的结果均为正确值。为减少计算时间，针对压电晶片选择单元数量为5724的网格级别，即最大网格尺寸为200 μm，最小网格尺寸为40 μm。

图4 不同网格条件下接收的声波

将激励信号频率设定为642.5 kHz，使用表3中不同级别网格，计算了水饱和多孔介质接收声波的首个波谷到达时刻，测试了模型计算值对网格的依赖性。图4所示为接收声波信号，信号首个波谷到达时刻随网格数量变化情况如图5所示。结果表明，当网格数量达到588 993后，首个波谷到达时刻的模型计算值已经不随网格数量的增加而变化。为了平衡计算精度和计算时间，选择单元数量为588 993的网格级别，对应的最大和最小网格尺寸分别为60 μm和20 μm。

图5 首个声波波谷到达时刻随网格数量的变化

从接收声波信号中准确提取首波到达时刻是获取多孔介质声速的前提。将整个多孔介质几何结构设定为同种材料（即赋予相同的声速值），结合多孔介质几何尺寸能够计算出首波到达时刻理论值。为了建立首波到达时刻提取方法，首先赋予多孔介质中各组分相同的声速和密度值，然后保持密度不变，在1 500 ～ 4 500 m/s范围内设定一系列材料声速值进行数值求解。经过数值求解得到样品的接收声波信号如图6所示，首先从接收信号中提取波形的首个波谷到达时刻，然后计算其与首波到达时刻理论值0的差值Δ=–0。图7所示为Δ随的变化情况。

可以将图7所示Δ随的变化趋势划分为三个区间，两个转折点对应的首个波谷到达时刻分别为1.256 μs和1.457 μs。针对三个区间内Δ随的变化规律，分别采用线性函数进行拟合，拟合结果如式（5）所示。

采用水饱和多孔介质时，接收声波的首个波谷到达时刻t的模型计算值为1.515 μs，结合上述线性拟合式［式（5）］，可以得到水饱和多孔介质的Δt为0.604 μs，因此其首波到达时刻的真实值t0为0.911 μs。已知多孔介质样品沿声波传播方向的长度L为1.624 mm，因此能够计算出水饱和多孔介质的声速Vp为1 782.7 m/s。在后续的模型求解数据分析中，采用相同的方法对所获取的首波到达时刻t进行校正，从而得到真实的首波时刻t0和多孔介质声速。

图7 Δt随首个波谷到达时刻变化

2 含水合物多孔介质声速特性

图8所示为悬浮模式、接触模式和胶结模式条件下含水合物多孔介质接收的声波信号。利用所接收的声波信号计算出首个波谷到达时刻，从而得到三种微观赋存模式下不同水合物饱和度的多孔介质样品的声速。含水合物多孔介质样品的声速随水合物饱和度的变化情况如图9所示。

分析图8和图9可知：

（1）三种赋存模式对应的含水合物多孔介质的声速均随水合物饱和度增大而增大。随着水合物饱和度的增大，更多的孔隙水转化为水合物，由于水合物比孔隙水具有更高的体积模量，因此多孔介质整体的体积模量和剪切模量增大，从而使多孔介质具有更高的声速。

（2）含悬浮和接触模式水合物多孔介质的声速随水合物饱和度的增大而增大的规律基本一致，声速增长速率稳定，与DUCHKOV等[17]采用模拟实验方法得到的声速与饱和度之间的关系一致。当水合物饱和度为22%时，两种赋存模式下多孔介质的声速分别达到2 090.1 m/s和2 117.3 m/s。参见图2和图3，接触模式可通过将悬浮模式下的水合物颗粒运移到与颗粒接触而得到，因此水合物饱和度变化所引起的多孔介质弹性模量的变化差异不大。

（3）含胶结模式水合物多孔介质中水合物饱和度由0增大到2%时，声速出现大幅增长现象，由1 782.7 m/s增长为1 956.6 m/s；水合物饱和度由2%增大到22%时，声速随水合物饱和度的增长速率趋于稳定，饱和度为22%时声速达到2 343.4 m/s。在相同水合物饱和度下，含胶结模式水合物多孔介质的声速显著高于含悬浮模式和接触模式水合物多孔介质[16]。胶结模式条件下的水合物使得多孔介质的体积模量和剪切模量大幅增长，因此多孔介质声速大幅提升。随着水合物饱和度进一步增大，水合物饱和度变化成为声速变化的主导因素，因此声速增长速率趋于稳定（与悬浮和接触模式类似）。

图8 接收的声波信号（时间长度为2 μs）：（a）悬浮模式；（b）接触模式；（c）胶结模式

图9 多孔介质声速随水合物饱和度变化

采用时间平均方程[34]、伍德方程[35]、权重方程[36]、Lee改进的Biot-Gassmann（Biot-Gassmann Theory modified by Lee, BGTL）模型[37]和等效介质理论模型[32,38]等对多孔介质的声速进行计算，计算时采用的介质声速、密度、剪切模量和体积模量等参数参考表1。将三种微观赋存模式对应的含水合物多孔介质声速数值模拟结果与理论模型计算结果进行了对比，如图10所示。其中基于等效介质理论模型A、B和C的声速计算结果分别与含悬浮、接触、胶结模式水合物多孔介质数值模拟结果进行了比较。

图10 声速数值模拟与理论模型计算结果对比：（a）悬浮模式；（b）接触模式；（c）胶结模式

分析图10可知：

（1）对于悬浮模式（图10a），声速数值模拟结果低于时间平均方程计算结果，高于伍德方程计算结果；权重方程（取1.9、取0.8）计算出的声速与数值模拟结果具有较高的一致性；BGTL模型（取0.9、取1.4）声速计算结果与数值模拟结果基本一致；等效介质理论模型A计算出的声速随水合物饱和度增长速率略低于数值模拟结果。

（2）对于接触模式（图10b），声速数值模拟结果低于时间平均方程计算结果，高于伍德方程计算结果；权重方程（取1.9、取1）计算出的声速与数值模拟结果具有较高的一致性；BGTL模型（取1、取1.5）计算结果与数值模拟结果一致性较好；等效介质理论模型B计算出的声速随水合物饱和度增长速率略低于数值模拟结果。

（3）对于胶结模式（图10c），声速数值模拟结果低于时间平均方程计算结果，高于伍德方程计算结果；权重方程（取1.7、取2.9）计算出的声速与数值模拟结果差异最小；BGTL模型（取1、取0.6）计算结果在水合物饱和度大于10%时与数值模拟结果几乎一致，但在小于10%饱和度区间内，BGTL模型计算结果显著高于数值模拟结果；等效介质理论模型C计算结果与数值模拟结果高度一致。

（4）对于三种赋存模式对应的含水合物多孔介质，在所研究的水合物饱和度区间内，基于时间平均方程的多孔介质声速预测值偏高，而伍德方程预测值偏低[39]；权重方程和BGTL模型对预测含悬浮和接触模式水合物多孔介质声速具有优势；等效介质理论模型C对预测含胶结模式多孔介质声速具有独特优势。

3 含水合物多孔介质声衰减特性

声波在多孔介质中传播时，在声阻抗不连续的界面处会产生折射与反射现象，从而导致声波能量衰减。本节讨论多孔介质中水合物饱和度和赋存模式对声衰减的影响。经过数值求解，得到针对含悬浮模式、接触模式和胶结模式水合物多孔介质所接收到的声波信号，如图11所示。研究声波信号的衰减时，将信号的振幅（以最大电压表示）作为考察指标。基于信号振幅的降低定义了用于量化声波能量衰减的参数a，a由公式（6）计算得到。

从图11所示接收声波信号中，提取三种赋存模式对应的含水合物多孔介质在不同水合物饱和度条件下的信号振幅，计算得到多孔介质声波衰减参数Ma。图12展示了含水合物多孔介质声波衰减参数Ma随水合物饱和度变化的情况。

图12 声波衰减参数Ma随水合物饱和度变化

分析图12可知：

（1）三种水合物赋存模式对应的含水合物多孔介质声衰减参数a随着水合物饱和度的增大都呈现出近似线性增长的规律[22]。与水饱和多孔介质相比较，水合物的存在导致声波在传播过程中遇到更多不连续的声阻抗界面，从而加剧声波能量的衰减。采用线性函数对两者关系进行拟合，得到式（7）。

（2）含悬浮和接触模式水合物多孔介质的声衰减参数a随水合物饱和度的变化规律基本一致，这与以上声速随水合物饱和度变化规律类似，分析其原因为悬浮和接触模式条件下的水合物对多孔介质内部声阻抗结构影响的差异性不显著。

（3）在相同的水合物饱和度条件下，含胶结模式水合物多孔介质的声衰减参数a显著低于悬浮和接触模式。与悬浮和接触模式相比较，胶结模式条件下多孔介质内部组分界面出现差异。主要差异在于：前两种模式中，水直接与石英砂接触，而在胶结模式中水合物与石英砂接触，同时水与石英砂不接触。考虑到水合物与石英砂之间的声阻抗差异小于水与石英砂，因此胶结模式条件下多孔介质声衰减参数a较低，即衰减较小。

4 结 论

基于电−力−声多物理场耦合建模方法，针对悬浮、接触和胶结三种水合物微观赋存模式，分别建立了含水合物多孔介质的有限元数值模型，分析了微观赋存模式和水合物饱和度对多孔介质声速及衰减的影响规律，对比了声速数值模拟与理论模型计算结果，建立了声波衰减参数与水合物饱和度之间的关系式。得到如下结论：

（1）对于三种水合物赋存模式，由于水合物相比孔隙水具有更高的弹性模量，多孔介质的声速随着水合物饱和度的增大而增大；水合物的存在导致声波在传播过程中遇到更多不连续的声阻抗界面，声衰减随着水合物饱和度的增大而近似线性增大。

（2）悬浮和接触赋存模式条件下，水合物饱和度对多孔介质的声速和衰减影响规律基本一致；对于相同的水合物饱和度，胶结模式条件下含水合物多孔介质具有更高的声速和更小的声衰减。

（3）通过合理选择参数值，利用权重方程与BGTL模型估算的含悬浮和接触模式水合物多孔介质的声速较为准确；通过等效介质理论模型C计算的含胶结模式水合物多孔介质的声速更为准确。

[1] LI X S, XU C G, ZHANG Y, et al. Investigation into gas production from natural gas hydrate: a review[J]. Applied energy, 2016, 172: 286-322. DOI: 10.1016/j.apenergy. 2016.03.101.

[2] YE J L, QIN X W, XIE W W, et al. The second natural gas hydrate production test in the South China Sea[J]. China geology, 2020, 3(2): 197-209. DOI: 10.31035/ cg2020043.

[3] LIU J W, LI X S. Recent advances on natural gas hydrate exploration and development in the South China Sea[J]. Energy & fuels, 2021, 35(9): 7528-7552. DOI: 10.1021/ acs.energyfuels.1c00494.

[4] 苏正, 刘丽华. 南海北部陆坡天然气水合物成藏特征研究进展分析[J]. 新能源进展, 2020, 8(1): 35-41. DOI: 10.3969/j.issn.2095-560X.2020.01.006.

[5] 周守为, 李清平, 吕鑫, 等. 天然气水合物开发研究方向的思考与建议[J]. 中国海上油气, 2019, 31(4): 1-8. DOI: 10.11935/j.issn.1673-1506.2019.04.001.

[6] LI J F, YE J L, QIN X W, et al. The first offshore natural gas hydrate production test in South China Sea[J]. China geology, 2018, 1(1): 5-16. DOI: 10.31035/cg2018003.

[7] KOH D Y, KANG H, LEE J W, et al. Energy-efficient natural gas hydrate production using gas exchange[J]. Applied energy, 2016, 162: 114-130. DOI: 10.1016/j. apenergy.2015.10.082.

[8] 钟广法, 张迪, 赵峦啸. 大洋钻探天然气水合物储层测井评价研究进展[J]. 天然气工业, 2020, 40(8): 25-44. DOI: 10.3787/j.issn.1000-0976.2020.08.002.

[9] SHANKAR U, RIEDEL M. Gas hydrate saturation in the Krishna-Godavari basin from P-wave velocity and electrical resistivity logs[J]. Marine and petroleum geology, 2011, 28(10): 1768-1778. DOI: 10.1016/j. marpetgeo.2010.09.008.

[10] KOU X, LI X S, WANG Y, et al. Effects of gas occurrence pattern on distribution and morphology characteristics of gas hydrates in porous media[J]. Energy, 2021, 226: 120401. DOI: 10.1016/j.energy.2021.120401.

[11] DAI S, SANTAMARINA J C, WAITE W F, et al. Hydrate morphology: physical properties of sands with patchy hydrate saturation[J]. Journal of geophysical research: solid earth, 2012, 117(B11): B11205. DOI: 10.1029/2012JB009667.

[12] 宁伏龙, 梁金强, 吴能友, 等. 中国天然气水合物赋存特征[J]. 天然气工业, 2020, 40(8): 1-24. DOI: 10.3787/ j.issn.1000-0976.2020.08.001.

[13] XING L C, QI S Y, XU Y, et al. Numerical study on complex conductivity characteristics of hydrate-bearing porous media[J]. Journal of natural gas science and engineering, 2021, 95: 104145. DOI: 10.1016/j.jngse. 2021.104145.

[14] WINTERS W J, WAITE W F, MASON D H, et al. Methane gas hydrate effect on sediment acoustic and strength properties[J]. Journal of petroleum science and engineering, 2007, 56(1/3): 127-135. DOI: 10.1016/j.petrol. 2006.02.003.

[15] YUN T S, FRANCISCA F M, SANTAMARINA J C, et al. Compressional and shear wave velocities in uncemented sediment containing gas hydrate[J]. Geophysical research letters, 2005, 32(10): L10609. DOI: 10.1029/2005GL022607.

[16] HU G W, YE Y G, ZHANG J, et al. Acoustic properties of hydrate-bearing unconsolidated sediments measured by the bender element technique[J]. Chinese journal of geophysics, 2012, 55(6): 635-647. DOI: 10.1002/cjg2. 1758.

[17] DUCHKOV A D, DUCHKOV A A, PERMYAKOV M E, et al. Acoustic properties of hydrate-bearing sand samples: laboratory measurements (setup, methods, and results)[J]. Russian geology geophysics, 2017, 58(6): 727-737. DOI: 10.1016/j.rgg.2016.09.029.

[18] PRIEST J A, BEST A I, CLAYTON C R I. A laboratory investigation into the seismic velocities of methane gas hydrate-bearing sand[J]. Journal of geophysical research: solid earth, 2005, 110(B4): B04102. DOI: 10.1029/ 2004JB003259.

[19] PRIEST J A, REES E V L, CLAYTON C R I. Influence of gas hydrate morphology on the seismic velocities of sands[J]. Journal of geophysical research: solid earth, 2009, 114(B11): B11205. DOI: 10.1029/2009JB006284.

[20] REN S R, LIU Y J, LIU Y X, et al. Acoustic velocity and electrical resistance of hydrate bearing sediments[J]. Journal of petroleum science and engineering, 2010, 70(1/2): 52-56. DOI: 10.1016/j.petrol.2009.09.001.

[21] SAHOO S K, NORTH L J, MARÍN-MORENO H, et al. Laboratory observations of frequency-dependent ultrasonic P-wave velocity and attenuation during methane hydrate formation in Berea sandstone[J]. Geophysical journal international, 2019, 219(1): 713-723. DOI: 10.1093/gji/ggz311.

[22] LIU Z C, NING F L, HU G W, et al. Characterization of seismic wave velocity and attenuation and interpretation of tetrahydrofuran hydrate-bearing sand using resonant column testing[J]. Marine and petroleum geology, 2020, 122: 104620. DOI: 10.1016/j.marpetgeo.2020.104620.

[23] 赵建鹏, 陈惠, 李宁, 等. 三维数字岩心技术岩石物理应用研究进展[J]. 地球物理学进展, 2020, 35(3): 1099-1108. DOI: 10.6038/pg2020DD0486.

[24] SAENGER E H, BOHLEN T. Finite-difference modeling of viscoelastic and anisotropic wave propagation using the rotated staggered grid[J]. Geophysics, 2004, 69(2): 583-591. DOI: 10.1190/1.1707078.

[25] DEL VALLE-GARCÍA R, SÁNCHEZ-SESMA F J. Rayleigh waves modeling using an elastic lattice model[J]. Geophysical research letters, 2003, 30(16): 1866. DOI: 10.1029/2003GL017600.

[26] MAKARYNSKA D, GUREVICH B, CIZ R, et al. Finite element modelling of the effective elastic properties of partially saturated rocks[J]. Computers & geosciences, 2008, 34(6): 647-657. DOI: 10.1016/j.cageo.2007.06.009.

[27] ZHU X, MCMECHAN G A. Numerical simulation of seismic responses of poroelastic reservoirs using Biot theory[J]. Geophysics, 1991, 56(3): 328-339. DOI: 10.1190/1.1443047.

[28] 孙建孟, 闫国亮, 姜黎明, 等. 基于数字岩心研究流体性质对裂缝性低渗透储层弹性参数的影响规律[J]. 中国石油大学学报(自然科学版), 2014, 38(3): 39-44. DOI: 10.3969/j.issn.1673-5005.2014.03.006.

[29] 刘向君, 朱洪林, 梁利喜. 基于微CT技术的砂岩数字岩石物理实验[J]. 地球物理学报, 2014, 57(4): 1133- 1140. DOI: 10.6038/cjg20140411.

[30] YANG J J, HE B S, ZHANG J Z. Multicomponent seismic forward modeling of gas hydrates beneath the seafloor[J]. Applied geophysics, 2014, 11(4): 418-428. DOI: 10.1007/s11770-014-0465-x.

[31] 卜庆涛, 刘圣彪, 胡高伟, 等. 含水合物沉积物声学特性——实验模拟与数值模拟的对比分析[J]. 海洋地质前沿, 2020, 36(9): 56-67. DOI: 10.16028/j.1009-2722. 2020.079.

[32] ECKER C. Seismic characterization of methane hydrate structures[M]. Stanford University ProQuest Dissertations Publishing, 1998, 9837195.

[33] LEE M W, COLLETT T S. Pore- and fracture-filling gas hydrate reservoirs in the gulf of Mexico gas hydrate joint industry project leg II green canyon 955 H well[J]. Marine and petroleum geology, 2012, 34(1): 62-71. DOI: 10.1016/j.marpetgeo.2011.08.002.

[34] WYLLIE M R J, GREGORY A R, GARDNER G H F. An experimental investigation of factors affecting elastic wave velocities in porous media[J]. Geophysics, 1958, 23(3): 459-493. DOI: 10.1190/1.1438493.

[35] PEARSON C F, HALLECK P M, MCGUIRE P L, et al. Natural gas hydrate deposits: a review of in situ properties[J]. The journal of physical chemistry, 1983, 87(21): 4180-4185. DOI: 10.1021/j100244a041.

[36] WOOD W T, STOFFA P L, SHIPLEY T H. Quantitative detection of methane hydrate through high-resolution seismic velocity analysis[J]. Journal of geophysical research: solid earth, 1994, 99(B5): 9681-9695. DOI: 10.1029/94jb00238.

[37] LEE M W. Biot–Gassmann theory for velocities of gas hydrate-bearing sediments[J]. Geophysics, 2002, 67(6): 1711-1719. DOI: 10.1190/1.1527072.

[38] HELGERUD M B, DVORKIN J, NUR A, et al. Elastic-wave velocity in marine sediments with gas hydrates: Effective medium modeling[J]. Geophysical research letters, 1999, 26(13): 2021-2024. DOI: 10.1029/1999GL900421.

[39] 雷亚妮, 吴时国, 孙金, 等. 海洋天然气水合物综合地球物理识别方法评述[J]. 中南大学学报(自然科学版), 2022, 53(3): 864-878. DOI: 10.11817/j.issn.1672-7207. 2022.03.009.

Sound Velocity and Attenuation Characteristics of Hydrate-Bearing Porous Media Based on Numerical Model of Electrical-Mechanical-Acoustic Multi-Physics-Field Coupling

XU Yuan1, ZHANG Huan-huan1, XING Lan-chang1, WEI Wei2, HAN Wei-feng2, WEI Zhou-tuo3, DENG Shao-gui3

(1. College of Information and Control Engineering, China University of Petroleum (East China), Qingdao 266580, Shandong, China; 2. Department of Alternative Energy, PetroChina Research Institute of Petroleum Exploration & Development, Langfang 065007, Hebei, China; 3. School of Geosciences, China University of Petroleum (East China), Qingdao 266580, Shandong, China)

There are still technical difficulties in accurately controlling and evaluating the micro-distribution mode and saturation of hydrate in physical simulation experiments. Limitations exist in the experimental technologies for investigating the acoustic characteristics of hydrate-bearing sediments and establishing interpretation models of reservoir parameters. Three-dimensional numerical models were established for hydrate-bearing porous media based on digital rock physics technology. For the three kinds of hydrate micro-distribution modes (suspension, contact, and cementation), finite-element models were established individually based on the method of electrical-mechanical-acoustic multi-physics-field coupling. The effects of micro-distribution mode and hydrate saturation on sound velocity and attenuation of porous media were examined. The results of sound velocity from the numerical and theoretical models were compared. The relationship between acoustic attenuation parameters and hydrate saturation was established. It was demonstrated that: (1) for the three micro-distribution modes, the sound velocity of porous media increased with the increase of hydrate saturation due to the higher elastic modulus of hydrate compared with pore water; the existence of hydrate caused the sound wave to encounter more discontinuous acoustic impedance interfaces in the process of propagation, thus the sound attenuation increased approximately linearly with the increase of hydrate saturation; (2) under the conditions of suspension and contact modes, the influence of hydrate saturation on sound velocity and attenuation in porous media was quite similar; under the condition of cementation mode, the sound velocity was higher and the sound attenuation was lower for the porous media with the same hydrate saturation; (3) based on a reasonable selection of parameters, the sound velocity of hydrate porous media with suspension and contact modes of hydrate estimated by the weighted equation and Biot-Gassmann Theory modified by Lee (BGTL) model was more accurate; the sound velocity of porous media containing cemented hydrate calculated by the effective medium theoretical model C was more accurate. This study may provide a numerical modeling method for obtaining the acoustic characteristics of hydrate-bearing sediments under complex conditions, and provide a theoretical support for the fine evaluation of hydrate reservoirs based on acoustic logging data.

natural gas hydrate; digital rock physics; micro-distribution mode; hydrate saturation; sound velocity; sound attenuation

2095-560X（2022）05-0400-10

TK01；P631；P315

A

10.3969/j.issn.2095-560X.2022.05.002

2022-04-29

2022-07-25

中石油重大科技项目（ZD2019-184-001）；中央高校基本科研业务费专项资金项目（20CX05005A）；中国石油科技创新基金项目（2018D-5007-0214）；山东省自然科学基金项目（ZR2019MEE095）

邢兰昌，E-mail：xinglc@upc.edu.cn

徐 源（1998-），男，硕士研究生，主要从事天然气水合物储层物性参数检测技术研究。

张欢欢（1997-），女，硕士研究生，主要从事天然气水合物储层物性参数检测技术研究。

邢兰昌（1983-），男，博士，副教授，主要从事天然气水合物、CO2等多相流相关检测理论与方法、多物理场耦合数值模拟方法、智能感知与检测技术、计算机测控系统研究工作。