水合物法C
O
2
封存（HBCS）作为一种潜力巨大的新型碳封存方法，旨在将大量二氧化碳以固态水合物的形式封存在海洋沉积物中。因其具有封存储量大、应用范围广和稳定性高的特点而逐渐引起各国研究人员的关注。有效评估水合物法C
O
2
海底封存的可行性需要精准刻画多孔介质中液态C
O
2
-C
O
2
水合物体系的热-流-化（T-H-C）多场耦合过程，其中包括传热、多相多组分渗流以及C
O
2
水合物的形成与分解。以往的研究已经从理论分析与室内试验证实了海底沉积物中封存C
O
2
的有效性，然而由于理论分析没有完全耦合相关物理过程，室内试验受制于尺度效应，评估场地尺度下水合物法C
O
2
海底封存的长期稳定性有赖于相关数值模拟程序的开发。
为此，我们开发了一套热-流-化（T-H-C）多场耦合模拟代码HBCS_SC v1.0，用于模拟水合物法C
O
2
封存过程中涉及的多相多组分流动、传热以及C
O
2
水合物的生成和分解过程。该模型能够精确描述多相流在多孔介质中的强非线性特性，同时考虑了C
O
2
在水中的溶解以及由于C
O
2
水合物生成和分解引起的盐度变化。与此前的其他模拟器相比，HBCS_SC代码具有以下功能和特点：（i）精确表征液态C
O
2
的热物理性质；（ii）考虑C
O
2
溶解；（iii）考虑动态盐度效应；（iv）支持两种类型的水合物反应模型；（v）采用的FVM离散化方法具备天然的守恒性和稳定性。
我们获取了南海W19站位地层的详尽地质数据，并以此为基础建立了注入-封存数值模型，用于模拟南海海底沉积物中短期C
O
2
注入（5年）和长期C
O
2
封存（100年）过程。我们首次对海底沉积物中C
O
2
水合物为期100年的长期稳定性进行了验证，并详细分析了100年内C
O
2
组分在不同相态下的质量演化。通过对压力、温度、水合物饱和度、液态C
O
2
饱和度、溶解态C
O
2
摩尔分数质量分数及盐度的时空演化分析，加深了对水合物法二氧化碳海底封存过程中的多物理场耦合机制的认识。为阐明地质因素对封存表现的控制作用，我们设计了全因子实验（DOE）敏感性分析，以探讨关键地质因素对液态C
O
2
注入与C
O
2
水合物封存效率的影响。
图1 自主研发程序HBCS_SC v1.0中热-流-化多场耦合机制示意图
水合物法C
O
2
封存是复杂的多物理场耦合过程，涉及多相多组分流动、传热以及C
O
2
水合物的生成与分解。为了精确表征这一过程，我们开发了数值模拟代码HBCS_SC v1.0，该代码能够描述上述热-流-化耦合过程（见图2）。HBCS_SC基于TOUGH+HYDRATE v1.5代码框架开发，该代码框架的准确性已在多个尺度下进行过充分验证。
我们修改了TOUGH+HYDRATE v1.5代码框架中的主变量转换法（PVSM），以适应液态C
O
2
注入和水合物法C
O
2
封存过程中可能发生的相变，并将液态C
O
2
和C
O
2
水合物的相关热物理性质引入代码框架。该代码能够模拟从实验室尺度到实际场地尺度的含C
O
2
水合物多孔介质系统行为。
基于有限体积法（FVM），该代码通过设置对称带状雅可比矩阵求解耦合非线性代数方程组，并采用牛顿迭代法计算每一时间步的各主变量。时间步长可根据每次迭代步的收敛情况动态调整，以确保数值计算的精度和稳定性。相关数值计算流程如图2所示，其中主变量转换法中涉及的不同相态组合之间的转换准则如图3所示。
图3 HBCS_SC v1.0的主变量转换法中不同相态组合的转化准则
2.2 南海泥质粉砂型沉积物中C
O
2
注入-封存数值模型构建
图4展示了C
O
2
注入-封存数值模型的示意图。根据HBCS_SC代码的前期试算，C
O
2
水合物稳定区（CHSZ）的相平衡底界位于Z=-126 m。为了达到C
O
2
水合物形成条件，液态C
O
2
注入段设计为Z=-90 m至Z=-120 m，垂直井总长度为30 m，完全位于相平衡稳定区范围内。注入段顶部设定在海底泥线以下90 m，以尽可能避免C
O
2
泄漏到海洋中。作为探索性研究，我们设计了一个相对较低的恒定C
O
2
注入速率（
Q
CO2
=4.5×10
-3
kg/s，即141.9 t/y）和较短的注入时间（
t
I
=5年），以尽可能避免松散沉积物的力学失稳。
图4 注入-封存数值模型示意图（CHSZ：C
O
2
水合物相平衡稳定区；CHUZ：C
O
2
水合物热力学分解区）
图5 注入阶段（0-5年）
P
，
T
，
S
H
，
S
CO2-L
，
X
A
CO2
和
X
S
的时空演化规律
模拟包括两个阶段：
C
O
2
注入阶段（0-5年）和封存阶段（5-100年）。
图5展示了C
O
2
注入阶段不同时间点压力（
P
）、温度（
T
）、水合物饱和度（
S
H
）、液态C
O
2
饱和度（
S
CO2-L
）、溶解态C
O
2
质量分数（
X
A
CO2
）及盐度（
X
S
）的时空分布演化规律。图中白色虚线表示液态C
O
2
注入区间。
由图5（a）可见，注入井周围形成了高压区，随着C
O
2
的持续注入，该高压区逐渐向径向扩展。由于初始压力梯度的存在，高压区下部压力高于上部。同时，注入井周围出现了沿径向传播的显著高温区[图5（b）]，这是由于液态C
O
2
注入的压缩放热引起的。受初始地热梯度影响，高温区同样表现出下部温度高于上部的分布。
从图5（c）可以看出，由于C
O
2
压缩放热，水合物未随着C
O
2
注入立即生成，而是首先在温度较低的注入井上方形成，随后逐渐向径向和下方扩展，形成低饱和度的C
O
2
水合物盖层。形成的C
O
2
水合物会占据孔隙空间，显著降低了地层的有效渗透率，从而限制C
O
2
向上迁移。如图5（d）所示，注入的C
O
2
随时间推移逐渐向径向运移，但未观察到明显的向上运移趋势。这可能由以下原因共同造成：（a）在海底下的 温度 和 压力 条件下，C
O
2
为高密度液态（970 kg/m
3
），其密度接近孔隙水，密度差引起的浮力作用有限；（b）液态 C
O
2
具有较高的粘度，影响其自身流动特性；（c）低渗透率的C
O
2
水合物盖层显著限制了C
O
2
的向上迁移；（d）泥质粉砂型沉积物渗透率较低（1-10 mD），阻碍了孔隙流体流动。在注入结束时，C
O
2
的最大径向迁移距离为10 m，最大向上运移距离为7 m。
图5（e）显示，由于初始C
O
2
浓度差较大，液态C
O
2
注入后大量溶解于孔隙水中。溶解区随C
O
2
注入逐渐扩展，近井区域的溶解度低于远井区域，这是因为近井区域温度较高。5年时，溶解区的最大径向长度与液态C
O
2
的最大径向运移距离一致（10 m），而其向上运移距离（10 m）大于液态 C
O
2
的最大向上运移距离（7 m）。C
O
2
溶解于孔隙水会导致
X
S
（孔隙水中NaCl质量分数）降低，因此盐度降低区域与C
O
2
溶解区范围基本一致[图5（f）]。此外，C
O
2
水合物生成时存在排盐效应，因此水合物生成区可观察到显著盐度增加，而盐度升高会抑制水合物的形成。
图6 封存阶段（5-100年）
P
，
T
，
S
H
，
S
CO2-L
，
X
A
CO2
和
X
S
的时空演化规律
图6展示了封存阶段不同时间点的
P
、
T
、
S
H
、
S
CO2-L
、
X
A
CO2
和
X
S
的空间分布。由于C
O
2
注入停止，注入井周围的压力和温度快速降低[图6（a），6（b）]。如图6（c）所示，温度的降低为封存阶段大量C
O
2
水合物的生成提供了前提条件。在适宜的温压条件下，C
O
2
水合物在高液态C
O
2
饱和度的区域逐渐向下生成，水合物饱和度显著增加，最大值达0.6。100年时，高饱和度的C
O
2
水合物稳定存在，证实了水合物法C
O
2
封存的长期稳定性。同时，由于水合物生成是放热降压的反应，C
O
2
水合物生成区可观察到显著的温度升高[图6（b），t=10年，20年]和压力降低[图 6（a）]。随后，由于与周围环境的热交换，水合物生成放热导致的高温区逐渐消失[图6（b），t=50年]，最终温度与地层初始温度趋同 [图6（b），t=100年]。
如图6（d）所示，上部致密水合物盖层显著限制了液态C
O
2
的向上运移，使其在封存过程中逐渐转化为水合物相C
O
2
和溶解相C
O
2
。溶解C
O
2
逐渐向上运移，100年时可达最大高度Z=-47 m，远高于液态C
O
2
的向上迁移距离[图6（b）]。溶解C
O
2
的迁移同样超出了C
O
2
水合物盖层的分布范围。这表明，在进行水合物法C
O
2
封存时，应充分关注溶解相C
O
2
的运移距离。由于封存期内C
O
2
水合物的大量生成，盐析效应表现得更显著[图6（f）]，高盐度区减缓了C
O
2
水合物的形成速率。100年时，高盐度区的最大上移高度可达Z=-55 m。
如图7所示，在注入阶段（SI）中，液态C
O
2
质量（
M
CO2-L
）和溶解态C
O
2
质量（
M
CO2-D
）随时间线性增加。在注入阶段，C
O
2
主要以液态形式存在，其次是溶解态，C
O
2
水合物的生成量极少。5年时，液态C
O
2
质量占比（
R
CO2-L
）为69.93%，溶解态C
O
2
质量占比（
R
CO2-D
）为29.24%，而水合物相C
O
2
质量占比（
R
CO2-H
）仅为0.83%，这表明C
O
2
水合物在注入阶段主要起到低渗盖层作用，而非作为主要封存介质。
在封存阶段（SS）中，水合物相C
O
2
质量占比从初始的0.83%迅速升高至最终的45.94%。与此同时，溶解态C
O
2
质量占比略微增加，从初始的29.24%升至最终的38.42%。液态C
O
2
质量占比在封存结束时（100年）时仅为15.64%。这表明，在封存阶段，C
O
2
水合物逐渐转变为主要封存介质，并继续作为低渗盖层阻碍其下方的液态C
O
2
的上移。
图8 （a）水合物相C
O
2
质量占比（
R
CO2-H
），（b）C
O
2
水合物分布面积（
A
CO2-H
）和（c）溶解态C
O
2
分布面积（
A
CO2-D
）的主效应图与帕累托图
图8 展示了全因子实验敏感性分析中水合物相C
O
2
质量占比（
R
CO2-H
），C
O
2
水合物分布面积（
A
CO2-H
）和溶解态C
O
2
分布面积（
A
CO2-D
）的主效应图与帕累托图。由图可知，6个因素均对水合物相C
O
2
质量占比有显著影响，且影响顺序是盐度>热导率>渗透率>渗透率+热导率>盐度+热导率>热导率+盐度；对于C
O
2
水合物分布面积，只有3个单因素具有显著影响，其影响顺序为渗透率>热导率>盐度；而对于溶解态C
O
2
分布面积则存在4个显著影响因素，其相对影响大小为渗透率>盐度>渗透率+盐度>热导率。这些发现揭示了不同关键参数组合的相对重要性，弥补了常规单因素敏感性分析的局限性。
（a）开发了热-流-化多场耦合数值模拟代码，该代码所引入的液态C
O
2
与C
O
2
水合物热物理性质均已经过实验验证；
（b）在C
O
2
注入阶段，液态C
O
2
从注入段沿径向向外运移，同时在局部形成高压高温区域。高温区的存在抑制了早期C
O
2
水合物的形成。5年时，仅有质量占比为0.8%的C
O
2
水合物在注入井上部生成；
（c）在C
O
2
封存阶段，C
O
2
水合物盖层持续向下生长，并消耗注入的液态C
O
2
。低渗C
O
2
水合物盖层的形成有效限制了液态C
O
2
的向上运移，而溶解C
O
2
则呈现出更大运移分布范围；
（d）在100年封存期内，水合物相C
O
2
的质量占比持续增加，并最终达到注入C
O
2
总量的45.9%。在C
O
2
注入阶段，溶解C
O
2
的质量占比显著增加，而在C
O
2
封存阶段仅略微增加至最终的38.4%。100年时，液态C
O
2
的总转化率达到84.3%，其中C
O
2
水合物是主要的封存介质；
（e）敏感性分析结果表明，较低的盐度和较高热导率有利于C
O
