二氧化碳（CO₂）排放是导致全球一系列气候问题的最主要影响因素^[1-4]。目前，为了降低 CO₂排放带来的负面影响，中国政府积极推进“碳达峰-碳中和”政策，制定了一系列“双碳”行动方案^[5-7]。对于传统油气行业，废弃气藏CO₂地质埋存是实现“碳中和”经济化的重要途径^[6]。废弃气藏经过长期开发后存在稳定的地质结构，同时大量的地层水可以为 CO₂提供充足的封存空间。

CO₂埋存是将空气中捕集的 CO₂进行提纯处理后注入深部地层，通过封闭的地质构造体、束缚气、溶解气及地球化学反应，实现 CO₂在地层中的永久埋存，也称为碳捕集与封存技术（CCS，Carbon Cap⁃ ture and Storage）^[8-10]。溶解气封存是 CO₂埋存最重要机理之一^[11-12]。CO₂注入深部地层后，通过扩散作用逐渐溶解到咸水中，引起咸水密度发生变化，从而引发局部动力学不稳定，产生重力驱动的指进现象。与单纯的扩散相比，该对流过程可以大大地增加CO₂在咸水中的溶解度。

废弃气藏中存在一定量的不具备经济效益的残余气，如氮气（N₂）、甲烷（CH₄）及硫化氢（H₂S）等。空气中捕集的CO₂提纯后往往也存在微量的杂质气体。CO₂埋存过程中，往往将这些杂质气体一起注入地层，从而降低提纯的经济成本。相关研究表明，杂质气体的存在会影响 CO₂在咸水层中的溶解度^[13]。目前，鲜有定量评价不同杂质气体对 CO₂溶解度影响的相关研究报道。针对该问题，研究了不同杂质气体对 CO₂扩散能力的影响，定量评价了不同混合体系的溶解速率。笔者首先建立了一种考虑逸度-活度的热力学模型，准确捕捉气水混合体系的相态行为；进而通过基于算子的线性化方法高效准确地求解质量守恒方程，计算不同杂质气体对 CO₂溶解速率的影响；最后基于荷兰鹿特丹西佩尔尼斯气藏，开展废弃气藏 CO₂埋存的可行性分析。该研究成果为废弃气藏及深部咸水层CO₂埋存提供了一定的技术支持。

1 热力学模型

热力学模型用来描述液相与非液相之间的热力学平衡状态。液相可溶解部分气体，气相为多组分混合物，主要以气态、超临界态或者凝析态存在。数值模拟的过程中，通常假设局部瞬时平衡，因此相态的计算独立于流动及传质过程。在多相体系求解过程中，牛顿迭代的过程需要准确求解热力学模型，即：

传统的 Peng-Robinson 状态方程仅适用于油气系统，无法准确捕捉CO₂在水中的溶解。因此，采用逸度-活度模型来求解热力学模型：气组分在气相中采用逸度方程表达，气组分在液相中采用活度表达，即：

在热力学平衡系统中需要满足$f_c^{\mathrm{g}}=f_c^{\mathrm{w}}$，则组分 c的相平衡常数为：

热平衡条件下，水组分的平衡常数通过Spycher 模型^[14] 表示：

采用逐次替换迭代法求解热力学模型，通过不断更新逸度系数获取新的K值^[15]。当热力学模型收敛后，即可求得相应的与质量守恒方程相关的热物性参数（如黏度、密度等）。运行 CO₂埋存模型的过程中，每个网格在牛顿迭代的过程中都需要运行该过程。研究表明，对于多组分复杂系统，求解热力学模型所消耗的时间超过总模拟时间的 50%。因此，对于多组分模拟，需要寻求一种提高效率且能满足精度的高效数值模拟方法。

2 多组分数学模型建立及离散求解过程

2.1 数学模型的建立

对于多组分混合物，其质量守恒方程为：

采用达西定律来描述重力及毛管力存在情况下两相体系的流动：

