天然气水合物储量大、能量密度高、绿色无污染，是“双碳”时代下极具开发潜力的新能源^[1]。目前，绿色、安全、高效开发天然气水合物已成为新能源研究的热点^[2–5]。然而，在开采过程中由于热量供应不足和焦耳−汤姆逊效应的影响，会二次生成水合物，这不仅会降低储层渗透率进而影响产气效率，还会降低储层的有效应力，影响水合物储层的稳定性，对开采安全构成威胁^[6–8]。因此，研究开采过程中水合物二次生成和控制水合物二次生成的降压方式对保障天然气水合物安全、高效开采具有重要意义^[9]。

大量数值模拟和试验研究表明，在降压开采天然气水合物的过程中，水合物二次生成区域主要集中在水合物分解前缘和井筒附近，并且和压力降低幅度有着极大的关联^{[6, 8, 10]}。降压开采压差越大，二次生成水合物的位置离井筒越近，二次生成水合物的速率越大^[8]。并且，在节流膨胀效应影响下，温度会降至0℃形成冰堵，影响产气效率^[11]。因此，为了在保证产气效率的同时控制水合物二次生成带来的不利影响，降压方式的设计就显得尤为重要^[12–13]。已有试验研究表明，逐步降压能有效避免一步快速降压水合物二次生成和结冰的现象^[14]。Yang Mingjun等人^[7]认为，造成水合物二次生成的主要因素是水合物分解过程中的热量供应不足，而逐步降压能够抑制温度快速降低，提高储层的最低温度，防止二次水合物和冰的生成，并且，压力降低幅度越小，水合物分解速率越高，产气过程越平稳^[15]。Lv Tao等人^[16]的试验结果也表明，逐步降压不但会延长开采进程，而且在开采过程中能提升整体温度，有助于防止二次生成水合物。然而，目前相关数值模拟研究大多关注不同降压方式和储层性质对产气产水的影响，对降压方式如何影响天然气水合物二次生成的研究仍较少^[17–18]。

鉴于此，笔者基于实地数据建立了储层尺度下直井降压开采天然气水合物的数值模拟模型，利用有限差分法求解模型，分析了开采过程中不同降压方式下水合物的二次生成区域、生成时间、时空演化特征以及生成量，探索了控制水合物二次生成的最佳降压方式，旨在为后续试采降压策略的设计提供一定参考。

1.   降压开采天然气水合物数值模拟模型

1.1   几何模型

SHSC-4站位的水合物藏可分成3层：1）水合物+水层；2）水合物+水+气层；3）水+气层^[19]。笔者基于该站位水合物藏的特征和直井的开采方式，建立了直井降压开采水合物圆柱体几何模型，如图1所示。参考该站位的测井数据和岩心分析结果设置模型的地质参数，结果见表1。

图  1  直井降压开采水合物几何模型示意

Figure  1.  Schematic diagram of the geological model used in simulation

下载: 全尺寸图片 幻灯片

表  1  储层初始基本地质参数^[19]

Table  1.  Basic reservoir parameters used for the numerical simulations

储层	孔隙度	渗透率/mD	水饱和度	气体饱和度	水合物饱和度
水合物+水层	0.35	2.9	0.660	0	0.34
水合物+水+气层	0.33	1.5	0.526	0.164	0.31
水+气层	0.32	7.4	0.922	0.078	0

下载: 导出CSV 

| 显示表格

由于没有上下覆层相关的公开数据，其孔隙度、渗透率的取值与其相连的水合物储层一致^[20]。此外，地温梯度44.3 ℃/km，海水的密度和盐度分别为1 035 kg/m³和0.0305，沉积物的密度2 650 kg/m³、本征渗透率10 mD、比热容1 000 J/（kg·℃）^[20–21]。模型共划分6 750个网格。其中，径向r方向上第1列网格的尺寸为0.1 m，用来表示生产井，剩余网格的尺寸随着r增大呈指数增长。

