滞后现象对CO2深部咸水层地质封存的影响

郭朝斌，张可霓，鲁维丰，凌璐璐，郭振东

（1.北京师范大学水科学研究院，北京 100875；2.神华鄂尔多斯煤制油分公司，内蒙古 鄂尔多斯伊旗 017209）

温室气体排放增加导致温室效应加剧，目前全球变暖已经成为国际社会关注的热点问题。CO2的封存问题就是在这样的背景下提出来的，CO2的捕获和封存技术目前已经成为一种国际公认的减排措施，被认为有巨大的减排潜力［1］。地下咸水层在世界各地广泛分布，而且规模大，是切实可行和最具发展前景的减排途径，被认为是长期地质封存的最有效方法之一［2］。

近年来国内外很多学者对CO2地质封存的数值模拟陆续展开了研究。R.Juanes等［3］通过对CO2捕获机制和相对渗透率的研究及数值模拟，表明交替注入水和CO2可以提高注入效率。注入的水使大面积CO2晕分裂，增大吸湿过程的时间和面积从而增强CO2的捕获和固定，但这种方法会导致井孔底部压力增大，受到盖层的封闭性能、管理和经济方面等的限制。理解大体积CO2晕的流体力学行为对提高咸水层中CO2地质封存技术非常重要，Christine Doughty［4］对美国加州圣华金河峡谷南部Kimberlina发电厂进行数学建模研究，利用TOUGH2对倾斜地层中CO2地质封存大规模初步试验中CO2晕的演变进行模拟，结果表明在4年注入周期后，注入的100×104t CO2中约20%为溶液相，25年后，CO2晕达到稳定，溶液相态比例上升到38%，62%为超临界状态，其中只有3%为自由态。随着时间继续，晕的演变仅表现为缓慢的溶解作用，并伴有饱和溶解CO2咸水的对流。国内对提高CO2封存能力、安全评估等方面研究中，许雅琴等［5］通过模拟提高低渗地层中渗透率等探讨如何提高咸水层CO2封存注入率，但对滞后效应在CO2封存方面的研究相对较少。

1 方法和理论

1.1 超临界状态CO2的吸湿和疏干过程

将超临界状态CO2封存到深部咸水层会涉及到诸多多相流问题，注入的CO2为非润湿相，水及少量溶解到水中的CO2作为润湿相。由于润湿相和非润湿相流体性质不同，运移过程中因接触角不同和残余气相的产生会对流体运移产生影响。注入时CO2晕连续增长，超临界状态CO2(表征为气相)作为非润湿相驱替含水层中润湿相(咸水)，即疏干过程。对于随后停止注入阶段，当CO2晕在浮力作用下向上运移(图1)，不同位置经历过程不同，在运移前端为疏干过程，在晕尾端，润湿相驱替非润湿相，为吸湿过程，有部分残余CO2被捕获，变为不可运动的残余气相，进而影响相对渗透率和毛细压力的计算。不同位置和过程中残余气相不同，使得相同饱和度下相对渗透率和毛细压力不相同，在特征曲线上表现出滞后现象。

图1 CO2向上运移过程晕后边缘残余示意图［3］Fig.1 Schematic diagram showing the trail of residual CO2as the plume migrates upward

1.2 数值模拟方法

数学模型和数值模拟在评估地质封存CO2可行性上具有重要的作用，它们是设计和实施CO2地质处置必须的工具［6］。可以评估预测地层的封存能力，评估CO2地质储存的可行性和可靠性，解释CO2在地层中的运移变化行为以及在注入到地层中后对其运移过程进行监测和分析。

在大多数通用多相流数值模拟软件中，相对渗透率和毛细压力的计算函数仅与当前饱和度有关，然而实验表明［7］，相对渗透率和毛细压力不仅和当前区域的饱和度有关，还和流体经历的过程及历史饱和度有关。尽管滞后效应特征曲线在石油工业已经被应用，但是在包气带、地热储藏工程、核废料地质处理以及CO2地质封存等这些两相流的问题中很少被应用到［8～9］。

目前在CO2地质封存方面，应用较为普遍的模拟软件为TOUGH2及2008年正式公开发布用于大规模并行计算的并行版本TOUGH2-MP［10］。目前成功应用最大的模型是东京湾的CO2地质处置模型，其单元数超过1 千万个［11］。TOUGH+ 是继 TOUGH2［12］之后的新一代TOUGH家族的模拟器，在继承了TOUGH2所有功能外，主要在数组处理效率方面提高模拟的计算效率，同时也提高了一些模拟功能。TOUGH+CO2是TOUGH+系列中模拟CO2封存的模块［13］，最新版本TOUGH+CO2中，在 iTOUGH2［14］中关于滞后效应模块的基础上改进更新了滞后效应模拟的功能。