对于可压缩岩石，孔隙度的变化通过岩石的压缩系数描述：

2.2 数值模型的建立及求解方法

为准确捕捉 CO₂在咸水层中的溶解过程，往往需要建立精细的网格模型。ELENIUS 通过一系列敏感性研究发现^[16]，油藏网格尺度在 1 m 时可以获得数值上收敛的结果。那么对于大尺度模型来讲，高精度会导致模型计算量大大增加，因此需要采用高效的数值模拟方法。采用基于算子的线性化方法（Operator-Based Linearization approach，简称 OBL）^[17]，即将质量守恒方程中仅取决于物理状态的物理量组合为一个整体算子。因此（9）式离散之后的算子形式为：

相应的算子表达式分别为：

传统的数值模拟需要计算每个时间步下每次牛顿迭代的物理属性及相应的导数，该过程往往计算效率较低。基于算子的线性化方法提供了一种新的高效线性化方法，如图1 所示。数值模拟运行之前，将连续的未知变量空间（本研究采用压力p及组分 c、总摩尔分数 z_c）进行离散。由于各个算子仅取决于状态变量，因此可以计算各个算子在未知变量空间的离散分布，即离散空间内未知变量交点的算子值。数值模拟的过程中，根据当前未知状态，在算子空间内线性插值得到算子值及相应的导数。该方法大大地简化了雅克比矩阵的构建过程，同时提高了数值模拟的效率。

2.3 数值模型参数

采用荷兰代尔夫特理工大学开发的基于 OBL 方法的开源软件 DARTS^[17-19]，建立 CO₂在咸水中溶解的模型。二维数值模型网格步长（dx，dz）=（ 0.5 m，0.5 m），网格总数为 200×100=2×10⁴。模型横向长度为 100 m，纵向厚度为 50 m，埋深为 2 475 m。孔隙度为0.25，渗透率为100 mD，温度为75℃，纵向温度梯度忽略不计。二维数值模型如图2 所示，模型初始饱和纯水，底部及侧面为封闭边界，顶部为饱和 CO₂的定浓度边界（饱和 CO₂的单一液相，z_c为常数）。图3为模型所用的相对渗透率曲线。

三维模型参考荷兰鹿特丹西佩尔尼斯气藏（图4），底部存在大量的水体，顶部为一气藏，主要含甲烷气。为简化模型，假设该气藏为纯甲烷气藏。工区长度为 1 395 m，宽度为 4 500 m，厚度为 161.8~221.7 m，埋深为 2 610~3 206 m。渗透率为 0.24~1 600 mD，孔隙度为1×10^-4~0.27，模型中部存在层状泥页岩。自2002年开始生产，气藏原始地质储量为6 226×10⁶ m³。为减少模拟时间，数值模型采用45.0 m×30.0 m×2.0 m 的网格，网格总数为 31×150×87= 404 550。首先采用衰竭式开发将油藏压力降至废弃压力，然后从底部水体注入 CO₂，监测 CO₂运移规律及相应的埋存量。相对渗透率模型使用 Corey模型^[20]，具体参数见表1。

3 结果分析

3.1 不同杂质气体对CO₂在咸水中溶解度的影响

CO₂在咸水中的溶解度受压力、温度及矿化度的影响。由不同温度及矿化度条件下CO₂溶解度随压力的变化（图5）可见，压力越高，CO₂越充分地溶解于咸水，溶解度相应地越高。但是不同温度条件下，CO₂溶解度差异较大。室温（25℃）时，CO₂以气态存在，随着压力的升高，CO₂压缩后与水作用变强，因此其在水中的溶解度迅速增加；一旦压力高于超临界压力，CO₂溶解度出现拐点，溶解度增速变缓。当温度升高到90℃时，CO₂溶解度急剧降低，主要是由于温度升高加速了水中 CO₂分子的运动能力，导致溶解的 CO₂从水中再次逸出。25℃下的溶解度约为90℃下的2倍（图5）。当压力高于超临界压力时，拐点消失，呈现平滑的过渡曲线，说明在高温条件下，压力的变化不会对CO₂-咸水系统产生剧烈的影响。矿化度是影响CO₂溶解度的另一重要因素。相同温度条件下，矿化度为2 mg/L时，CO₂在咸水中的溶解度大大降低，且随着压力的升高，降低幅度越大（图5）。

