神狐海域水合物储层地质沉降数值模拟研究

摘要

我国南海神狐海域赋存大量天然气水合物资源，能够有效解决我国未来的能源紧缺问题。但由于其储层地质特征复杂，储层物性条件差，开采难度大，开发成本高，商业化开采之前仍需要进行大量的基础研究工作。本研究通过数值模拟手段研究南海神狐海域水合物，根据实际地质参数建立数值模拟模型，采用单一垂直井开采手段，研究不同生产参数对开采的影响，并对地质沉降进行评估，结果表明，井底压力越低，储层打开程度越高，产气速度越快，但是地层沉降越明显，为避免地质灾害，应合理控制生产压力。

关键词

天然气水合物; 数值模拟; 地层; 神狐海域

0 引言

天然气水合物是甲烷被水分子包裹的一种结晶物质^［1］。其具有清洁高效无污染的特点，蕴藏量十分丰富，其有机碳含量相当于化石能源的两倍，是未来理想的替代能源^［2］。

目前水合物的主要开发手段有降压^［3］，注热^［4］，注入抑制剂^［5］，CO₂置换^［6］，固体法^［7］以及应用上述多种开采方式的联合法。上述所有方法中，由于降压法具有高能效，工艺简单，操作和维护成本低，是最经济有效的开采手段，因此该开采方法受到了重视，降压法目前广泛应用于世界各国的现场开采试验。

2017年，我国在南海神狐海域进行试采工作^［8］，位于我国南海西沙海槽与东沙群岛之间海域，有良好的天然水合物成藏地质条件，神狐海域天然气水合物分布区域较广，沉积类型为泥质沉积，储层物性条件较差，连续试开采60 d，最高产量达3.5×10⁴ m³/d，累计产气30.8×10⁴ m³，取得重要科研成果。

由于降压法开采水合物受到了世界各国的重视，且研究难度和技术实现难度较低，另外有成熟的商业数值模拟软件模拟降压开采，因此有关降压研究的比较完善。然而，虽然降压法开采水合物研究完善，但研究都没有考虑地质力学的影响，均假设地层稳定，无沉降，无出砂，开采状态十分理想。目前考虑地质力学的水合物开采研究尚处于初级阶段，有关水合物地质力学的研究主要包括出砂和沉降。对于出砂的研究，Yu等^［9］利用CMG-STARS，根据神狐地质参数建立出砂模型，结果表明生产压差越大，地层渗透率越高，虽然产气效果越好，但出砂更为严重，另外提出了防砂手段并对防砂效果进行评价，最终优选合理的生产参数。Zhu等^［10］基于TOUGH+HYDRATE+Biot代码开发出一套模拟水合物出砂的程序，通过该程序研究了神狐水合物储层条件下的出砂状况，分析了出砂量和出砂状况，证实了防砂的必要性。对于沉降的研究Sun等^［11］通过结合TOUGH+HYDRATE和FLAC3D两款软件，研究了神狐储层条件下长期开采过程中的沉降特征，结果表明渗透率越高、井底压力越低，产气量越高，但地层整体沉降越大。Li等^［12］应用CMG-STARS根据神狐地质参数建立THMC四场耦合模型，初步研究了降压开采过程中地层位移规律，结果表明，井筒所在位置的顶部地层沉降最大，随着开采的进行，地层沉降速度逐渐降低。另外本文研究了井底压力，地层渗透率和岩石杨氏模量对地层位移的影响，结果表明渗透率越高，井底压力和杨氏模量越低，地层沉降越大。总而言之，由于能够模拟水合物地质力学的软件较少（传统TOUGH+HYDRTE软件不具备固体力学模块），且考虑地质力学的代码编程难度较高，目前的模拟研究相对较少。

综上所述，本文通过单一垂直井进行长期降压生产。研究不同降压生产措施下的水合物开采动态，并对地层沉降进行预测。

1 数值模拟

本文中的数值模拟研究采用CMG-STARS数值模拟软件。对于水合物的研究，虽然和TOUGH+HYDRATE相比公开发表的成果较少，但CMG-STARS能够有效模拟多组分多相流，且能够添加化学反应，能够有效模拟水合物的开采过程，模拟水合物方面，该软件模拟得到了有效的验证。

1.1 数学模型

（1）质量守恒方程。应用CMG-STARS数值模拟软件，能够实现水合物开采的模拟计算，对于三相三组分的水合物开采模型，各组分质量守恒方程为：