1.3 TOUGH+CO2中相对渗透率相关理论

滞后效应特征曲线不仅取决于当前饱和度，还和历史饱和度以及经历的过程即疏干过程(非润湿相替代润湿相)或吸湿过程(润湿相替换非润湿相)相关［15～17］。

CO2深部咸水层封存中主要涉及多相流系统为H2O－CO2－NaCl系统，超临界状态CO2在TOUGH+CO2中表征为非水相，即气相，溶解于水中的CO2表征为溶液相；计算系统中流体相的相对渗透率的常用方法为 van Genuchten［18］公式。

考虑相对渗透率滞后现象时，TOUGH+CO2采用Parker和Lenhard根据无滞后效应van Genuchten方程修改得到的公式［14］，

式中:krl——液相相对渗透率，无量纲；

krg——气相相对渗透率，无量纲；

式中:Sl——液相饱和度，无量纲；

Slr——液相残余饱和度，无量纲；

Sls——饱和水饱和度，无量纲；

Sgr——残余气相饱和度，无量纲。

图2 历史饱和度对相对渗透率的影响Fig.2 Effect of the history saturation on relative permeability

1.4 TOUGH+CO2中毛细压力相关理论

在计算H2O－CO2系统中流体相的毛细压力时经常用到的方法为 van Genuchten［18］函数。

式中:Pcap——毛细压力(Pa)；

P0——进入毛细压力(Pa)。

考虑毛细压力滞后现象时采用基于van Genuchten方程修改得到的公式［14］，

式中:p——流体流动过程，疏干过程(d)或吸湿过程(w)；

mp——函数适配参数 mp=(nγ－1)/nγ，n、γ 为适配参数。

2 应用

2.1 模型基本参数

通过对案例模型的模拟，可以分析相对渗透率和毛细压力滞后效应对CO2咸水层地质封存的影响，为了研究方便，假设模型为理想条件模型。

进行CO2深部咸水层封存模拟时，需要给出的初始条件有压力、温度、盐度以及流体中CO2的质量分数。模型模拟的咸水层在地下1200m深，模拟范围为1km×1km，模型厚度为100m，模型上下边界为低渗透地层，其中盖层为5m厚泥岩，注入层为85m厚砂岩，下伏层为10m厚泥岩(图3)。模型共剖分成64000个网格，其中注入井附近网格进行加密处理。四周边界条件视为无流量边界，对于封闭系统，即储层被完全封闭，其CO2储存量由在储层中压力积聚引起孔隙体积的膨胀和咸水密度的增加来得到。

图3 模型范围XZ剖面示意图Fig.3 Profile of the model in the XZ direction

此模型中只考虑等温过程，初始温度为37℃，模型的顶部边界压力为1.20×107Pa，初始压力分布根据静水压力平衡计算得到。超临界状态下CO2的性质根据实验数据加入到程序，其他关键物理参数见表1。

表1 模型采用关键物理参数Table 1 Key properties of rocks in the model

方案1:不考虑滞后现象对CO2封存的影响，在计算中，相对渗透率和毛细压力的计算方式均选择van Genuchten函数。模型区域中介质的残余气相饱和度Sgr为0.177，整个含水层中的残余气相饱和度均为固定Sgr。在数值模拟中，根据固定值Sgr，计算相对渗透率特征曲线，在驱替过程中(包括疏干和吸湿过程)，相对渗透率均沿此曲线变化，同一流体饱和度下对应相同相对渗透率。毛细压力的计算亦是如此。

方案2:考虑滞后现象对CO2封存的影响。相对渗透率和毛细压力计算函数选择修改后的公式。模型区域中残余气相饱和度的最大可能值Sgrmax为0.177。在数值计算中，根据驱替过程转换点的气相饱和度和经历的过程计算相应的残余气相饱和度，然后通过插值得到扫描曲线，与主线不同。残余液相饱和度Slr为0.2，krgmax为 1。

CO2注入以0.926 kg/s(约3×104t/a)恒定速率注入4年，为观察停止注入后的运移状态，在第5～50年停止注入，模拟时间总共为50年。

2.2 结果分析

2.2.1 对封存量的影响

在模拟时间为50年时，99%的CO2存储于注入层中，只有1%存在盖层和底部地层中(表2)。因浮力作用CO2向上运移，由于盖层渗透率低，封闭性好，CO2大部分被封存在注入层中。所以两个方案中封存的CO2总量基本一致，但CO2存在相态比例不同。

表2 模拟时间为50年时各层CO2含量Table 2 Mass of CO2in various forms at 50 years (kg)

2.2.2 注入CO2的存在形式及其分布

