封存CO2裹挟可挥发性污染物的迁移模型分析

杨勇，刘永忠,2，于博，丁天

（1西安交通大学化学工程系，陕西 西安 710049；2热流科学与工程教育部重点实验室，陕西 西安 710049）

封存CO2裹挟可挥发性污染物的迁移模型分析

杨勇1，刘永忠1,2，于博1，丁天1

（1西安交通大学化学工程系，陕西 西安 710049；2热流科学与工程教育部重点实验室，陕西 西安 710049）

为了研究封存CO2注入过程可挥发性污染物的释放和迁移特性，本文建立了CO2-水-残余油多相流驱替过程中的污染物迁移模型。采用数学模拟方法分析了多相驱替过程和污染物迁移过程，并研究了相间传质特性、初始油相分布和CO2注入速率等对污染物迁移过程的影响。研究表明：CO2驱替过程将促使可挥发污染物进入CO2相并随之在地层中迁移，逐渐形成相间传质区域。相间传质区域的演化反映了污染物释放和CO2裹挟污染物迁移的特性。可挥发污染物的传质系数越大，相间传质区域越窄，油相饱和度衰减越快；油相初始饱和度较大时，其饱和度衰减相对缓慢，对应的相间传质区域也较窄。当CO2注入速率增大时，相间传质区域增大，油相饱和度衰减变快。本文模型可用于不同地质储层环境下封存CO2时可挥发性污染物的迁移特性分析，并用于封存CO2的风险分析与评价。

CO2；封存；多相流；驱替过程；相间传质；污染物迁移

引 言

将CO2注入废弃的油气藏等地质储层中是缓解温室效应的重要技术之一[1-2]。这些储层中通常充满了二次或三次采油采气后残余的油相和水相。当CO2注入这些地层时，超临界CO2将驱替宿主流体，并裹挟残余油相中的可挥发污染物进入CO2相。随后，可挥发污染物将随注入的CO2在地下环境中迁移，将可能通过贯穿地层的开采井、废弃井[3]或地质缺陷及断裂[4]向饮用水层及地表迁移，对人类生态环境造成潜在的危害。显然，超临界态CO2注入地层引起污染物的释放和迁移增大了CO2封存的不确定性和风险，CO2裹挟污染物的迁移和泄漏将加剧对生态环境的危害。因此，封存CO2裹挟地下污染物迁移的问题对于CO2地下封存风险评估具有重要意义。

封存CO2裹挟污染物迁移问题需要同时考虑多相流动和污染物相间传质过程。目前，人们对CO2多相驱替中污染物的释放和迁移特性缺乏直观和深入的理解，对CO2裹挟污染物迁移的特性也缺乏定量描述。然而，研究人员对封存CO2在地层中多相流动问题已有许多研究工作。Class等[3]提出了CO2地质封存过程的若干基准问题，给出了基本的物理模型以便于不同数值方法的比较和验证。Liu等[5-6]采用清晰和非清晰界面追踪方法描述了CO2注入到盐水层中的封存和泄漏问题，特别分析了两相界面的演变和CO2的泄漏特性。Wei等[7-8]研究了CO2在煤层中封存的特性，探讨了注入的CO2与煤层和煤气的相互作用。Sohrabi等[9-10]分析了CO2注入到油藏中的多相驱替过程，讨论了CO2封存和石油开采回收的关系。与此相应，地下环境中污染物的控制和修复技术也引起广泛关注[11]。Al-Raoush[12]研究了地层孔隙中非水溶剂相微团的分布特征。Powers等[13]指出合理地描述地下污染物源的溶解释放过程需采用非平衡传质模型。Rivett等[14]系统研究了可挥发有机物在地下污染物源中的分布特征和迁移规律，总括了污染物溶解、挥发和迁移的相关模型。Nguyen等[15]通过三维数值模型描述了气相抽提法修复土壤环境中污染物的传质过程，给出了污染物随气相萃取过程的迁移特征。这些污染物分析和模型描述为地下环境污染物的迁移提供了分析依据。