1.2   控制方程

天然气水合物开采过程是一个伴随相变的多相多组分多场耦合的非定常流动，且天然存在的水合物大多为甲烷水合物。因此，作以下假设：1）气相为100%甲烷；2）流动符合达西定律；3）不考虑冰相和水合物相的流动；4）开采过程中沉积物骨架不发生变形。

1.2.1   质量守恒方程

根据质量守恒定律，单位时间、单位体积的质量增量等于单位时间、单位体积的质量流入减去质量流出，由此可得水合物相、气相以及液相的质量守恒方程分别为：

式中：$\phi$为孔隙度；$S$为饱和度；$\rho$为密度，kg/m³；${X}^{\kappa }$为质量分数；$m$为单位时间、单位体积分解或生成的质量kg/（m³·s）；$p$为压力，Pa；${K}_{0}$为大气压下的绝对渗透率，m²；${K}_\mathrm{r}$为相对渗透率，m²；$\mu$为黏度，Pa·s；$q$为注入或生产速率，kg/（m³·s）；$\mathit{g}$为重力加速度，m/s²；下标${\mathrm{H，G，A和R}}$分别表示水合物相、气相、液相和岩石相。

水合物相、气相和液相的饱和度满足，

液相和气相的压力满足，

式中，${p}_{{\mathrm{cap}}}$为毛细管压力，Pa。

毛细管压力由van Genuchten模型^[22]来描述：

式中：${p}_{0}$为大气压力，Pa；${S}_{{\mathrm{mxA}}}$为液体最大饱和度；${S}_{{\mathrm{irA}}}$为束缚水饱和度；λ为van Genuchten指数。

采用修正后的Stone模型描述相对渗透率^[23]：

式中，${S}_{{\mathrm{irG}}}$为残余气体饱和度。

假设多孔介质骨架的绝对渗透率不变，渗透率与水合物饱和度的关系为：

式中：${\phi }_{0}$为无水合物相时的孔隙率；$n$为渗透率衰减指数，通常由孔隙结构决定。

1.2.2   能量守恒方程

热传导方式考虑导热、热辐射和对流，能量守恒方程为：

式中：$U$为比内能，J/kg；${c}_{{\mathrm{R}}}$为干岩的比热容，J/（kg·K）；$T$为温度，K；$\kappa$为热导率，W/（m·K）；${f}_{\sigma }$为表面辐射系数；${\sigma }_{0}$为玻尔兹曼常数；$H$为焓，J/kg；$v$为渗流速度，m/s；${Q}_\text{diss}$为水合物分解吸收的能量，J。

根据Selim的研究，水合物分解的反应热可表示为^[24]：

不考虑压力对水和水合物焓的影响，气体、水、水合物和多孔介质的焓可表示为：

式中，${\sigma }_{{\mathrm{G}}}$为甲烷节流系数。

1.2.3   相平衡方程以及动力学方程

在降压或升温条件下，天然气水合物分解会生成气体和液态水。相应的，在适当条件下气体和液态水也会再次生成天然气水合物。描述上述过程的化学方程式为：

式中，${N}_{{\mathrm{H}}}$为水合物的水和系数。

根据式（13）可得：

式中，$M$为摩尔质量，kg/mol。

利用Moridis基于前人研究数据推导出的回归表达式来描述水合物处于相平衡时的温压条件^[25]。

$T \geqslant 273.20 \text{ K}$时，相平衡压力为：

$T < 273.20 \text{ K}$时，相平衡压力为：

式中，${p}_{{\mathrm{pe}}}$为温度$T$下的相平衡压力，MPa。

单位时间、单位体积水合物分解或二次生成水合物的质量^[26–27]为：

