水合物生成导致沉积物孔隙结构和渗透率变化的低场核磁共振观测

张永超，刘昌岭，刘乐乐，陈鹏飞，张准，孟庆国

1.自然资源部天然气水合物重点实验室，中国地质调查局青岛海洋地质研究所，青岛 266071 2.青岛海洋科学与技术国家实验室海洋矿产资源评价与探测技术功能实验室，青岛 2660713.中国地质大学（武汉）工程学院，武汉 430074

天然气水合物是水和甲烷在低温高压环境下形成的类冰状物质，具有分布范围广、资源量大等特点[1-5]。由于水合物主要赋存于沉积物的孔隙空间中，深刻理解其微观孔隙结构特征对于水合物的勘探和开发十分重要[6-7]。由于天然样品获取和保存困难，多数针对含水合物沉积物物性的研究是通过实验室人工合成样品完成的[8]。大量实验结果表明，天然气水合物在沉积物孔隙中呈现非均质分布，且会由于沉积物物性和制样条件的影响呈现不同的赋存特征[9]。天然气水合物生成过程中水合物饱和度的增加，会导致水合物的赋存特征（例如赋存位置、赋存形态和赋存量）发生变化，同时使得样品孔隙尺寸、孔隙分布规律、有效孔隙度、渗透率等宏微观物性特征发生明显改变[10-14]。

核磁共振（nuclear magnetic resonance,NMR）技术是含水合物沉积物微观探测的有效手段，具有响应快、分辨率高、无损伤的特点；其基本测量原理是：不同种类和相态物质中目标质子的含量、晶格状态不同，因而发生核磁共振作用的驰豫响应有所差异[12,15]。Kleinberg等[16]较早地将NMR方法用于气体水合物生成过程的研究中，发现气体水合物更倾向于在尺寸较大的孔隙内生成，气体水合物的生成会导致沉积物渗透率的减小。陈合龙等[17]通过1H谱NMR监测了二氧化碳水合物在石英砂中的生成和分解过程，并依据水合物的驰豫特征推断二氧化碳水合物在孔隙中的赋存形态主要为悬浮型或胶结型[18]。低场核磁共振（low-field nuclear magnetic resonance,LFNMR）技术和NMR技术的测量原理相同，但其灵活度更高、成本更低、测量也更为迅速。Ge等[19]利用LFNMR的横向弛豫时间（T2）特征谱分析了甲烷水合物的生成和分解过程，发现水合物生成过程中T2在大于 10ms的区间内信号下降明显；水合物的饱和度和分解速率与样品孔隙度的相关性较强。Ji等[20]利用相似的方法研究了甲烷水合物分布状态对有效渗透率的影响，发现水相渗透率与水合物饱和度之间的关系受水合物非均质分布情况影响很大。需要指出的是，1H谱NMR和LFNMR在进行甲烷水合物体系的监测时会受到甲烷气中1H信号干扰，尤其是在甲烷压力较高的情况下影响更为严重[21]，这会导致其测量结果不能完全体现水相含量的变化。此外，现有基于NMR或LFNMR的水合物物性研究尚存在两个方面的不足：一是缺少含水合物沉积物微观孔隙结构特征的定量分析和表征，这些特征包括孔隙尺寸、孔隙尺寸分布等；二是诸多渗透率预测模型（包括SDR模型[16,22]、Kozeny-Carman模型[23]、分形模型[24]）均是基于含水合物沉积物的基础物性对储层的渗透率特征进行预测，哪种模型能够更好地体现沉积物孔隙结构与渗透率间的关系尚不明确。

利用氙气进行水合物测试能够有效避免1H信号干扰问题，同时生成速度更快、更易操作。本研究选取氙气水合物为研究对象，分析了3组人造砂样在水合物生成过程中的孔隙结构特征演化规律，具体特征参数包括相饱和度、孔隙尺寸、孔隙尺寸分布等；以实验获得的孔隙结构参数为基础，进一步对比了3种渗透率预测模型（SDR模型、Kozeny-Carman模型、分形模型）在水相渗透率预测方面的表现。所得研究结果能够为含天然气水合物储层评价和开发方案制定提供理论支持。

1 实验材料与方法