本文根据地下储层孔隙中污染物分布特征，在多相流模型基础上建立CO2裹挟污染物迁移的模型。借助模型来分析污染物释放和CO2裹挟污染物迁移特性，并研究了相间传质特性和CO2注入速率等因素对污染物迁移规律的影响，以期深入认识多相流动过程中的污染物迁移特性，为CO2封存和泄漏风险控制提供理论分析基础。

1 封存CO2裹挟污染物的迁移模型

1.1 物理模型

参考CO2封存的基准物理模型[3]，目标地质储层的上下边界为冠岩层所限定，储层中初始由水相和残余油相所饱和，CO2以一定速度通过注入井进入到地质储层中进行封存。多相驱替过程以及残余油相中污染物的迁移过程如图 1所示。注入的CO2相驱替储层中宿主流体进入地层孔隙结构中。根据Al-Raoush等[12-13]对不同非水溶剂相微团的表征，残余油相可能以单粒、柱状、大球和指形等4种典型形状分布在储层中。随着多相驱替过程进行，CO2相与残余油相相互接触诱使残余油相中的可挥发性污染物释放进入CO2相。

实际上，在多相流动体系中，两相或多相存在的区域发生相间传质，随之各相的组成也会发生变化。本文目的是研究油相中的可挥发性污染物释放进入CO2相并随之迁移的过程。地质储层中初始均匀分布着残余油相和水相，可认为污染物在水相和固体骨架上的传质过程始终处于平衡状态。而超临界CO2与水相之间溶解作用和矿物化学反应是长期封存的主要机理，在本文中注入短期过程暂不考虑这些作用对污染物迁移的影响。

1.2 数学模型

1.2.1 多相驱替过程的控制方程 在上述CO2-水-油多相体系中，各相的流动满足质量守恒方程和Darcy定律，即

图1 封存CO2裹挟可挥发污染物的迁移过程Fig.1 Migration of volatile contaminant during supercritical CO2injection

式中，φ为地层的孔隙率；S是饱和度；ρ是密度，kg·m−3；t为时间，s；Qα是源项，kg·m−3·s−1；κ是本征渗透率，m2；rkα是相对渗透率；μ是动力黏度，Pa·s；p是相压力，Pa；g为重力加速度，m·s−2；z是高度，m；v是各相速度，m·s−1。下角标α代表CO2相、水相和油相。残余油相保持停滞状态，即o0=v。

各相的饱和度满足

本文采用Brooks-Corey模型[6,16]来描述各相间的毛细压力、相对渗透率和饱和度的关系，使上述控制方程组封闭。

1.2.2 相间传质过程 本文采用非平衡界面传质模型描述污染物的相间传质的动力学特征。在此仅考虑可挥发性污染物向CO2相的释放过程，则式（1）中，。假定油相是单一组分的污染物源相，其进入主体流体相的相间传质通量为

式中，fk表示相间传质系数，m·s−1；goa表示C O2-油相两相接触的比界面面积，m−1；分别是可挥发性污染物在CO2相中的浓度和平衡浓度，kg·m−3。

比界面面积取决于孔隙结构和两相接触特征等因素[13,17]，可表示为

其中，(A/V)b是油相微团的比面积，取决于油相在地层孔隙中的几何分布特征；f是油相微团向着流体相的曝露表面积分数，其取值范围在1%～50%，对当前CO2-水-油相三相体系可表示为

式（4）、式（5）和式（6）描述了相间传质过程受到界面传质速率系数、地层孔隙结构、油相微团结构和各相饱和度的影响。

1.2.3 污染物迁移方程 可挥发性污染物在CO2相中的对流-弥散方程可表示为

式中，gD为污染物的弥散张量[13,18]；I代表单位矩阵；是污染物在CO2相的自由分子扩散系数，m2·s−1；τ是多孔介质的曲折因子；La和Ta表示纵向和横向的弥散系数，m。

1.3 初始条件和边界条件

在CO2注入之前，图 1所示的地层包含了水相和作为污染物源的残余油相，且水相压力为静水压力。CO2沿注入井进入地层，地层上下为冠岩层，即对各相均不可渗透。地层远端为压力边界。则上述模型的初始条件表示为

边界条件为