式中：${K}_{0}$为水合物反应速率常数，kg /(m²·Pa·s)；$\Delta {E}_\mathrm{a}$为水合物反应活化能，J/mol；$R$为气体常数，8.314 J/(mol·K)；${F}_{{\mathrm{A}}}$为面积调整因子；$A$为参与反应的表面积，m²；${f}_\mathrm{eq}$为反应平衡逸度，Pa，可用反应平衡压力代为计算；${f}_{v}$为气相逸度，Pa，可用气体压力代为计算。

1.3   初始条件和边界条件

由于储层具有渗透性，孔隙水可以与海水交换，因此将孔隙水压力看作静水压力，其计算公式：

式中：${p}_\mathrm{pw}$为孔隙水压力，$\text{MPa}$；${p}_{{\mathrm{atm}}}$为标准大气压，$\text{MPa}$；${\rho }_\mathrm{sw}$为海水密度，kg/m³；$H$为储层至海平面的距离，m。储层的初始温度分布由地温梯度确定。上、下覆层边界满足狄利克雷边界条件，其温度、压力等物性参数保持恒定。

1.4   降压方式设置

一步降压法（one-step depressurization，OD）的生产压力分别为2，3，4和5 MPa。逐步降压法（step-wise depressurization，SD）设置30 d降压一次，降压幅度从0.25 MPa到4.00 MPa，依次增加0.25 MPa，最低降至2.00 MPa，开采1 800 d。

2.   结果与讨论

2.1   模型求解和模型验证

为求解模型，利用有限差分法对模型进行离散，得到一系列非线性方程组后，采用Newton-Raphson迭代法进行求解^[28–29]。通过对比开采60 d时的数值模拟结果与试采结果来验证模型的可靠性^[30]，如图2所示。由图2可以看出，模拟累计产气量与试采累计产气量基本一致，且变化趋势较为吻合，但开采前期试采累计产气量高于模拟累计产气量，这可能是因为在实际情况中，钻井施工或者储层改造导致生产井附近渗透率较高造成的^[30–31]。综上所述，可以认为该模型具有一定的可靠性。

图  2  模拟生产情况与试采情况的对比

Figure  2.  Simulation results of the gas production compared with the field test data

下载: 全尺寸图片 幻灯片

2.2   水合物二次生成区域

基于第1节构建的数值模型，针对1.4节设置的不同降压方式开展数值模拟，获得了不同降压方式下生产不同时间时储层中水合物的饱和度，如图3所示。由图3可以看出：随着开采进行，部分区域颜色逐渐加深，表明此处有水合物二次生成；不同降压方式下二次生成水合物的位置主要集中在水合物层交界处；水合物+水层的水合物分解前缘的二次生成与直井呈30°倾斜分布，具有中间厚，两边窄的特点，其长度和宽度最大可达23和5 m；水合物层交界处二次生成水合物的范围成水平分布，长度和宽度最大可达100和5 m；开采后期，两处二次生成水合物的范围有在r=17 m处逐渐相交的趋势。对比图3（a）和图3（b）可以发现，一步降压法的生产压力越低，水合物+水层水合物分解前缘处的水合物二次生成范围越大，而水合物层交界处的二次生成范围没有明显区别。对比图3（c）和图3（d）可以发现，逐步降压法的降压幅度越低，水合物+水层水合物分解前缘处的二次生成范围越小，而水合物交界处的二次生成范围越大。对比图3（a）和图3（d）可以发现，当逐步降压法的降压幅度达到4 MPa时，其二次生成位置和区域与一步降压法类似。相较于一步降压法，逐步降压法可以缩小水合物分解前缘处的二次生成范围，但降压幅度过小时，反而会扩大水合物层交界处的二次生成范围。

图  3  水合物饱和度分布图

Figure  3.  Spatial distribution of the hydrate saturation

下载: 全尺寸图片 幻灯片