由不同杂质气体对 CO₂在纯水中（矿化度为 0） 溶解度的影响（图6）可见，低压状态下，杂质气体对 CO₂溶解度影响较小，可忽略不计；随着压力的升高，杂质气体的存在使CO₂的溶解度逐渐降低，尤其是CO₂处于超临界状态之后。在相同浓度下，H₂S的存在对 CO₂的溶解度影响较小，其次为 N₂，CH₄降低 CO₂溶解度幅度最大，降低幅度约为 12%。从结果可以看出，少量的杂质气体对 CO₂溶解度影响幅度较小，因此在实际实施过程中，允许少量的杂质气体，从而可以降低CO₂提纯的经济成本。

图5与图6分别对比了改进的EOS与NIST^[21] 和 Mao/Duan 模型^[22]。新的模型可以相对准确地预测 CO₂溶解度随不同影响因素的变化。图5表明，在高矿化度及高压下，该模型与参考值存在一定的差距。AAVATSMARK 等研究指出，在高压条件下， Mao/Duan 模型低估了超临界 CO₂的热力学参数^[23]。本研究所选模型压力不超过30 MPa，该范围内不同模型差异较小，因此并未对此做详细探讨。

3.2 CO₂在咸水层中的指进现象

CO₂注入地层水后涉及复杂的动力学机理：扩散与对流。由于气水密度差异，CO₂会超覆到顶部，从而形成 CO₂富气区。一部分 CO₂逐渐溶解到下部咸水层，导致密度分层及局部密度差异（扩散过程）。溶解了CO₂的咸水密度增大，会触发重力不稳定效应，从而产生有重力主导的对流现象，即溶解了 CO₂的重咸水向下运移，未溶解 CO₂的轻咸水向上运移。该对流过程可以大大提高 CO₂的溶解速率。

由纯 CO₂及不同混合体系下 CO₂指进发育模式 （图7）可见，在重力主导的对流过程中，CO₂羽流逐渐扩展并向下移动，逐渐与周围的羽流混合形成大的羽流。在 100 a 时，纯 CO₂指进速度最快，羽流已经到达模型底部，说明 CO₂的溶解能力最强。溶解少量 H₂S 时，CO₂羽流也接近模型底部，但羽流相对细小。溶解少量 N₂时，移动最快的 CO₂羽流并未到达模型底部，上部羽流相对较粗，但彼此间混合速度慢。当CO₂羽流中含有少量CH₄时，其降低了CO₂ 在水中的溶解度，减缓了CO₂的指进过程，因此羽流移动速度最慢。

将CO₂的溶解速率定义为单位时间内单位面积 CO₂溶解到单一液相中的质量，其表达式为:

由纯 CO₂及不同混合体系下 CO₂溶解速率随时间的变化（图8）可以看出，不同混合体系的成分影响对流触发时间（即溶解速率开始增加时的拐点）。初始时刻，传质过程主要由浓度差引起的扩散主导，CO₂摩尔浓度呈现层状分布。随着扩散的进行，层与层之间的浓度差逐渐变小，扩散作用变弱，表现为溶解速率逐渐降低（黑色箭头）。由于扩散作用，导致局部咸水密度差异，从而引发由重力主导的对流，CO₂溶解速率开始逐渐增加（绿色箭头）。由于 CH₄的存在，降低了 CO₂在水中的溶解度，因此少量CH₄气体的存在会延长对流触发时间；少量H₂S 的存在几乎不影响对流触发时间，N₂次之。由图8 可以看出，CO₂+5%CH₄混合体系对流触发时间是纯 CO₂的 2.5 倍。当对流作用主导传质过程时，CO₂羽流逐渐发育并与相邻羽流相互融合。该过程中， CO₂溶解速率相对较为稳定（灰色箭头）；同时可以看出，不同杂质气体对CO₂溶解速度影响相对较小。由于模型底部为封闭边界，当 CO₂羽流到达底部边界时，开始横向扩展，模型内 CO₂浓度逐渐升高，咸水密度差异降低，对流作用变弱，呈现出CO₂溶解速率降低的规律（深黄色箭头）。总体来说，少量杂质气体的存在影响 CO₂指进过程，延长了对流触发时间，但是对最终溶解量影响不大。