式中，topp是储层顶部的静水压力，Pa；H为目标储层的高度，m；ginjv是注入速度，取决于CO2的注入通量。

2 模型参数与验证

2.1 模型参数

本文模型采用的地层特征参数数据和操作参数如表 1所示[3,19-21]。表1中的地质条件来源于CO2封存相关基准问题的物理模型[3]，该基准问题的物理模型广泛应用于数值模型比较和CO2封存过程基本特征分析。表 1中污染物相关的特征参数包括油相饱和度、密度、污染物的溶解度和传质系数。残余油相的饱和度取决于地下储层中在孔隙结构中体积分数，通常在10%～35%。而残余油相的组成成分复杂，主要由重质石油烃和多环芳烃构成。在可挥发性石油类化合物中，苯系物和多环芳烃被认为是危害最大的，而苯系物和双环芳烃通常包含在石油组成的C6～C11范围中。双环芳烃萘（C8）的挥发性低于苯系物的挥发性，但略高于其他重质石油烃，并且分子量也较为接近其他重质烃。因此，为了简化本文的分析模型，采用双环芳烃萘（C8）作为模型示例污染物，污染物的相关特性采用双环芳烃萘的物性参数[19-21]。

本文采用COMSOL Multiphysics 3.5aTM软件对数学模型进行求解，以获得污染物的迁移特性。

表1 地层条件和特征参数Table 1 Formation and operational data

2.2 模型验证

在地质封存条件下，地下污染物迁移的实测和现场数据非常匮乏。为了验证本文CO2裹挟污染物迁移模型，将残余油相的初始饱和度置零，模型将简化为CO2-水相的两相流模型。基于CO2封存的基准问题物理模型[3,24]，可验证本文的两相流模型。在考虑毛细压力作用条件下，不同数值模型模拟的CO2相饱和度分布结果比较如图 2所示。其中CO2相分布对应的是相似变量r2t−1的计算结果。图 2表明CO2相饱和度的分布保持一致。由此，将本文模型拓展用于分析CO2-水-油相体系的可挥发性污染物的多相流动特征是可行可靠的。

图2 多相流模型验证和比较Fig.2 Comparison of multiphase flow model

3 封存CO2裹挟污染物迁移特性分析

3.1 相间传质区域

图 3给出CO2注入360 d后各相的分布和CO2中污染物浓度分布，油相初始残余饱和度Soin为0.2。

由图可见，随着超临界CO2的持续注入，在多相驱替和重力的共同作用下，CO2相逐渐上浮。入口端的油相饱和度则随着可挥发污染物进入CO2相而逐渐衰减到残余饱和度值。可挥发污染物在CO2相中也逐渐累积，并在CO2相前沿区域达到其平衡浓度。

图 3中，通过比较CO2相、油相饱和度和污染物浓度分布可见，在CO2相覆盖的区域中，相间传质仅发生在局部区域，该区域称为相间传质区域（interphase mass transfer region，IMTR）。靠近注入口端，油相饱和度达到其残余饱和度值，即油相中的可挥发性污染物耗尽或受孔隙结构限制不能进入CO2相。在CO2相分布前沿，污染物在其中逐渐累积达到了其平衡浓度，即传质推动力减为零。当前CO2注入360 d的条件下，仅在25～100 m的局部范围内发生相间传质。

图 4给出目标储层中部H=15 m处模拟计算结果随着相似变量（r2t−1）的变化特性。不同注入时间下，模型计算结果随相似变量的变化特征具有一致性，即在任意时刻下，CO2饱和度、油相饱和度和污染物浓度都具有与图 4相似的分布特征。显然，在图 4中，可将CO2相覆盖范围分成3个区域：（1）纯CO2相区域，这一区域靠近入口，其中的CO2相中不含可挥发污染物，并且作为污染物源的油相达到其残余饱和度值；（2）相间传质区域，可挥发污染物进入CO2相的中间区域；（3）可挥发性污染物饱和区域，该区域靠近CO2相前沿，其中可挥发污染物在CO2相中累积到其饱和浓度，而油相保持初始条件下的分布特征。随着CO2的持续注入，相间传质区域将随着CO2的流动而向前推移。因此，相间传质区域的变化反映了污染物释放特征和CO2裹挟污染物的迁移特性。