水： $$\frac{\partial }{\partial t}({\rho }_{\mathrm{w}}\phi S_{\mathrm{w}})=-\nabla ({\rho }_{\mathrm{w}}{v}_{\mathrm{w}})+{\dot{m}}_{\mathrm{w}}+q_{\mathrm{w}}$$

（1）

甲烷： $$\frac{\partial }{\partial t}({\rho }_{\mathrm{g}}\phi S_{\mathrm{g}})=-\nabla ({\rho }_{\mathrm{g}}{v}_{\mathrm{g}})+{\dot{m}}_{\mathrm{g}}+q_{\mathrm{g}}$$

（2）

水合物： $$\frac{\partial }{\partial t}({\rho }_{\mathrm{h}}\phi S_{\mathrm{h}})=-{\dot{m}}_{\mathrm{h}}$$

（3）

式中：ρ_w、ρ_g、ρ_h——水、甲烷、水合物的密度，kg/m³；φ——介质的孔隙度，小数；S_w、S_g、S_h——水、甲烷、水合物的饱和度，小数；ʋ_w、ʋ_g——水、甲烷的体积流速，m/s；ρ_w、ρ_g、ρ_h——单位时间单位体积内水合物分解产生水、甲烷及消耗水合物的质量，kg·m^-3·s^-1；q_w、q_g——单位时间单位体积内注入或产出的水、气质量，kg·m^-3·s^-1。

（2）能量守恒方程应用CMG-STARS数值模拟软件能够是一款专业模拟原油热采的软件，能够有效模拟。不同生产工况下储层温度变化，各其能量守恒方程为：

$$\begin{array}{l}\frac{\partial }{\partial t}(C_{\mathrm{e}\mathrm{f}\mathrm{f}}T)=\nabla ({\lambda }_{\mathrm{e}\mathrm{f}\mathrm{f}}\nabla T)+\nabla \cdot \left[\left(v_{\mathrm{w}}{\rho }_{\mathrm{w}}{C}_{\mathrm{w}}+v_{\mathrm{g}}{\rho }_{\mathrm{g}}{C}_{\mathrm{g}}\right)T\right]-n_{\mathrm{h}}\mathrm{\Delta }{H}_{\mathrm{h}}+q_{\mathrm{w}}{C}_{\mathrm{w}}T+\\ q_{\mathrm{g}}{C}_{\mathrm{g}}T\end{array}$$

（4）

其中，方程左边两项可以表示为：

$$C_{\mathrm{e}\mathrm{f}\mathrm{f}}=(\mathrm{1}-\phi ){\rho }_{\mathrm{R}}{C}_{\mathrm{R}}+\phi S_{\mathrm{w}}{\rho }_{\mathrm{w}}{C}_{\mathrm{w}}+\phi S_{\mathrm{g}}{\rho }_{\mathrm{g}}{C}_{\mathrm{g}}+\phi S_{\mathrm{h}}{\rho }_{\mathrm{h}}{C}_{\mathrm{h}}$$

（5）

$${\lambda }_{\mathrm{e}\mathrm{f}\mathrm{f}}=(\mathrm{1}-\phi ){\lambda }_{\mathrm{R}}+\phi S_{\mathrm{w}}{\lambda }_{\mathrm{w}}+\phi S_{\mathrm{g}}{\lambda }_{\mathrm{g}}+\phi S_{\mathrm{h}}{\lambda }_{\mathrm{h}}$$

（6）

式中：C_g、C_h、C_w、C_R——甲烷气、水合物、水、岩石的比热容，J·kg^-1·K^-1；λ_g、λ_h、λ_w、λ_R——甲烷气、水合物、水、岩石的导热系数，W·m^-1·K^-1；ρ_R——岩石的密度，kg/m³；ΔH——每摩尔水合物分解时吸收的热量，J/mol。

（3）动量守恒方程。在CMG-STARS中，可以使用地质力学模块来计算一些参数，例如应力（σ），应变（ε）和位移（u）。当忽略动量变化时，力平衡方程可表示为^［12-13］：

$$\nabla \cdot \left[C:\frac{\mathrm{1}}{\mathrm{2}}\left(\nabla u+{(\nabla u)}^T\right)\right]={\rho }_{\mathrm{r}}B-\nabla \cdot [(\alpha p-\eta \mathrm{\Delta }T)I]$$

（7）

式中：σ——总应力，Pa；B——颗粒体积力，m/s²；ρ_r——沉积颗粒密度，kg/m³；C——切向刚度张量，Pa；α——毕渥数，无因次；p——压力，Pa；I——单位矩阵，无因次；u——位移，m；ΔT——温差，℃；η——热弹性系数，Pa/℃。