4 废弃气藏CO₂埋存可行性分析

西佩尔尼斯气藏部署 PW W2，PW W4 及 PW W5三口生产井，采取降压开发，平均日产量分别为 10.2×10⁴，15.5×10⁴，54.3×10⁴ m³ /d。PW W1 为注水井，日注水量为111 m³ /d。经现场评价，当地层压力降至 3.5 MPa 时，该气藏将失去经济效益。图9a 和 9c 为原始状态下含气饱和度分布及压力分布。模拟至 2024 年，全区压力降至废弃压力，此后停止生产，含气饱和度及压力分布见图9b和9d。由于气体的膨胀，全区含气饱和度相对较高，但油藏压力较低，此时总体残余含气已无经济价值。废弃时，整体累积采收率达到95%。

自 2025年开始，在底部水层设置一口 CO₂注入井。经现场报告分析，常温下当注入压力达到 29.82 MPa 时可能诱发岩石破裂。现场 CO₂注入温度为40℃，同时考虑到焦耳-汤普森冷却效应，近井地带所能承受的注入压力低于常温下的破裂压力。研究表明，当CO₂冷却参数为6℃/MPa时，该模型运行 10 a 后，近井端压力降低 0.9 MPa，温度降低5.4℃，因此，CO₂实际的注入温度仅为 34.6℃，此时实际压力降为 5.3 MPa，实际的岩石破裂压力仅为 24.52 MPa。将此破裂压力作为模型的约束压力，第 10 a 时，CO₂注气量接近 0，总的注气量为 14.69×10⁶ t，此时关闭注入井，监测CO₂的运移规律。

由 CO₂羽流随时间的运移规律（图10）可见，第 2 a 时（图10b），CO₂未到达气区，因此呈现 CO₂驱替地层水，地层水驱替甲烷气的现象。随着 CO₂的不断注入，第 10 a 时，CO₂与气驱剩余气混合，含气面积大大增加（图10c）。关井后 30 a 时（图10d），CO₂ 主要以游离气的形式存在，并出现一部分残余气与少量溶解气。经计算，关井 30 a 后，游离气占 74.2%，残余气占19.1%，溶解气仅占6.7%。

5 结论

相同温度及矿化度条件下，压力越高，CO₂溶解度越高；但溶解度随温度及矿化度的升高而降低； 少量 H₂S 对 CO₂溶解度影响较小，N₂次之，CH₄降低 CO₂溶解度的幅度最大。杂质气体的存在会影响 CO₂羽流扩散及对流作用。少量的杂质气体会减缓 CO₂羽流在咸水中的移动速度，延长对流触发时间。主要是由于杂质气体的存在一定程度上降低了CO₂ 的溶解度。从长时间尺度来看，杂质气体对总的 CO₂溶解量影响不大。废弃气藏开展 CO₂埋存潜力巨大，是CO₂地质封存的理想场所之一。短时间内， CO₂在水体中主要以游离气及束缚气的形式存在，但随着时间的推移，溶解量会逐渐增多。

符号解释

c——组分，c=1，2，···，n_c；

$\stackrel{-}{c}$——CO₂在水相中的平均物质的量浓度，mol/mol；

c_r ——岩石的压缩系数，1/Pa；