图4 注入360 d时各相饱和度和污染物浓度随着相似变量的分布Fig.4 Phase saturation and contaminant concentrations as a function of similarity variable on 360 d

3.2 相间传质系数对相间传质区域的影响

相间传质速率通常取决于污染物组分的挥发性、孔隙中污染物源的几何形态以及多相流动特性。不同体系传质系数的取值差异也很大，如多孔介质中甲苯传质系数的范围为10−7～10−6m·s−1[25]。为了分析相间传质系数对传质区域的影响，本文选取相间传质系数的若干数量级来分析其对可挥发性污染物迁移特性的影响。当CO2注入360 d后，不同传质系数对应的污染物迁移特性如图 5 所示。由图5可见，当传质速率系数kf值的数量级从10−6m·s−1减小到10−10m·s−1时，未达到污染物平衡浓度的区域增大，传质区域变大，油相的衰减速度也逐渐变慢。

3.3 油相初始饱和度对相间传质区域的影响

残余油相的初始饱和度代表着地层孔隙中的污染物源的数量。图 6给出不同初始饱和度条件下油相饱和度和CO2相中可挥发污染物浓度的变化特性。为了便于比较，本文对油相饱和度作均一化处理，即采用(So−Sor)/(Soin−Sor)表征油相饱和度变化。由图可见，较高的油相初始饱和度对应的界面传质区域相对较窄，即整体传质速率较快，但可挥发性污染物含量的增大也意味着污染物源的衰减速度降低。此外，初始油相饱和度较大时，CO2相中污染物曾达到其平衡浓度的区域也扩大。

3.4 CO2注入速率对相间传质区域的影响

CO2注入流率决定了不同位置的流体流速和流体饱和度分布，从而影响了CO2相和油相的相互接触和相间传质区域的推进。图 7给出了不同CO2注入速率下，油相饱和度和污染物浓度的分布特性。注入速率不同，CO2的前沿位置不同，鉴于柱坐标下的质量守恒条件，污染物浓度采用量纲1距离或相似变量对结果进行分析[图 7（b）]。显然，不同的CO2注入速率对应污染物迁移特性具有相似性，即油相饱和度随量纲1距离的变化一致。然而，注入流量增加时相间传质区域将增大，加速油相的衰减。进入CO2相中可挥发污染物的绝对量将增加。图7（b）也表明，CO2注入流量的增大会扩大界面传质区域和达到残余油相饱和度的区域。

图5 CO2注入360 d不同传质系数对应的油相和污染物浓度分布Fig.5 Effect of interphase mass transfer coefficient on distribution of oil phase and contaminant on 360 d

图6 油相初始饱和度对相间传质的影响Fig.6 Effect of initial oil saturation on interphase mass transfer region on 360 d

4 结 论

CO2注入废弃油气藏等地质储层中是缓解温室效应的重要技术之一。然而，注入CO2裹挟污染物在地下环境中迁移将可能对人类生态环境造成潜在的危害。为了研究封存CO2裹挟地下污染物的迁移过程，本文建立了超临界CO2注入水相和油相所饱和的地质储层中引起的可挥发性污染物迁移模型。该模型可以有效地描述孔隙结构中两相驱替过程的污染物释放（溶解）和迁移特性。数学模拟计算结果表明：

（1）CO2注入过程中，可挥发性污染物的释放和迁移将CO2所覆盖的区域分为纯CO2相区域、相间传质区域和污染物饱和区域等3个区域。

（2）相间传质区域随着CO2的持续注入和油相衰减向前稳定推移。相间传质区域的变化表明了污染物释放特性和CO2裹挟污染物迁移的演化特性。

（3）较大的传质系数缩小了相间传质区域的范围，但加快了油相污染物源的衰减；而当油相初始饱和度较高时，也会导致相间传质区域较窄，但对应油相污染物源的衰减速度也较缓慢；较大的CO2注入流量将增大相间传质区域，并增加了油相污染物源的衰减速度。

