0.   引言

天然气水合物是一种笼形晶状化合物（图 1），具有能量密度高、储量巨大等特点，其中97%以上分布在海底沉积物中^[1-2]。据估算，全球天然气水合物中含碳量相当于全球已探明传统化石燃料总碳量的两倍，被普遍认为是21世纪最具潜力的新型清洁能源之一^[2]。目前世界很多国家都开展了水合物勘探和试采工作，例如美国JIP计划^[3]、日本MHR & DP项目^[4]，挪威GANS项目^[5]。另外，俄罗斯、德国、韩国、印度等国家都相继开展了水合物相关的研究项目和开发计划^[6]。中国在水合物领域的研究起步较晚，但取得了重要成果。在南海北部陆坡圈定了6个水合物成矿远景区，19个成矿区带，25个有利区块和24个钻探目标区^[7-8]。特别是在神狐、东沙等海域钻探获取了水合物实物样品，发现了多种类型的高品位水合物储层，锁定了试采目标^[9]。并于2017年和2020年成功完成了两次水合物试采工作，取得了突破性进展^[10-11]。

图  1  水合物实物与分子结构图^[2]

Figure  1.  Physical and molecular structural diagrams of hydrate

下载: 全尺寸图片 幻灯片

目前水合物领域的研究主要集中在资源储量勘探评估、水合物开发技术和水合物地层稳定性评价3个方面。其中水合物地层稳定性评价是水合物安全可持续生产的前提条件，特别是考虑开采影响下水合物储层的动态响应及海底斜坡稳定性问题。通过海洋工程地质调查、室内模型试验和数值模拟研究发现，含水合物储层的海底斜坡失稳破坏，往往与水合物分解后产生的超孔隙压力有关^[12-14]。因此，近年来专家学者对与水合物分解相关的超孔压计算模型进行了深入研究。例如，Nixon等^[15]提出一个数学模型，通过假定水合物的分解比例，量化了不排水条件下水合物分解产生的超孔压。在此基础上，Xiong等^[16]考虑甲烷的溶解度对孔隙压力的影响，对模型进行了改进。Zhang等^[17]通过RK方程计算了真实甲烷气体的压缩系数，并将此系数引入到上述模型。Chen等^[18]利用此模型，通过考虑水合物分解过程中压力和温度的平衡，计算了水合物分解的最大超孔隙压力。此外，刘锋等^[19]、Grozic^[20]、Zhu等^[21]进行了相关研究。然而，这些研究中所采用的孔隙压力模型很少考虑水合物储层的热-水-化(THC)耦合作用，导致无法准确获取瞬态超孔隙压力及其演化特征，因此很难准确评价海底斜坡的实时动态稳定性。

本文考虑水合物注热开采的影响，发展了一种热-水-化（THC）耦合作用下海底斜坡动态稳定性简化评价方法，并自主开发了相应的数值计算程序，建立了水合物饱和度与超孔隙压力间的内在关系；进而获取了多物理场耦合作用下水合物分解引发的瞬态超孔压演化特征及储层的动态响应；并利用改进的极限平衡法对海底斜坡稳定性进行动态评价；最后以南海北部陆坡水合物试采区典型斜坡为目标，研究了水合物注热开采影响下水合物储层动态响应及实时稳定性演化。

1.   考虑THC耦合作用的超孔压确定方法

1.1   THC耦合方程

含水合物沉积物可视为由多组分（土骨架、水、水合物）和多相（固、液、气）组成的动态系统（图 2）。根据地质调查^[22]，海底水合物大多储存于相对封闭的沉积地层中，其形成由沉积物沉积、孔隙流体流动及底部甲烷供给等因素决定，且这类储层中一般缺乏明显的横向变化，因此只需考虑竖向的流体运移和能量迁移。基于此，采用Xu等^[23]提出的一维简化THC耦合模型描述海底水合物地层的赋存特性。该模型通过求解质量、能量和甲烷守恒方程，给出了水合物储层边界的计算公式，并推导了初始地层温度–压力分布和水合物聚集的表达式。