dx，dz——x和z方向的网格步长，m；

D——埋藏深度，m；

D_cj ——组分c在j相的扩散系数，m² /d；

f_c——组分c在气液两相中的逸度；

$f_c^{\mathrm{g}}，f_c^{\mathrm{w}}$——组分c在气相和液相中的逸度和活度；

F——CO₂溶解速率，kg/（m² ·a）；

g——加速度常数，m/s²；

h——亨利溶解常数；

h_c——组分c的亨利溶解常数；

j——体系内的相数，j=1，2，···，n_p，可表示非润湿相（超临界CO₂气相）和润湿相（水相）；

J_cj——组分c在j相内菲克扩散流量，m³ /d；

K——渗透率张量，mD；

K_c——组分c的相平衡常数；

$K_{{\mathrm{H}}_2\mathrm{O}}^{}$——水组分的平衡常数；

$K_{{\mathrm{H}}_2\mathrm{O}}^0$——水组分在参考压力为1×10⁵ Pa下的平衡常数；

K_rge，K_rwe——气相和水相相对渗透率端点值；

K_rj ——j相相对渗透率；

l——离散单元界面；

L（i）——离散单元界面集合；

m——咸水质量浓度，mg/L；

n_c——体系内组分数；

n_g，n_w——气相和水相渗透率模型指数；

n_p——体系内相数；

p——系统压力，MPa；

p_c——毛管力，Pa；

p_j ——j相压力，Pa；

p_n——非润湿相压力，Pa；

p_ref ——参考压力，Pa；

p_w——润湿相压力，Pa；

q_j ——j相体积流量，m³ /d；

R——气体摩尔常数，J/（kg·K）；

S_gr——残余气饱和度；

S_j——j相饱和度；

S_wc——束缚水饱和度；

t——时间，d；

T——温度，K；

u_j ——达西速度，m/d；

v_j——j相在多相体系中的摩尔分数；

V——单元体体积，m³；

$V_{{\mathrm{H}}_2\mathrm{O}}^{}$——水的摩尔体积，L/mol；

x_c，y_c——组分c在液相和气相中的摩尔分数；

x_cg——组分c在气相中的摩尔分数；

x_cj ——组分c在j相中的摩尔分数；

x_cw——组分c在水相中的摩尔分数；

x_g——气相总摩尔分数；

$x_{{\mathrm{H}}_2\mathrm{O}}^{}$——水组分在水相中的摩尔分数；

x_w——液相总摩尔分数；

$y_{{\mathrm{H}}_2\mathrm{O}}^{}$——水组分在气相中的摩尔分数；

z——某一组分总摩尔分数；

z_c——多相体系中组分c总摩尔分数；

α_c（ω）——当前时间步的累积项算子，mol/m³；

α_c（ω_n）——前一时间步的累积项算子，mol/m³；

β_c，β_cj（ω）——流动项算子，mol/（m³ ·mPa·s）；

γ_j（ ω）——饱和度状态算子，1/MPa；

κ_c——活度系数；

μ_j——j相黏度，mPa•s；

ρ_j ——j相摩尔密度，mol/m³；

ϕ——孔隙度；

ϕ₀——油藏的原始孔隙度；

${\Phi }_{{\mathrm{H}}_2\mathrm{O}}$——水组分在气相中的逸度系数；

$\mathrm{\Delta }{\mathrm{\Phi }}_j^l$——离散单元界面l上的压力差，Pa；

Γ^l ——离散单元界面l传导系数，mD/m³；

Γ^l_d——离散单元界面l的扩散传质系数，1/m³；

χ_cj（ω）——扩散状态算子，mol/（m·d）；

Δχ_cj ——离散单元界面上组分 c 在 j 相中扩散状态算子差，mol/（m·d）；

ψ_c——组分c在气相中的逸度系数；

ω——表示物理状态的量，表示压力、饱和度及温度；

ω_n——前一时间步的物理状态。

参考文献