（4）相平衡方程。根据前人研究结果^［14-15］，模拟天然气水合物相平衡的取值如下：k₁为1.617408 × 10⁹，k₂和k₃均为0，k₄为-1414.909123，k₅为-105.25，具体见式（8）：

$$P_{\mathrm{e}}=\mathrm{1.617408}\times {\mathrm{10}}^{\mathrm{9}}{e}^{\frac{-\mathrm{1414.909123}}{T_{\mathrm{e}}+\mathrm{105}}}$$

（8）

式中：P_e——平衡压力；T_e——平衡温度。

（5）分解动力学方程。由于本文假设水合物为I型水合物，水合物分解产生气，水，反应见式（9）。

CH4·nhH2O$$\underset{}{\rightleftarrows }$$CH4＋nhH2O

（9）

式中：n_h——水合物中水分子化学计量系数，取5.75。

参考Kim-Bishnoi模型，根据前人结果^［15］，假定多孔介质中水合物赋存形状为球形，球型颗粒的表面积为A_hs，多孔介质孔隙度为φ，水合物饱和度为S_h，则多孔介质中单位体积水合物有效分解区域面积见式（10）：

$$\frac{d{n}_{\mathrm{h}}}{dt}=k_{\mathrm{d}}{A}_{\mathrm{h}\mathrm{s}}{S}_{\mathrm{w}}{S}_{\mathrm{h}}{\phi }^{\mathrm{2}}{P}_{\mathrm{e}\mathrm{q}}{\mathrm{e}}^{-\frac{\mathrm{\Delta }E}{RT}}(\mathrm{1}-k)=\frac{k_{\mathrm{d}}{A}_{\mathrm{h}\mathrm{s}}}{{\rho }_{\mathrm{w}}{\rho }_{\mathrm{h}}}\cdot \phi S_{\mathrm{w}}{\rho }_{\mathrm{w}}\cdot \phi S_{\mathrm{h}}{\rho }_{\mathrm{h}}\cdot P_{\mathrm{e}\mathrm{q}}{\mathrm{e}}^{-\frac{\mathrm{\Delta }E}{RT}}(\mathrm{1}-k)$$

（10）

上式即CMG中天然气水合物分解速率公式，右边第一项为反应频率因子。

（6）绝对渗透率模型。本研究采用CMG-STARS自带的Kozeny-Carman方程描述水合物开采。Kozeny-Carman方程是第一个用于解决储层孔隙度和渗透率之间关系的方程，最早是由Kozeny于1927年和Carman于1937年提出^［16］，该模型的提出主要是为了研究地下水在多孔介质中的运移规律。经过近百年的发展，该模型广泛用于工程领域。在CMG-STARS中，绝对渗透率模型见式（11）：

$$K=K_{\mathrm{0}}$$$${\left(\frac{{\phi }_{\mathrm{f}}}{{\phi }_{\mathrm{v}}}\right)}^{ck}{\left[\frac{\mathrm{1}-{\phi }_{\mathrm{v}}}{\mathrm{1}-{\phi }_{\mathrm{f}}}\right]}^{\mathrm{2}}$$

（11）

其中：

$${\phi }_{\mathrm{f}}=$$$${\phi }_{\mathrm{0}}$$$$\left(\mathrm{1}-S_{\mathrm{h}}\right)$$

（12）

式中：K——当前渗透率，10^-3 μm²；K₀——初始有效渗透率，10^-3 μm²；ck——渗透率经验系数，常数；$\phi$_v——初始有效孔隙度，小数；$\phi$_f——当前有效孔隙度，小数；$\phi$₀——不含水合物孔隙度，小数；S_h——当前水合物饱和度，小数。

（7）相对渗透率和毛细压力模型。本文所用各相相对渗透率由OPM模型进行计算，该模型是TOUGH+HYDRATE中相对渗透率计算的基本模型^［17］。其通过束缚水饱和度，残余气饱和度以及水相和气相的相对渗透率衰减指数估算相对渗透率，见式（13）：

$$\begin{array}{l}{k}_{\mathrm{r}\mathrm{w}}={\left(S_w^{\mathrm{*}}\right)}^{n_w}\\ k_{\mathrm{r}\mathrm{g}}={\left(S_g^{\mathrm{*}}\right)}^{n_{\mathrm{g}}}\\ S_w^{\mathrm{*}}=\left(S_{\mathrm{w}}-S_{\mathrm{i}\mathrm{r}\mathrm{w}}\right)/\left(\mathrm{1}-S_{\mathrm{i}\mathrm{r}\mathrm{w}}\right)\\ S_g^{\mathrm{*}}=\left(S_g-S_{\mathrm{i}\mathrm{r}\mathrm{g}}\right)/\left(\mathrm{1}-S_{\mathrm{i}\mathrm{r}\mathrm{w}}\right)\end{array}$$