图  2  水合物沉积物多组分多相孔隙介质示意图

Figure  2.  Schematic diagram of multi-component multi-phase porous sediments containing gas hydrate

下载: 全尺寸图片 幻灯片

首先，根据质量守恒给出一维简化情况下流体运移的控制方程：

式中，k为渗透率，ρ_l，μ_l分别为液体相的密度和动态黏度，P为地层压力，z为空间坐标，g为重力加速度，q_f为总质量通量。

其次，水合物储存系统的能量迁移可通过能量守恒方程描述，其简化模式为

式中，λ为有效导热系数，T为地层温度，C_pl为水合物的比热容，q_e为总能量通量。

最后，储存系统中的甲烷运移可由甲烷介质的质量守恒方程描述：

式中，ϕ为沉积物孔隙度，K_m为扩散系数，S_l为液体相饱和度，M_l为液体相中甲烷浓度，q_m为与深度相关的甲烷通量。

以上3组基本控制方程间的耦合作用，其具体的建立过程为对式（1），（2）进行积分，可得到关于压力剖面P和地层温度剖面T的函数表达式：

式中，ξ为参数，ξ=q_e/(q_fC_pl)，P₀，T₀为海底面处的压力和温度。联立式（4），（5），可以得到任意深度z处的压力与温度之间的关系式。

水合物的形成需要适当的温压条件，地层温度剖面和水合物相平衡曲线共同确定了水合物的理论稳定区间，如图 3所示。Sloan等^[24]通过试验测定了甲烷水合物的相平衡曲线，表达为

图  3  水合物理论稳定区域与生成区域的空间关系示意图^[26]

Figure  3.  Schematic diagram of spatial relationship between hydrate stability zone and hydrate formation zone

下载: 全尺寸图片 幻灯片

式中，T_diss为水合物分解的相平衡温度。

将式（4）~（6）联立，可得到水合物理论稳定区域的上边界z_st和下边界z_sb：

式中，上边界z_st和下边界z_sb见图 3右侧标注。另外，水合物的生成不仅需要适宜的温度和压力，还取决于孔隙中甲烷的可用量。只有当孔隙中甲烷含量大于或等于甲烷溶解度时，才可能在孔隙中形成水合物，如图 3曲线（c）和曲线（d）所示。在水合物生成区域内部，多余的甲烷会与水结合形成水合物，因此孔隙甲烷含量等于甲烷溶解度，可由下式描述：

式中，M_sl为任意温度、压力及孔隙中盐浓度影响下的甲烷溶解度，M_{sl, diss}(T_diss, P, S)为甲烷溶解度关于相平衡温度T_diss、压力P和盐度S的函数，由试验测定^[25]，β=14.4 ℃为经验常数。

而在水合物生成区域上部，孔隙中甲烷含量小于其溶解度，甲烷介质以溶质赋存于孔隙中，可由式（3）描述。为确定水合物生成区域的上边界z_ot和下边界z_ob(图 3右侧标注)，在0→z_ot区间内对式（3）进行积分，并辅以式(4)，(5)，(9)可得

式中，M₀为海底面处的甲烷含量，q_mt为z_ot处的甲烷通量。

综合考虑式（7），（8），（10），（11），同时满足水合物理论稳定区域条件和水合物生成区域条件即可得出水合物真实赋存区域，其上边界THZ和下边界BHZ如图 3所示。

进一步，考虑水合物形成的时间尺度和空间分布与甲烷供应有关，可给出确定整个水合物储层初始水合物饱和度的表达式^[23]。甲烷水合物填充孔隙度ϕ到任意平均饱和度S_h所需的时间t₀近似为

式中，t为时间，∂S_h/∂t为水合物的积累速率，