本文实验采用自主研发的水合物专用LFNMR测试系统（图1）完成。反应釜腔室由铝合金材料制成的环腔和PEEK材料构成的可拆卸样品塞组成，工作压力上限为15 MPa。反应釜控温仪的控温范围为−10～30 ℃，控温精度为0.1 ℃。实验中使用FC-40型电子氟化液（美国3M公司提供）控制反应釜的温度和环腔压力。实验采用3组由不规则石英砂制成的人造砂样完成（编号为样品1、样品2、样品3）。每种石英砂首先压实填充入热缩管中，两端用PEEK样品塞封堵后，利用热缩管的加热变形性质进行样品密封。3种样品的基础物性参数如表1所示。样品长度的差异对测量结果没有影响。实验采用去离子水和纯度为99.9%的氙气进行水合物样品生成。

在水合物生成实验前，先通过真空饱和法对样品进行水饱和，然后进行初次LFNMR扫描。通过称重法测得的每组样品的初始含水饱和度分别为44.8%、37.5%、42.1%。该状态下的初始饱和度可认为是沉积物的绝对孔隙度。初次扫描后，使用氙气驱替饱和样品，使样品的含水饱和度达到50%～60%。采用过量气法[9]进行氙气水合物生成模拟，模拟的温度和初始压力分别设定为278 K和2.0 MPa。实验过程中实时监测反应釜压力变化，待反应釜压力基本不变时结束实验。LFNMR测量采用CMPG（Car-Purcell-Meiboom-Gill）脉冲序列完成，测量的谱仪频率为21 MHZ，磁场强度为0.5 T。脉冲重复间隔时间和回波间隔时间为5000ms和0.4ms。为保证成像质量，测量的扫描次数和回波数设定为8和10000。单次LFNMR扫描的总耗时约为58 s，实验中每间隔30min进行一次扫描。实验设定下，所得T2驰豫时间对应的孔隙尺寸分辨率约为10 nm。

图1 水合物用LFNMR测试系统设备示意图Fig.1 Schematic diagram of the LFNMR experimental system used for hydrate testing

表1 测量沉积物样品基础物性参数Table 1 Physical properties of sedimentary samples for testing

2 实验结果

图2 饱和水状态下的样品T2特征谱Fig.2 T2 spectrum of the tested samples in the water-saturated state

图3 水合物生成过程中反应釜内的压力变化Fig.3 The changes of pressure in the reactor during the formation of hydrate

图2展示了3组样品在饱和水状态下的LFNMR测量结果，图中纵坐标为水中1H信号振幅，横坐标为T2弛豫时间，该结果反映了样品中所有孔隙空间的弛豫响应。依据LFNMR测量原理，孔隙中流体的驰豫表现受3个方面因素的影响，即自由驰豫、表面驰豫和扩散驰豫。在沉积物孔隙较小、气水共存的情况下，自由驰豫和扩散驰豫的作用可以忽略，介质表面T2驰豫时间可以表征为[15,17]

式中，Vp/Sp为孔隙体积与表面积的比值；在应用中，可以将孔隙视作球形，Vp/Sp=r/3，r为孔隙半径； ρ2为矿物横向表面弛豫率。由公式1可知，在沉积物表面弛豫率不变的条件下，T2驰豫时间与孔隙半径满足正比例关系，即测得的T2驰豫时间越长，对应孔隙的等效半径越大。因而由图2可知，样品1和样品3的最大孔隙半径大于样品2，同时样品2和样品3的最小孔隙半径略小于样品1。图3展示了氙气水合物生成过程中反应釜内压力的变化。由图可知，随着孔隙中水合物的生成，反应釜内压力逐渐下降；样品3的压力下降速度变化不大，样品1、样品2在水合物生成前期下降速率较快，在生成后期压降速率逐渐变慢，表明样品1和样品2在水合物生成前期的生长速率较快；3组样品的压力均最终稳定在0.4 MPa左右，表明该条件下样品内水合物的生成基本达到稳定。

图4展示了氙气水合物生成0～390min时间段内3组样品的T2特征谱连续变化，这些T2特征谱反映了沉积物孔隙水中1H的驰豫信号分布。图中纵坐标为1H的信号强度，下横坐标为T2弛豫时间，上横坐标为根据公式1计算得到的孔隙（被水所占据的孔隙）半径。在孔隙半径的换算中，3组样品的ρ2为30 μm/s[25]。由图4可知，样品1、样品2和样品3中被水占据孔隙的半径分别分布于0.1～96.5，0.06～48.5、0.03～127.6 μm；水合物生成所引起的水相体积变化主要发生在孔隙大于1 μm的孔隙空间中，而小于1 μm孔隙中的水很难生成水合物。造成该现象的原因可能是小孔隙中的水多为束缚水状态，水的活度更低，不利于液态水向固态水合物的转变；或是小孔隙中的水与氙气之间的接触关系较差，缺少生成水合物的条件[26]。需要说明的是，T2特征谱仅能反映样品中被水占据孔隙体积的整体分布情况；在水合物生成过程中，样品孔隙的体积会由于水合物的生成而发生改变，因而通过图4无法求取某特定孔径范围内孔隙中水含量的变化。同时，上文提及的孔径值均是在球形孔隙假设条件下得到的，Vp/Sp取值不同时上述估算的孔径值可能有所差别。