（13）

式中：k_rw、k_rg——水相、气相相对渗透率，无因次；S_w、S_g——水相、气相饱和度，无因次；S_irw、S_irg——束缚水、残余气饱和度，无因次；n_w、n_g——水相、气相相对渗透率衰减指数，无因次。

毛细压力同样借鉴TOUGH+HYDRATE，如式（14）所示：

$$\begin{array}{l}{P}_{\mathrm{o}}=-P_{\mathrm{c}\mathrm{o}}{\left[(S^{\mathrm{*}}{)}^{-\frac{\mathrm{1}}{{\lambda }_{\mathrm{c}}}}-\mathrm{1}\right]}^{\mathrm{1}-{\lambda }_{\mathrm{c}}}\\ S^{\mathrm{*}}=(S_{\mathrm{w}}-S_{\mathrm{i}\mathrm{r}\mathrm{w}})/(\mathrm{1}-S_{\mathrm{i}\mathrm{r}\mathrm{w}})\end{array}$$

（14）

式中：P_o——毛细压力，kPa；λ_c——毛细管迂曲度，无因次；P_co——初始压力，kPa。

综上，CMG-STARS能够较为全面地考虑水合物开采过程中的问题，且有效性得到了多次验证^［14-15，18-20］，因此能够有效研究水合物开采。

1.2 模型构造和网格设计

在本研究中，建立了一个3D笛卡尔储层模型，其长度为500 m（x坐标），厚度为327 m（z坐标），用于数值模拟。根据现场地质资料，2017年进行的真实水合物生产试验，假设整个储层模型位于1495~1572 m的深度，水深1266 m。沿着z坐标，网状区域从上到下分为3个区域：上盖层1294~1495 m，水合物层1495~1572 m和下盖层1572~1622 m。目前已经通过测井和岩心数据等方法确定了不同层位的储层物性参数，即A层（水合物+水）：1495~1530 m；B层（气+水+水合物）：1530~1545 m；和C层（气+水）：1545~1572 m。在B层射孔生产，井底压力设定为5.4 MPa。见表1、图1。

表1 主要储层物性参数

参数	A层	B层	C层	上下盖层
层厚/m	35	15	27	200/50
孔隙度	0.35	0.33	0.32	0.3
绝对渗透率/md	2.9	1.5	7.4	0.5
水饱和度	0.66	0.5175	0.922	1
气饱和度	0	0.1725	0.078	0
水合物饱和度	0.34	0.31	0	0

图1 神狐海域水合物储层地质模型示意

2 模拟结果和讨论

2.1 不同降压幅度对开采影响

对于单一垂直井降压开采法，井底压力是水合物开采过程中核心生产操作参数。通过控制不同的井底压力来控制产气。通常而言，井底压力越低，生产压差越大，流体驱动力越高，产气速率和累积产气量均较高。然而，过低的井底压力容易引发诸多的问题，例如由于储层压力过低，相平衡温度低于冰点，水合物大量分解导致冰的生成，另外局部高压低温的环境会二次生成水合物，均堵塞孔隙，影响开采，另一方面，过低的井底压力会大幅度改变力学状态，造成生产井变形、出砂沉降等一系列问题，因此实际试采中的井底压力不会太低。本研究中，暂不考虑出砂和井筒问题，地质力学方面仅考虑沉降问题，研究了1.0 MPa到基础案例的井底压力。本研究假设实际试采过程中井底压力为5.4 MPa，该压力为基础案例。

图2显示了不同井底压力下的产气速度和累积产气量随时间变化曲线。结果表明，井底压力越低，产气速度越快，累积产气量越高。这是由于一方面，根据水合物的相平衡曲线，相同温度下压力越低，越有利于水合物的分解；另一方面，井底压力越低，生产压差越大，流体驱动力越高，促进了气体由地层向井筒的运移。对于1.0 MPa的生产压力，2500 d时累积产气量约为基础案例的2倍，整体产气速度更快，产气峰值更高，达到了1.5×10⁵  m³/d，但和商业化开采仍有较大的差距。

（a） 产气速度

（b） 累积产气量

图2 不同井底压力下的产气速度和累积产气量随时间变化曲线