水合物二次生成是一个动态的过程，图3展示了开采过程中二次水合物不断累积的区域（生成速率＞分解速率），但忽略了水合物分解后又二次生成的情况。因此，基于各网格开采60，180，360和900 d时水合物饱和度差的分布（见图4），进一步探究水合物二次生成区域和位置。由图4可以看出，除了水合物分解前缘和水合物层交界处之外，水合物+水+气层水合物分解前缘也会二次生成水合物，此处二次生成水合物的范围呈竖直分布，长度和宽度最大可达15和5 m。对比图4（a）和图4（b）可以发现，一步降压法的生产压力越低，水合物+水+气层水合物分解前缘处的二次生成范围越大。对比图4（c）和图4（d）可以发现，逐步降压法的降压幅度越小，水合物+水+气层水合物分解前缘处的生成范围越小，且离井筒位置越近。这是因为降压幅度越小，水合物完全分解半径越小，相应二次生成水合物的位置就离井筒越近。对比图4（a）和图4（d）可以发现，尽管其水合物+水+气层水合物分解前缘处二次生成的位置类似，但逐步降压法下二次生成的范围更小，二次生成时间更晚。这是因为逐步降压可以避免采用一步降压时水合物大量分解吸热导致的分解区域温度急剧降低，因而提高了储层平均温度，破坏了二次生成水合物的条件。总的来说，不同降压方式对二次生成水合物的位置没有显著影响，但会影响二次生成范围和二次生成时间。

图  4  不同时间节点水合物饱和度差的分布

Figure  4.  Spatial distribution of the hydrate saturation difference between different time

下载: 全尺寸图片 幻灯片

2.3   水合物二次生成时空演化特征

在图5颜色加深的水合物二次生成区域设置了垂直方向（水合物+水层水合物分解前缘）和水平方向（水合物层交界处）的监测点。模拟监测点在开采前和开采60，180，360，900和1800 d时水合物的饱和度，结果如图6所示。

图  5  监测点设置

Figure  5.  Schematic diagram of monitoring points distribution

下载: 全尺寸图片 幻灯片

由图6（a）可以看出：在时间演化上，水平方向在开采360 d后开始累积二次生成的水合物，主要集中在360 ～900 d累积，开采900 d后累积较少；在空间演化上，离井筒越近，水平方向上二次生成的水合物越少。这是因为水合物+水+气层处水合物分解生成的气体和游离气受浮力影响向上运移，为水合物层交界处提供了充足的气源。由于储层渗透率较低，压降对远离井筒区域的影响较小，其温压条件可以使水合物稳定存在，因而远离井筒处二次生成的水合物更多。然而，当逐步降压法的降压幅度为0.25 MPa时，开采1 200 d后靠近井筒处二次生成的水合物更多。这是因为在压差的影响下，气体会向靠近井筒处运移，较小的水合物分解半径使气体在此处聚集。其次，较低的降压幅度使储层整体压力较高，有利于水合物的二次生成。

如图6（b）所示，在时间演化上，垂直方向上不同降压方式的二次水合物累积时间不同。对一步降压法来说，开采60 d前垂直方向上二次生成的水合物就开始累积，累计时间集中在0～900 d。其生产压力越低，二次生成水合物累积越多，垂直方向上二次生成差异越大。逐步降压法下二次生成水合物的累积时间取决于降压幅度，降压幅度越小，二次生成水合物的累积时间越晚且二次生成的越少。空间演化上，离井筒越近，二次生成的越多。注意到，当一步降压法的生产压力为2 MPa或逐步降压的降压幅度为4 MPa时，开采360和900 d时在−55～−60 m有少量水合物分解，但热量供应不足导致二次生成水合物再次累积。这说明较大的压差虽然可以促进水合物分解，但同时也会在开采后期促进水合物二次生成。总的来说，不同降压方式对垂直方向上水合物二次生成的时空演化特征有明显影响，对水平方向影响较小。逐步降压相较一步降压可以推迟并减少垂直方向上水合物的二次生成。

图  6  监测点水合物饱和度分布时空演化特征分布

Figure  6.  Curves of hydrate saturation variation at monitoring points over time

下载: 全尺寸图片 幻灯片

2.4   水合物二次生成量