本文模型可用于分析不同地质储层环境下可挥发性污染物随CO2注入的迁移特性和污染物的演化分布特性，为CO2地下封存的风险分析提供基础模型。需要指出的是，在该模型的基础上，进一步考虑多相流驱替过程的非等温效应、污染物溶解度和相间传质的动态变化特性、CO2与宿主流体之间的溶解以及生物、矿物化学反应等作用将可使模型分析更接近实际。这些工作都值得进一步深入研究。

图7 CO2注入速率对相间传质区域的影响Fig.7 Effect of injection rate on interphase mass transfer region

[1] Zhao X,Liao X. Evaluation method of CO2sequestration and enhanced oil recovery in an oil reservoir, as applied to the Changqing oilfields, China [J].Energy Fuels, 2012, 26 (8): 5350-5354.

[2] Mohammad S A,Gasem K A M. Multiphase analysis for high-pressure adsorption of CO2/water mixtures on wet coals [J].Energy Fuels, 2012, 26 (6): 3470-3480.

[3] Class H, Ebigbo A, Helmig R,et al. A benchmark study on problems related to CO2storage in geologic formations [J].Computat. Geosci., 2009, 13 (4): 409-434.

[4] Pruess K. Leakage of CO2from geologic storage: role of secondary accumulation at shallow depth [J].Int. J. Greenh. Gas Con., 2008, 2 (1): 37-46.

[5] Liu Y Z, Wang L,Yu B. Sharp front capturing method for carbon dioxide plume propagation during injection into a deep confined aquifer [J].Energy Fuels,2010, 24: 1431-1440.

[6] Liu Yongzhong (刘永忠), Wang Le (王乐), Zhang Jialiu (张甲六). Prediction on leakage and rate characteristics of CO2storage [J].CIESC Journal(化工学报), 2012, 63 (4): 1226-1233.

[7] Wei X, Massarotto P, Wang G,et al. CO2sequestration in coals and enhanced coalbed methane recovery: new numerical approach [J].Fuel,2010, 89 (5): 1110-1118.

[8] Chen D, Pan Z, Liu J,et al. Modeling and simulation of moisture effect on gas storage and transport in coal seams [J].Energy Fuels,2012, 26 (3): 1695-1706.

[9] Sohrabi M, Kechut N I, Riazi M,et al. Coreflooding studies to investigate the potential of carbonated water injection as an injection strategy for improved oil recovery and CO2storage [J].Transp. Porous Med., 2011, 91 (1): 101-121.

[10] Jahangiri H R,Zhang D. Ensemble based co-optimization of carbon dioxide sequestration and enhanced oil recovery [J].Int. J. Greenh. Gas Con., 2012, 8: 22-33.

[11] Yang Mingxing (杨明星), Ye Xueyan (冶雪艳), Yang Yuesuo (杨悦锁),et al. Residual and migration diversities of petroleum hydrocarbons in vadose zone [J].CIESC Journal(化工学报), 2013, 64 (11): 4231-4239.

[12] Al-Raoush R I. Impact of wettability on pore-scale characteristics of residual nonaqueous phase liquids [J].Environ. Sci. Technol., 2009, 43 (13): 4796-4801.

[13] Powers S E, Loureiro C O, Abriola L M,et al. Theoretical study of the significance of nonequilibrium dissolution of nonaqueous phase liquids in subsurface systems [J].Water Resour. Res., 1991, 27 (4): 463-477.

[14] Rivett M O, Wealthall G P, Dearden R A,et al. Review of unsaturated-zone transport and attenuation of volatile organic compound (VOC) plumes leached from shallow source zones [J].J. Contam. Hydrol., 2011, 123 (3/4): 130-156.

[15] Nguyen V T, Zhao L,Zytner R G. Three-dimensional numerical model for soil vapor extraction [J].J. Contam. Hydrol., 2013, 147: 82-95.

[16] Kolditz O, Bauer S, Beyer C,et al. A systematic benchmarking approach for geologic CO2injection and storage [J].Environmental Earth Sciences, 2012, 67 (2): 613-632.