图4 水合物生成过程中T2特征谱变化Fig.4 The changes of T2 spectra during the formation of hydrate in the tested samples

3 讨论与分析

3.1 相体积特征变化

在沉积物颗粒空间形态保持不变的条件下，样品中的含水饱和度（Sw）可根据特定时刻的T2与初始水饱和状态下的T2之间的关系[19]得到（公式2），

式中，Ai和Aini分别表示水合物不同生成时间和初始水饱和状态下的T2曲线与X轴包裹的区域面积。在氙气水合物的生成过程中，部分氙气和水会转化为固态的水合物。根据质量守恒原则，参与水合物生成的水在水合物生成前后的质量保持不变，可以得到如下公式：

式中，氙气的分子质量MXe为131.3，水的分子质量

Mw为 18.0，水的密度（ ρw）为1.0 g/m3；氙气水合物在278 K的温度下的密度（ ρh） 和平均水合物指数（Dhi）采用Takeya等[27]的实验测量结果，分别为1.8 g/m3和5.94。在不同的实验条件下得到的氙气水合物密度和水合指数与所用数值可能有所差别。将上述参数代入公式3，计算得到1体积的水生成氙气水合物后的体积Vr约为1.24体积。据此，水合物饱和度（Sh）和含气饱和度（Sg）的计算公式为：

由上述公式计算得到的氙气水合物生成过程中三相饱和度的结果如图5所示。由图可知，在氙气水合物生成的过程中，水合物饱和度快速上升，含水饱和度快速下降，氙气饱和度也有下降但幅度较小。样品1和样品2在水合物生成的前期水合物饱和度上升较快，后期逐渐放缓；而样品3的水合物饱和度整体上升速率较为均一，后期略有放缓；该结果与图3中的压力变化结果具有较好的一致性。样品1、样品2和样品3在390min时达到的最大水合物饱和度分别为65.4%、52.1%和49.1%。理论上，水合物的生成速率受诸多因素的影响[28]。在本研究中供气充足、初始相平衡压力差一致的条件下，水合物的生成还会受到初始含水分布、介质连通状态的影响，如样品3中相饱和度的变化特征可能是由样品孔隙的连通性较差所致。

图5 水合物生成过程中相饱和度的变化Fig.5 The changes of phase saturation during the formation of hydrate in the tested samples

依据压力变化数据和水合物饱和度数据，可计算得到水合物生成速率与相平衡压差之间的关系（图6）。该处水合物生成速率计算为水合物体积增长量与相邻测量时间间隔（30min）的比值；相平衡压差为实际压力与氙气水合物在测量温度下的相平衡压力（约为0.4 MPa[27]）间的差值。由图6可知，随着水合物生成，样品体系的相平衡压差逐渐减小，水合物的生成速率也逐渐降低；水合物生成速率与相平衡压差之间为正相关关系。Kim等[29-30]提出的Kim-Bishnoi分解动力学模型认为，水合物的分解速率与相平衡压差、分解比表面积成正比。从图6来看，相平衡压差确为水合物生成的主要驱动力之一。水合物的生成速率同时可能受到水合物比表面积的影响，但从T2图谱中无法获取水合物与流体相之间的接触关系，故无法进一步判定水合物生成速率与比表面积间的定量响应。

3.2 孔隙结构特征变化

由于水合物在沉积物孔隙中以固相形式存在，水合物生成过程中产生的水合物会占据流体相空间，引起沉积物孔隙结构的变化。1H谱LFNMR测量能够反映水相孔隙空间信息，但无法区分气体和水合物。因而，下文所进行的孔隙特征分析均是针对被水占据的水相孔隙空间，所得结果并不能完全反映所有流体相空间的结构变化。图7展示了水合物生成过程中水相孔隙空间的最大/小孔隙半径、平均孔隙半径的变化。由图可知，水合物生成过程中水相孔隙空间的最大孔隙半径和平均孔隙半径逐渐减小；最小孔隙半径有先减小后增大的趋势，但整体的振幅小于0.1 μm；水合物生成的整个过程中，3组样品最大孔隙半径的降幅均为75%左右，平均孔隙半径的降幅分别为77.8%、74.2%、58.1%；最大孔隙半径与水合物生成时间之间呈近似对数关系，生成前期最大孔隙半径下降速率较快，后期逐渐变慢。这些孔隙尺寸的变化规律，可以为其他物性预测提供重要的基础参数[31-32]。