其中，ρ_h，M_h为水合物密度和水合物中甲烷的质量分数。综上，由式（1）~（13）构建了一套完整的控制方程，可用于计算地层初始状态下的温度、压力剖面和水合物饱和度剖面。

1.2   超孔压计算

水合物储层注热开采会破坏地层中温度的动态平衡，一旦地层中的温度T超过水合物的相平衡温度T_diss，水合物会产生分解。地层温度剖面T的演化可以用一维热传导方程描述：

式中，κ为地层有效热扩散系数，假设为常量，这是因为水合物的导热系数接近于水（相差 < 10%）^[27]，s(z，t)为源项，用以考虑可能的热量增加和消耗。需要注意的是，在地层升温过程中（见图 4（a），A→B），一旦温度超出水合物的稳定区域，水合物会发生分解。这其中包含水合物的相变过程，因此需要消耗一定的潜热。这意味着当地层温度T_{i, j}（下标i代表时间步，j代表某地层深度）达到T_diss后，随着热量的增加，温度会保持T_diss不变，这会延迟温度在地层中的传播（见图 4（b））。

图  4  注热开采导致的地层温度变化路径示意图

Figure  4.  Schematic diagram of change path of stratum temperature caused by hydrate heat injection exploitation

下载: 全尺寸图片 幻灯片

为了计算潜热的消耗，采用热积分法^[28]进行简化。具体实施过程：首先利用有限差分法求解式(14)，计算在任意时间步i时地层深度j的温度T_{i, j}。在水合物分解区域将此温度修正为T_diss，进而会产生一个温度增量dQ_{i, j}=C_ph(T_ij-T_diss)，此温度增量会根据时间步和深度更新。将每一时间步所产生的热量累计求和，可以计算出地层深度j增加的总增量∑dQ_{i, j}。一旦累积的总热量超过水合物的分解潜热，则意味着相变完成，地层温度重新按照式（14）升高。在此过程中，位于深度j处，时刻i时的水合物饱和度S_{h i, j}为

式中，S_{h, j}为深度j处的初始饱和度，L_h为水合物的分解潜热。

根据式（15）可得出在地层温度变化过程中水合物的分解量，

最后，利用气体非逃逸方法计算温度升高过程中瞬态超孔隙压力及其演化。此方法由Nixon等^[15]首先提出，后续不少学者对该方法进行了完善^[16-21]。这里在前人工作的基础上，充分考虑水合物赋存系统的THC耦合作用，获取水合物实时分解量(式(16))，并引入到该方法中，实现升温后地层中的真实超孔隙压力演化。其中做了如下假设：①依据中国南海水合物试采区的地质调查结果，其水合物储层沉积物以细粒为主，渗透系数较小^[9]，因此忽略水合物分解后气体的渗流；②水合物赋存区的孔隙水在水合物分解前已处于甲烷饱和状态，因此忽略分解后生成甲烷在孔隙水中的溶解；③相对于气体的压缩量其余介质的压缩量很小，忽略不计。

基于以上假设，单位水合物分解后沉积物的体积变化等于生成气体的体积∆V_unit，即

式中，V_g为水合物分解产生的气体体积，V_g=164.6P_atmT_nat/(T_atmP_nat)，P，T分别为地层压力和温度，下标atm和nat分别为标准条件和地层真实条件，V_wh为水合物分解产生水的体积。故沉积物总体积变化为

根据有效应力原理，超孔隙压力与沉积物变形有关，

式中，u_{e i, j}为时刻i时深度j处的超孔隙压力，ε_v为沉积物体应变，E_s为压缩模量，其余参数的物理意义同前。综上，可以得到地层升温后考虑THC耦合的超孔隙压力及其演化特征。

需要指出的是，在超孔压计算中做了一些合理性假设，忽略了一部分气体渗流和其余介质的压缩量，故计算出的超孔压会出现大于实际值的情况。另外，水合物分解产生的孔压会增大地层压力P_nat，进而反向影响水合物的分解，使水合物分解减缓，上述超孔压计算方法为单向耦合过程，即不考虑这种超孔压的反向耦合效应，使得海底斜坡的稳定性评价偏于安全。