[17] Niessner J, Hassanizadeh S M. Modeling kinetic interphase mass transfer for two-phase flow in porous media including fluid-fluid interfacial area [J].Transp. Porous Med., 2009, 80 (2): 329-344.

[18] Yang L, Zhang X. Modeling of contaminant transport in underground coal gasification [J].Energy Fuels,2009, 23: 193-201.

[19] McHugh M, Paulaitis M E. Solid solubilities of naphthalene and biphenyl in supercritical carbon dioxide [J].Journal of Chemical & Engineering Data, 1980, 25 (4): 326-329.

[20] Higashi H, Iwai Y, Arai Y. Solubilities and diffusion coefficients of high boiling compounds in supercritical carbon dioxide [J].Chem. Eng. Sci., 2001, 56: 3027-3044.

[21] Owabor C N, Ogbeide S E,Susu A A. Estimation of transport and degradation parameters for naphthalene and anthracene: influence of mass transfer on kinetics [J].Environ. Monit. Assess., 2009, 169 (1-4): 607-617.

[22] Panagiotopoulos A Z, Reid R C. New mixing rule for cubic equations of state for highly polar, asymmetric systems//Chao K C, Robinson R L. ACS Symposium Series, Vol 300. Equations of State-Theories and Applications[C]. American Chemical Society, 1986: 571-582.

[23] Jossi J A, Stiel L I,Thodos G. The viscosity of pure substances in the dense gaseous and liquid phases [J].AIChE J., 1962, 8 (1): 59-63.

[24] Pruess K, Garcia J, Kovscek T,et al. Code intercomparison builds confidence in numerical simulation models for geologic disposal of CO2[J].Energy,2004, 29 (9/10): 1431-1444.

[25] Anwar A H M F, Tien T H, Inoue Y,et al. Mass transfer correlation for nonaqueous phase liquid volatilization in porous media [J].Environ. Sci. Technol., 2003, 37 (7): 1277-1283.

Analysis of migration model for volatile contaminants induced by CO2geological storage

YANG Yong1, LIU Yongzhong1,2, YU Bo1, DING Tian1
(1Department of Chemical Engineering,Xi’an Jiaotong University,Xi’an710049,Shaanxi,China;2Key Laboratory of Thermo-Fluid Science and Engineering,Ministry of Education,Xi’an710049,Shaanxi,China)

In order to investigate characteristics of volatile contaminant release and migration during CO2geological sequestration, a contaminant migration model with multiphase displacement process in a CO2-aqueous phase-residual oil phase system was proposed. Numerical simulations were performed to figure out the multiphase displacement and the contaminant migration, and further to explore the effects of interphase mass transfer, initial profiles of the residual oil and injection rates on the volatile contaminant migration. The results indicate that the volatile contaminant enters into the supercritical CO2phase due to the multiphase displacement process, and migrates with the mobile phase in the subsurface formation. An interphase mass transfer region (IMTR) gradually forms. The evolution of IMTR directly reflects the release characteristics of contaminant and migration behaviors of CO2coerced with contaminant. A larger interphase mass transfer coefficient results in a narrower IMTR and a faster decay rate of the oil phase. When the initial saturation of the oil phase is larger, the IMTR becomes much narrower and the oil phase decays much slower. With increasing CO2injection rate, IMTR increases and the oil phase decays rapidly. The proposed model and the results intuitively describe the migration characteristics of thevolatile contaminants in a variety of geological reservoirs during CO2geological sequestration. The proposed model can also be used for the risk analysis and evaluations of CO2storage.

CO2; sequestration; multiphase flow; displacement process; interphase mass transfer; contaminant migration

LIU Yongzhong, yzliu@mail.xjtu.edu.cn

10.11949/j.issn.0438-1157.20150069

TQ 021.1

：A

：0438—1157（2015）10—4155—08

2015-01-16收到初稿，2015-03-28收到修改稿。

联系人：刘永忠。

：杨勇（1986—），男，博士研究生。

国家自然科学基金项目（21176198）；教育部博士点基金项目（20120201110071）。

Received date: 2015-01-16.

Foundation item: supported by the National Natural Science Foundation of China (21176198) and Doctoral Fund of Ministry of Education of China (20120201110071).