为定量阐释水合物生成过程中水相孔隙尺寸分布特征的变化，本文引入分形系数对其进行分析。分形系数是分形理论中最重要的参数，可用于描述不规则、高度复杂体系的形态分布特征[33-35]。本研究从分形系数的基本定义出发，推导了基于NMRT2时间确定分形系数的方法，具体展示如下：

根据分形理论，多孔介质的累计孔隙分布服从[34]

N(＞r)为 半径大于r的孔隙的数目，rmax为最大孔隙半径。对公式5取微分，得到孔隙分布函数f(r)[35]。

将沉积物孔隙视作球状，样品孔隙空间的体积可以计算为

图6 水合物生成速率与相平衡压差之间的关系曲线Fig.6 Relationship between the hydrate formation rate and the phase equilibrium pressure difference

由于T2时间反映了水相体积信号的强度，因而半径为r0的 孔隙体积（Vr0）可以表示为

式中，Ar0和Asum分 别代表半径为r0和所有孔隙T2的信号面积。由此，半径为r0的孔隙的数目为

将公式5两端取对数，得到

将公式1、9代入公式10，得到

式中，C为常数， W′为T2曲线中对应半径＞r0的孔隙所占的体积分数，T2为半径为r0的孔隙所对应的T2振幅。拟合W′与T2在双对数坐标系中曲线的斜率，可确定该状态下孔隙的分形系数。

图7 水合物生成过程中水相孔隙尺寸参数变化Fig.7 The changes of pore-size parameters in the water-phase pores during the formation of hydrate

图8 水合物生长过程中水相孔隙分形系数的变化Fig.8 The changes of fractal dimension in the water-phase pores during the formation of hydrate

图8展示了氙气水合物生成过程中水相孔隙分形系数变化。为了更好地展现不同生成时间孔隙分形系数的差别，图8中删减了测得曲线的数量，相邻曲线的时间间隔为60min。由图可知，在水合物生成的过程中，3组样品中被水所占据孔隙的分形系数呈现逐渐增大的趋势；样品1、样品2和样品3的分形系数变化范围分别为2.197～2.694、2.174～2.808、2.046～2.532。分形系数的变化是由水合物生长的非均质性造成的：水合物主要在沉积物中的较大孔隙中生成，造成这些大孔隙的有效孔隙空间减小；而原有较小孔隙空间体积保持不变，其体积相对比例升高。这种变化反映在图8中表现为：随水合物生成时间的增加，横坐标右侧曲线值减小，而横坐标左侧区域曲线值增大，最终导致整体曲线斜率增加。

3.3 水相渗透率变化

渗透率是水合物开发研究中的关键参数，控制着孔隙空间中的流体流动，会影响水合物开采过程中的压力和产气效率[36]。受限于探测手段和实验方法，现有针对含水合物沉积物渗透率的研究多是针对单相的水相渗透率[35,37]。准确预测含水合物沉积物渗透率的关键是提取能够体现沉积物微观孔隙结构特征的有效参数并建立这些参数与渗透率之间的定量关系。常用的含水合物储层渗透率预测方法包括Schlumberger-Doll Research（SDR）模型[16,22]和Kozeny-Carman（KC）模型[23]。分形方法也可用于含水合物储层的渗透率预测[38-40]。笔者在之前的研究中，推导了含水合物储层的水相渗透率预测分形模型[24]，该模型考虑了水合物饱和度和微观孔隙结构对水相渗透率的影响。本研究将基于前文LFNMR物性分析结果，对比上述3种渗透率模型在水合物生成过程中水相渗透率预测上的表现，所用公式和参数含义见表2。需要说明，由于无法通过LFNMR方法获取气相孔隙空间的分布数据，本部分计算中气体所占据的孔隙体积假设为死体积，在水合物的生成过程中保持不变（图5）且不参与流动；SDR模型中ρ2的取值为30 μm/s；KC模型中的有效孔隙度和分形模型中的分形系数的取值参见章节3.2；为方便对比，对3种模型的水相渗透率结果进行了归一化处理。