图4～5为两种方案中封存CO2相态比例变化。气相是指整个地层中超临界状态CO2，包含可移动的自由态CO2和捕获封存不可移动气相CO2。随着模拟时间的增加，CO2和咸水层接触时间增加，气相CO2逐渐溶解于溶液中成为溶液相。

图4 无滞后效应模型整个地层中各相态CO2质量分数Fig.4 Mass fraction of CO2in the non-hysteretic case over the entire model

图5 滞后效应模型整个地层中各相态CO2质量分数Fig.5 Mass fraction of CO2in the hysteretic case over the entire model

模拟到50年时，无滞后方案中气相形态逐渐减少到70%，滞后方案中气相保持约45%的比例。产生这种差异的原因是无滞后方案中，残余气相饱和度为0.177，即含水层介质都默认捕获此数值残余气相饱和度，只有气相饱和度高于此残余气相饱和度时才运移。而在滞后方案中，是可能存在的最大气相饱和度，即每个网格元素中根据经历的过程不同，转换点饱和度不同，残余气相饱和度不同。因此滞后效应方案中捕获的残余气相总是小于等于无滞后效应方案中的残余气相，因此表现出在模拟时间到50年时，无滞后效应方案中封存的气相CO2比例大于滞后效应中的气相比例。同样，封存的CO2在含水层中的分布也不同，图6为在CO2注入阶段，即在疏干过程中气相CO2饱和度的分布。在疏干过程中，两种方案的相对渗透率特征曲线均沿主线变化，不涉及滞后现象，但因为残余气相饱和度的不同，相应特征曲线则不同。以注入点附近区域为例，无滞后效应的气相krg总小于滞后效应方案的krg，因此在晕图上表现出无滞后效应方案中气相CO2运移缓慢(图7)。

第4年末CO2停止注入后，在晕的尾端发生吸湿过程。如图7，无滞后方案捕获更多的残余气相，故CO2晕分布集中在运移过程中，在滞后方案中部分气相CO2则快速运移到盖层下，由于盖层的封闭性好，气相CO2晕逐渐向周围扩散。

图8表示数值计算中相对渗透率特征曲线。在无滞后方案中，特征曲线沿着疏干过程逆向变化，即相互驱替过程只有一条曲线。而在滞后方案中，相对渗透率和毛细压力不仅与当前区域的饱和度有关，还和经历的过程和历史饱和度有关。例如，在注入点网格区域，疏干过程中Sg增加，相应的krg增加，krl减少到0。在随后停止注入的吸湿过程中，Sg逐渐减少，krg减少，krl逐渐增大。但在同一Sg时，两个过程的相对渗透率不同，产生滞后的现象。同样，在中间某网格处，相对渗透率变化趋势相同，只是转换点不同，进而相对渗透率变化曲线不同。

图6 CO2注入过程中无滞后效应方案(上)和滞后效应方案(下)中气相饱和度分布图Fig.6 CO2plume evolution during the injection period in the non-hysteretic case(top)and the hysteretic case(bottom)

图7 停止注入后无滞后效应方案(上)和滞后效应(下)方案中气相饱和度随时间变化分布图Fig.7 CO2plume evolution during the post-injection period in the non-hysteretic case(top)and the hysteretic case(bottom)

图8 模型中注入点网格相对渗透率变化曲线Fig.8 Relative permeability path for the injection grid

图9～10对应两种方案中毛细压力特征曲线，选取模型中注入点网格和中间网格两个不同位置得到其毛细特征曲线。滞后方案为一次完整疏干和吸湿过程，在中间网格处，气相饱和度从0逐渐增大到转换点饱和度，然后沿扫描曲线减少到相应残余气相饱和度，并未沿疏干过程主线变化，导致在同一气相饱和度下毛细压力不同。无滞后模型中，毛细压力均沿主线变化。在实际情况中，仅用实验测试得到的残余气相饱和度赋值给含水层来进行数值模拟，在封存能力、安全评估等方面会造成误差。

图9 滞后效应模型中不同位置毛细压力变化曲线Fig.9 Capillary pressure path for grids in the hysteretic

图10 无滞后效应模型不同位置毛细压力变化曲线Fig.10 Capillary pressure path for grids in the non-hysteretic case

3 结论

在CO2咸水层封存的数值模拟计算中，相对渗透率和毛细压力计算函数是重要计算参数之一。通过设计对应无滞后效应和考虑相对渗透率、毛细压力滞后的两种方案，运用数值模拟软件TOUGH+CO2进行模拟，结果表明不同相对渗透率和毛细压力计算函数的选择表现出较大差异。滞后效应方案能更详细描述每个网格中的相对渗透率和毛细压力变化情况。在盖层封闭性较好时，滞后效应对CO2封存总量没有影响，但相态比例不同。