1.3   方法验证

为了检验上述超孔压计算方法的正确性，复刻了Reagan等^[29]基于TOUGH+HYDRATE的模拟案例。TOUGH+HYDRATE为目前关于水合物分解模拟较为普遍的商业软件之一，并得到了广泛应用^[30-31]。此案例描述了浅层水合物赋存的典型场景，其中水深为320 m，海底温度为0.4℃，地温梯度为3℃/100m，模拟了海底水合物层100 a内，以每年0.03℃上升工况下水合物饱和度的动态响应。

对比文献[29]中设定升温工况下的计算结果，可见不同时刻（0，10，100 a）水合物饱和度的演化曲线具有较好的一致性，如图 5所示。根据式（19）可知，超孔压计算结果的合理性取决于水合物饱和度变化量∆S_{hi, j}的正确性，上述饱和度演化曲线的一致性很好的检验了本文提出的超孔压计算方法的正确性。

图  5  水合物饱和度剖面演化曲线对比图

Figure  5.  Comparison of evolution curves of hydrate saturation profile

下载: 全尺寸图片 幻灯片

2.   海底斜坡动态稳定性评价方法

极限平衡方法(LEM)被广泛用于海底斜坡的稳定性评价，尤其对于长度远大于沉积层厚度的无限坡。海洋工程地质调查表明，大多数水合物储层与海底面平行分布，贮藏于一定厚度的不透水覆盖层下方^[10]，因此将海底斜坡概化为由上覆层、水合物储层和下卧层组成的无限斜坡（图 6（a））。注热开采条件下，水平井附近的水合物开始分解，将水合物的分解前缘视为潜在滑动面，并以此对海底斜坡进行稳定性评价，如图 6所示。

图  6  水合物开采影响下典型海底斜坡稳定性评价示意图^[32]

Figure  6.  Schematic diagram of stability evaluation of typical submarine slope under influences of hydrate exploitation

下载: 全尺寸图片 幻灯片

根据极限平衡方法，海底斜坡的稳定性安全系数等于滑动面上的抗剪强度和剪应力的比值，

式中，${\vec \tau _{\text{r}}}$，${\vec \tau _{\text{d}}}$分别为滑动面处沉积物的抗剪强度和剪应力，dl为单元滑动面长度，通过建立垂直于l轴的力系平衡，可以导出：

式中，φ′为有效内摩擦角，${\sigma '_0}$为上覆有效应力${\sigma '_0}$= $\sigma - {u_{\text{e}}}$，$\sigma$为总上覆压力$\sigma$=$\gamma '$Hcos²$\alpha$，$\gamma '$为沉积物有效重度，H为滑面深度，c′为沉积物有效黏聚力。水合物分解会引起沉积物胶结强度弱化，进而导致黏聚力降低^{[31, 33]}，这里假设黏聚力随着水合物的分解呈线性折减：