图3显示了不同井底压力下的气水比随时间变化曲线，结果表明，井底压力越低，气水比越高，这与不含自由气的III类水合物储层的气水比有所差异，不含自由气的III类水合物储层条件下，井底压力越低，产气速度越快，累积产气量越高，但气水比越低。这是由于所有产气均源自水合物分解，随着井底压力的降低，产气增加幅度小于产水增加幅度，因此井底压力越低，气水比越低。而本研究考虑了储层自由气，在开采初期，此时产气主要来源于自由气，井底压力越低，自由气产气速度的增加幅度大于产水的增加幅度，在开采中后期，此时的产气主要来源于分解气，随着井底压力的降低，分解气产气速度的增加幅度大于产水速度，因此整体气水比较高。由此可见，在现有储层物性条件下，井底压力越低，气水比越高，产气效果越好。

图3 不同井底压力下的气水比随时间变化曲线

图4显示了不同井底压力下沉降量分布，可以看出，储层顶部的沉降量最大，这是由于降压生产以及水合物的开采，地层压实，地层岩石位移向下，导致地层沉降，储层顶部的沉降为每一层的沉降量之和，因此储层顶部沉降量最大。井底压力越低，最大沉降量越高，对海底稳定性的威胁越大。另外，通过图5可以看出，地层沉降速度初期较快，到了中后期沉降速度较慢，该趋势和前人研究相同^［11-12］。

图4 不同生产井井底压力下地层沉降随时间变化垂向分布

图5 不同井底压力下最大沉降量随时间变化曲线

针对神狐海域水合物试采区域的储层物性，根据上述研究，结果表明，降低井底压力能够有效促进水合物分解，有效提升产气速度和气水比，但与此同时地层具有更大的沉降量，增加了潜在的地质灾害风险，因此需要慎重选择合理的井底压力。

2.2 不同打开程度研究

由于上述所有研究均仅有B层沟通生产井，而A层未能沟通生产井，打开程度较低，打开程度是指生产井和储层沟通的长度和整个储层厚度的，其中裸眼完井的打开程度为1，射孔完井的打开程度根据射孔层段厚度决定。本模拟方案假定B层全部打开。

B层开采的主要优势在于存在自由气，直接降压开采的初期产气量较高，短期开采成效大，由于位于储层下部，压实性较好，地层稳定性较高，但由于其低渗透率的特性，导致长期开采效果较差，为了获得更高的开采效益，本节拟对高打开程度进行研究。

图6显示了储层不同打开程度下的产气曲线，可以看出，2500 d时AB层完全打开时的累积产气量约为基础案例的4倍，另外可以看出AB层全射孔的产气峰值和平均气量均高于基础案例，但是产气峰值出现时间更晚，这是由于生产井直接沟通了A层，增加了产能，同时由于A层不含自有气，全部气体来源于水合物，因此初期产气速率较低，但由于B层产气的补充，使其初期仍可以保持一个较高产气速度的水平。

图6 储层不同打开程度的产气曲线

图7显示了不同储层打开程度下的气水比，可以看出，AB层射孔虽然产气速度提升较大，但初期气水比（0~1000 d）较低，这是由于储层打开了A层，导致产水增加的缘故，又因为A层不含自由气，产气来源与水合物的分解，初期A层水合物分解产气量远小于A层中水产出量，因此初期气水比较低，而到了开采中后期，A层水合物逐渐分解，流动性增加，含气饱和度上升，气相相对渗透率增加，产水速度的下降幅度大于产气速度，因此气水比逐渐增加。由此可见，储层储层全射孔生产能够有效增加长期开采效果。

图7 不同储层打开程度气水比

图8显示了不同储层打开程度下的海底最大沉降量，可以看出，储层打开程度越大，地层沉降越大，开采初期地质沉降速度较快，开采后期地质沉降速度有所减缓。至开采2500 d时，AB层射孔的最大沉降量增长了16%，整体而言，开采和井底压力的影响趋势相同，即一方面促进了产气，另一方面，增大了潜在的地质灾害风险。

图8 不同储层打开程度下海底最大沉降量

3 结论

（1）研究了现有打开程度下，不同生产井井底压力对开采的影响，结果表明，井底压力越低，生产压差越高，驱动力越大，累积产气量越高，气水比越高，产气效果越好，但地层沉降量越大，地质灾害风险增加，因此需要合理设计井底压力，避免地质灾害。

（2）研究了增加打开程度对开采影响，结果表明，射孔层位越多，产气效果越好，而问题同样是沉降量过大，地质灾害风险增加。

（3）根据现有模拟结果，完全打开储层和采用最低井底压力开采具有最大的产气量，但与此同时地层沉降量最大，容易引发地质灾害，因此不建议采用，本文不再具体研究该方案。