表2 不同渗透率计算模型的计算公式和参数含义Table 2 Formulas and parameters employed in different permeability prediction models

图9 基于不同渗透率预测模型得到的水合物生成过程中水相渗透率计算结果Fig.9 Calculation results of water-phase permeability during the formation of hydrate obtained from different permeability prediction models

依据不同预测模型得到的水相渗透率计算结果如图9所示。由图可知，3种方法预测的水相渗透率随水合物饱和度的变化趋势基本一致：随着水合物饱和度的增加，水相渗透率逐渐减小；水合物生成初期，水相渗透率迅速下降，水合物生成后期，水相渗透率减小逐渐变缓。不同之处在于，基于SDR模型和分形模型的水相渗透率预测结果在水合物生成前期的下降速率明显小于KC模型；KC模型在水合物饱和度0～20%范围内下降更为明显，随后基本保持不变。3组样品的计算结果在SDR模型和KC模型中区别不大，而在分形模型中，样品3的水相渗透率在水合物生成前期下降速率更快，样品1的水相渗透率在水合物生成后期下降速率更慢，3组样品的归一化水相渗透率具有明显差别。从模型对比角度来看，分形模型更能体现不同样品中孔隙结构演化的差别。而在SDR模型和KC模型中，水合物生成过程中的孔隙结构变化主要通过T2特征谱的对数平均值和平均有效孔隙度来体现，这种处理不可避免地会引起孔隙结构因素对渗透率的影响被过度简化。因而，通过LFNMR手段获取含水合物沉积物的孔隙结构特征，并利用分形方法进行相渗透率计算，在水合物储层的渗透率预测中可能具有更好的适用性。

3.4 存在的问题及建议

本研究利用LFNMR方法对氙气水合物的生成过程进行了连续监测，并分析了该过程中沉积物的孔隙结构和水相渗透率的演化规律。但由于LFNMR方法无法区分沉积物体系中的气体相和水合物相，研究仅能得到被水占据孔隙空间的具体结构参数，并推测气相和水合物相的整体含量变化。从这个意义上来说，本文的研究尚不完整。基于现有认识，笔者提出以下建议，以期为领域内其他学者的研究提供参考和帮助。

（1）采用LFNMR和计算机断层扫描（CT）联合测量的方法，能够实现含水合物沉积物中各相的区分，并可对水合物的赋存形态特征进行描述，弥补单一测量方法的不足。

（2）在考虑沉积物有效孔隙的结构特征变化的基础上，建立能够体现沉积物中气相、水相的流动表现的两相渗透率模型，对水合物资源的高效开发有更直接的指导意义。

（3）本研究中发现在实验室条件下小孔隙中的水很难参与水合物的生成，并通过公式（1）估算得到该水合物生成孔隙下限约为1 μm；但在天然地质条件下具有更小孔隙半径的泥质粉砂储层也可以作为水合物赋存的重要介质[43-44]。这两者之间的矛盾可能与小孔隙对孔隙水的束缚作用以及水合物形成的时间尺度有关，相关研究还有待进一步开展。

（4）此外，本研究在基于T2谱进行物性特征反演时，沉积物固体界面的表面弛豫率设定为均一值，而在真实情况下水合物-水界面的表面弛豫率与砂-水界面的表面弛豫率可能并不相同[25]。但纯水合物-水界面的表面弛豫率很难直接通过实验方法测得。在未来的工作中，基于更准确的含水合物沉积物界面表面弛豫率表征能够得到更好的物性分析结果。

4 结论

（1）水合物不同生成阶段的水合物生成速率具有明显差别：生成前期水合物生成速率较快，后期逐渐减慢。粒度分布不同的样品中的水合物生成速率变化有所差别。水合物生成过程中，水合物饱和度和含水饱和度变化明显，气相饱和度变化不大。水合物的生成位置主要为沉积物中较大孔隙，半径较小的孔隙中很难生成水合物。

（2）水合物在孔隙中的非均质生长行为，造成沉积物最大孔隙半径和平均孔隙半径随水合物含量增加而减小，同时引起孔隙尺寸分布规律的变化。依据T2谱确定的沉积物孔隙分形系数结果表明，分形系数随水合物含量的增加而增大。

（3）沉积物孔隙尺寸和分布的变化进一步导致沉积物渗透率的变化，水相渗透率随水合物饱和度的增加呈现先迅速下降后缓慢下降。相较于SDR模型和KC模型，分形模型在含水合物沉积物的渗透率预测方面更能体现孔隙结构变化对水相渗透率的影响。