不同相对渗透率和毛细压力计算函数的选择表现出CO2晕运移分布不同，对实际封存方案设计，封存能力评估以及环境风险评价等都会带来不同影响，因此在CO2咸水层封存数值模拟中，尤其是停止注入之后应引入与滞后现象相关的模拟计算。

［1］郝艳军，杨顶辉.二氧化碳地质封存问题和地震监测研究进展［J］.地球物理学进展，2012，27(6):2369 － 2383.［HAO Y J，YANG D H.Research progress of carbon dioxide capture and geological sequestration problem and seismic monitoring research［J］.Progress in Geophys，2012，27(6):2369－2383.(in Chinese)］

［2］张炜，吕鹏.二氧化碳地质封存中对流混合过程的研究进展［J］.水文地质工程地质，2013，40(2):101－107.［ZHANG W，LV P.Density-driven convection in carbon dioxide geological storage:a review ［J］.Hydrogeology ＆ Engineering Geology，2013，40(2):101－107.(in Chinese)］

［3］Juanes R，Spiteri E J，Orr F M，et al.Impact of relative permeability hysteresis on geological CO2storage［J］.Water Resource Research，2006，42(12):1－13.

［4］Doughty C.Investigation of CO2plume behavior for a large-scale pilot test of geologic carbon storage in a saline formation ［J］.Transport in porous media，2010，82(1):49 －76.

［5］许雅琴，张可霓，王洋.基于数值模拟探讨提高咸水层CO2封存注入率的途径［J］.岩土力学，2012，33(12):3825－3832.［XU Y Q，ZHANG K N，WANG Y.Numerical investigation for enhancing injectivity of CO2storage in saline aquifers［J］.Rock and Soil Mechanics，2012，33(12):3825 －3832.(in Chinese)］

［6］张晓宇，成建梅，刘军，等.CO2地质处置研究进展［J］.水文地质工程地质，2006，33(4):85－89.［ZHANG X Y，CHENG J M，LIU J，et al.An overview of underground sequestration of carbon dioxide［J］.Hydrogeology ＆ Engineering Geology，2006，33(4):85 －89.(in Chinese)］

［7］周琦，朱学谦，刘传喜.毛细管力滞后曲线和相渗透率滞后曲线及处理［J］.油气采收率技术，1999，6(3):47－50.［ZHOU Q，ZHU X Q，LIU C X.Curves of hysteretic phenomenon in relative permeability and capillary pressure［J］.Oil＆ Gas Recovery Technology，1999，6(3):47－50.(in Chinese)］

［8］Finsterle S.Multiphase inverse modeling:Review and iTOUGH2 applications ［J］.Vadose Zone Journal，2004，3(3):747－762.

［9］Zhang Y Q，Oldenburg C M，Finsterle S，et al.Systemlevel modeling for economic evaluation of geological CO2storage in gas reservoirs［J］.Energy Conversion and Management，2007，48(6):1827 －1833.

［10］Zhang K，Doughty C，Wu Y S，et al.Efficient parallel simulation of CO2geologic sequestration in saline aquifers［C］//Lawrence Berkeley National Laboratory.SPE Reservoir Simulation Symposium.Berkeley:Society of Petroleum Engineers，2007.

［11］Yamamotoa H，Zhang K，Karasakib K，et al.Largescale numerical simulation of CO2geologic storage and its impact on regional groundwater flow:A hypothetical case study at Tokyo Bay［J］.Japan Energy Procedia，2009，1(1):1871－1878.

［12］Pruess K，Oldenburg C，Moridis G.TOUGH2 User’s Guide［R］. LBNL-43134.Berkeley:Lawrence Berkeley National Laboratory，1999.

［13］Zhang K，Moridis G，Pruess K.TOUGH+CO2:A multiphase fluid-flow simulator for CO2geologic sequestration in saline aquifers［J］.Computers ＆Geosciences，2011，37(6):714 －723.

［14］Finsterle S.iTOUGH2 user’s guide［R］.LBNL－40040.Berkeley:Lawrence Berkeley National Laboratory，1999.

［15］Doughty C. Modeling geologic storage of carbon dioxide Comparison of non-hysteretic and hysteretic characteristic curves［J］.Energy Conversion and Management，2007，48(6):1768 －1781.

［16］王大纯，张人权，史毅虹，等.水文地质学基础［M］.北京:地质出版社，1995.

［17］薛禹群.地下水动力学［M］.北京:地质出版社，1997.

［18］Van Genuchten M Th，A closed-form equation for predicting the hydraulic conductivity of unsaturated soils［J］.Soil Science Society of America Journal，1980，44(5):892－898.