式中，${c'_i}$为i时刻的黏聚力，${c'_1}$，${c'_2}$分别为初始黏聚力和最终黏聚力。

将式（21），（22）代入式（20）中，可得出海底斜坡不同深度、不同时刻的瞬态稳定性安全系数F_{s i, j}为

3.   注热开采下海底斜坡稳定性实例分析

3.1   研究区域

选取南海北部陆坡白云凹陷作为研究区，区域整体位于珠江口盆地珠二坳陷内，走向近东西，面积大于2.0×10⁴ km²，是南海北部最大的一个深水凹陷，其北部与番禺低隆起相接，南端是珠江口盆地的南部隆起带^[34]。区域内发育有大量现代海底峡谷（见图 7），沉积物失稳频发。中国地质调查局于2019年10月—2020年4月在此区域进行了第二次天然气水合物试采，实现了连续产气30 d，总产气量86.14×10⁴ m³，大大提高了日产气量和产气总量^[11]。此次试采首次提出了利用水平井开采水合物的工作思路，可有效克服直井开采大压差导致二次水合物和冰生成的风险。由于复杂地质环境影响，开采引起的水合物分解可能在海底斜坡内部产生潜在滑动面，即使在坡度很小的情况下也会引发斜坡破坏。特别是考虑到海底浅表层沉积物强度低^[35-36]，易在外部作用干扰下发生土体失稳，进而可能导致海底滑坡等地质灾害，威胁海洋开采平台、锚固基础和海底管线等水下设施的安全^[37-39]。这里主要探讨水平井注热开采工况下水合物储层的超孔压演化过程及海底斜坡动态稳定性评价。

图  7  研究区位置及区域海底地形条件^{[11, 40]}

Figure  7.  Location of study area and topography of regional seabed

下载: 全尺寸图片 幻灯片

为服务于第二次水合物试采井的轨迹优选，在试采水平井靶点位置附近进行了先导井测试工作，并利用测井得到的P波速度和S波速度预测了水合物饱和度S_h和游离气饱和度S_g分布（图 8）^[41]。分析可见，该区域水合物储层位于海底面以下深度207.8~253.4 m，厚度为45.6 m，平均孔隙度为37.3%，平均水合物饱和度为31.0%^[41]，全井段岩性主要是泥质粉砂，渗透系数较小。另外，海洋地质调查显示此区域水深1000~1500 m，平均地形坡度为3°，海底面温度为3.3~3.7℃，地热梯度为45~67 ℃/km^[9]。基于开采井原位测试数据，构建二维海底斜坡模型，由上覆层、水合物层和下卧层组成，如图 6（a）所示。为达到最大传热效率，将水平井设置在水合物储层中间位置，即H=230.6 m处。设定由于注热影响，地层温度以0.1 K/d的速度匀速增加，总升高温度为10 K，当温度达到设定值后保持不变。模型所用具体参数见表 1。

图  8  先导井测井响应柱状图^[41]

Figure  8.  Histogram logging response of pilot well

下载: 全尺寸图片 幻灯片

表  1  模型参数^{[9, 11, 26, 41-42]}

Table  1.  Model parameters

参数		取值	单位
坡体参数	z (深度)	—	m
	α (坡度)	3	°
	H(水平井位置)	230.6	m
地层参数	T0(海底温度)	3.5	℃
	G/(地热梯度)	56	K/km
	Sh (初始饱和度)	31	%
	k (渗透率)	2.38	mD
	κ (有效热扩散系数)	4.5×10−7	m2/s
	λ (有效导热系数)	3.66	W/(m·K)
	c1′ (初始黏聚力)	800	kPa
	c2′ (最终黏聚力)	30	kPa
	φ′ (内摩擦角)	8.03	°
	ρ (体密度)	1650	kg/m3
	ϕ (孔隙度)	0.37	m3/m3
通用参数	ρl/(液相密度)	1024	kg/m3
	ρh (水合物密度)	930	kg/m3
	Cpl (液体比热)	4.18×103	J/(kg·K)
	Cph/(水合物比热)	2.16×103	J/(kg·K)
	μl (液体动态黏度)	8.87×10−4	kg/(m·s)
	S (液体盐度)	3.5	%
	Lh/(水合物潜热)	56552-16.8/T	J/kg

下载: 导出CSV 

| 显示表格

3.2   水合物储层动态响应分析

图 9为模拟工况下地层温度剖面的演化曲线，其中THZ和BHZ分别表示水合物储层的上顶面和下底面。初始地温曲线（0 a）和水合物相平衡曲线（T_diss）所围成的阴影部分为水合物的理论稳定区间（下底面为z_sb）。可以看到，由于水平井中热量的注入，水合物储层中部的温度开始升高，并逐渐向上下扩散，一旦温度超过水合物的相平衡温度，水合物开始分解（图 10）。可以观察到，在曲线T_diss附近地层的温度出现局部的降低，这是因为水合物分解是吸热过程，消耗了一定量的潜热，潜热的消耗会造成温度传导的延迟。

图  9  地层温度剖面演化曲线

Figure  9.  Evolution curves of stratum temperature profile

下载: 全尺寸图片 幻灯片

图  10  升温过程中水合物储层动态响应

Figure  10.  Dynamic responses of hydrate reservoir during heating process

下载: 全尺寸图片 幻灯片

图 10（a）为水合物饱和度的演化曲线。

可以看到，随着热量的注入储层中部水合物最先开始分解，饱和度逐渐降为0，在地层中形成一条分解带，且随着热量的传递，上下两侧的水合物也开始分解，分解带不断延伸。图 10（b）为水合物分解前缘的演化曲线，可以更加直观地看到水合物的分解及传播过程，其中红色曲线代表水合物分解前缘，黑色曲线为水合物储层的上下界面。可以看到，分解前缘最先产生于储层中部，并不断向两侧移动，直到与水合物储层上下界面重合，意味着水合物全部分解；在储层中部，水合物分解前缘的演化非常迅速，而在储层两侧则相对缓慢。这是因为储层中部靠近水平井热源，此处的水合物分解反应更加剧烈，而如上所述，水合物分解会导致地层中热量的延迟传导，故储层两侧的水合物分解需要更长的时间。

3.3   海底斜坡动态稳定性分析

图 11显示海底斜坡安全系数沿深度的变化曲线。可以看到，初始状态下海底斜坡处于稳定状态，且由于水合物的胶结作用，增加了沉积物的强度，使得初始状态下水合物储层处的安全系数明显高于地层其它位置。但随着水合物分解，海底斜坡的安全系数迅速减小，并达到临界状态。一方面是由于胶结作用的丧失，沉积物强度迅速降低，另一方面水合物分解产生的超孔隙压力减小了地层的有效应力，随着水合物分解区域的不断延伸，不稳定区域也不断扩大。

图  11  海底斜坡安全系数沿深度的变化

Figure  11.  Variation of safety factor of submarine slope along depth

下载: 全尺寸图片 幻灯片

图 12为海底斜坡最小安全系数的演化曲线。可以看到，经过96 d的升温后，地层中的温度达到分解条件，水合物开始分解，此时海底斜坡安全系数开始下降，在开采358 d后安全系数降为1.0，意味着海底斜坡达到临界状态。此时为了后期的持续安全开采，需要调整开采方案。例如：调整施工参数，配合降压法，既可以有效降低地层中的孔隙压力，又可以提高开采效率。

图  12  海底斜坡最小安全系数随时间的变化曲线

Figure  12.  Variation of minimum safety factor of submarine slope with time

下载: 全尺寸图片 幻灯片

4.   结论

（1）考虑注热开采影响下海底水合物储层的THC耦合作用，推导了水合物分解引起的瞬态超孔隙压力求解格式，并基于极限平衡分析理论建立了海底斜坡动态稳定性简化评价方法。

（2）水平井注热开采影响下，水合物分解前缘最先形成于储层中部，并不断向上、下两侧延伸；由于水合物潜热消耗的影响，导致热量延迟传导，使储层上下两侧的水合物分解需要更长时间。

（3）水合物试采区典型海底斜坡初步分析表明，由于水合物的胶结作用，初始状态下水合物储层处的安全系数偏高；升温96 d后水合物开始分解，海底斜坡安全系数迅速降低，并于358 d后到达临界状态；为保证水合物的持续安全开采，提出调整开采方案。

（4）本文基于无限坡理论分析了水合物注热开采影响下海底斜坡的动态稳定性问题，而实际含水合物层的海底斜坡失稳破坏或滑坡模式，可能表现为弧形、盘式或拉裂刺穿等破坏形式^[13]，未来应加强不同破坏模式的海底斜坡动态稳定性评价方法的深入研究。