为探究不同降压方式下储层中二次生成水合物量随开采进行的变化趋势，采用下式计算二次生成水合物的量。

式中：${m}_\text{sec}$为二次水合物生成量，kg；${v}_{i}$为某一网格的体积，m³；${\phi }_{i}$为某一网格孔隙度；${\Delta \rho }_{i}$为某一网格的水合物密度变化量，kg/m³；${\Delta S}_{i}$为某一网格的水合物饱和度变化量；$q$为参与计算天然气水合物二次生成的网格数量。

基于数值模拟结果，结合式（19）可以计算开采过程中不同降压方式下二次生成水合物的量，结果如图7所示。由图7可以看出：不同降压方式下二次生成水合物量的变化趋势基本一致，但在峰值和峰值到达时间上存在差异；开采初期，储层的温度和压力未受到较大扰动，有利于水合物二次生成，使初期二次生成水合物的量较大；压降进一步传播后，二次生成水合物的量开始减小，但随着水合物不断分解吸热，二次生成的水合物开始逐步累积，这与图6（a）中的二次生成水合物积累的时间一致；开采600 d后，二次生成水合物的量趋于平稳，说明此时储层内部水合物的分解和生成达到了稳定状态；开采1200 d后，二次生成水合物的量缓慢增大，这可能是由于储层整体渗透率较低，到了开采后期热量供应不足导致的^[20]。

图  7  不同降压方式下二次水合物生成量随时间的变化

Figure  7.  Curves of secondary hydrate production over time under different depressurization modes

下载: 全尺寸图片 幻灯片

采用一步降压法时，生产压力越高，二次生成水合物的量就越大，且更先达到峰值。采用逐步降压法时，降压幅度越小，开采前期二次生成水合物的量越大，后期二次生成水合物的量越小。然而，当降压幅度达到2 MPa时，进一步增大降压幅度对二次生成水合物量的影响不大。这可能是因为降压幅度过大时，降压时长基本相同造成的。逐步降压法的降压幅度大于1 MPa时，在开采过程中二次生成水合物的量基本都比一步降压法少。总的来说，不同降压方式对开采过程中二次生成水合物量的变化趋势没有太大影响，但会影响二次生成水合物量的峰值和峰值到达时间。

3.   结　论

1）二次生成水合物的位置主要在水合物分解前缘和水合物层交界处，分别呈倾斜分布和水平分布。一步降压法的生产压力越低，分解前缘的二次生成水合物的范围越大。逐步降压法能够推迟和减小水合物分解前缘二次生成水合物的范围，降压幅度越小效果越好，但最终二次生成水合物的范围受当前生产压力的影响。

2）时间演化上，水合物层交界处二次生成的水合物在开采360 d后开始积累，而分解前缘处二次生成水合物的累积时间取决于降压方式：一步降压法从开采60 d前就开始积累，而逐步降压法能够推迟累积时间，降压幅度越小，推迟越久。空间演化上，水合物层交界处二次生成水合物的位置离井筒越远，压降传播影响越小，二次生成的水合物越多；水合物分解前缘处离井筒越近，二次生成的水合物越多。

3）二次生成水合物的量在开采初期先减小后增大，中期略有下降后平稳，后期由于热量供应不足导致增大。一步降压法下，二次生成水合物的量与降压幅度呈反比。逐步降压法下，开采前期二次生成水合物的量与降压幅度成反比，后期二次生成水合物的量与降压幅度成正比。

4）综合考虑水合物二次生成的区域和生成量，建议设计降压方式时优先选择逐步降压，降压幅度在1～2 MPa。此外，降压时长也会影响水合物的二次生成，后续可考虑研究降压时长的设计。

5）模拟时假设开采过程中沉积物骨架不发生变形，未考虑降压开采过程中由于孔隙变化引起的力学响应，后续应建立耦合地应力场的模型，以研究储层力学响应和水合物二次生成之间的相互影响。