基于TOUGHREACT-MP的苏北盆地盐城组咸水层CO2矿物封存数值模拟

莫绍星, 龙星皎, 李 瀛, 郑 菲, 施小清, 张可霓, 赵 良

1.表生地球化学教育部重点实验室/南京大学地球科学与工程学院, 南京 210023 2.北京师范大学水科学研究院, 北京 100875

基于TOUGHREACT-MP的苏北盆地盐城组咸水层CO2矿物封存数值模拟

莫绍星1, 龙星皎1, 李 瀛1, 郑 菲1, 施小清1, 张可霓2, 赵 良1

1.表生地球化学教育部重点实验室/南京大学地球科学与工程学院, 南京 210023 2.北京师范大学水科学研究院, 北京 100875

基于TOUGHREACT并行版，建立了苏北盆地盐城组下段砂岩层在碳封存中的CO2-水-岩反应二维径向模型，并评估储层的矿物封存潜力，探讨分析了网格剖分精度和储层各向异性对矿物封存模拟的影响。模拟结果表明：在CO2封存过程中，片钠铝石、方解石和菱铁矿沉淀较明显，在CO2注入5 000 a后矿物封存的比例高达34.0%；网格剖分精度的不同对矿物封存反应路径没有影响，但粗网格模型计算得到的矿物封存量偏高；降低储层的kv(垂向渗透率)在短期内会促进CO2的水平运移，有利于溶解和矿物封存；但随着时间延长，降低kv对对流作用的抑制开始体现，表现为1 000 a后中等kv值的模型计算出的矿物封存量高于较大kv和较小kv值的模型。

CO2的矿物封存；网格剖分精度；各向异性；苏北盆地

0 引言

CO2的地质封存是减少碳排放最有效的措施之一，地下深部咸水层因其分布广泛、人类尚不能利用而成为最有前景的CO2地质封存的目标层[1]。目前国外已开展了广泛的研究，取得了显著成果[1]。虽然我国的地质封存潜力非常可观[2]，但是相关研究还比较滞后，迫切需要加强。

在咸水层封存中，CO2的封存形式有构造封存、残余气体封存、溶解封存和矿物封存，矿物封存是最理想的封存形式[1]。由于CO2咸水层封存涉及的物理和化学过程非常缓慢，数值模拟方法是重要的研究手段，许多模拟软件已成功运用于相关研究中[3-5]。在数值模拟中，空间离散会引起误差，这些误差与网格剖分精度有关[6]。例如当网格剖分太粗时，可能会高估CO2的迁移距离[7]，难以精确刻画CO2的迁移过程[8-9]，模型短期内会由于数值弥散而高估了CO2的溶解量[10-12]，在长期中则由于抑制了对流作用而低估CO2的溶解量[11-13]。目前对网格剖分精度的影响研究较少。Audigane等[13]的研究表明，提高剖分精度对模拟结果的影响较小；但由于受计算量巨大所限，他们采用的办法是只增加注入井附近区域的剖分精度，而降低远离注入井区域的剖分精度。

注入到储层中的CO2在浮力作用下将向上运移，当遇到低渗透性盖层后将向四周运移[7, 13]。在各向异性储层中，kv(垂向渗透率)的降低，将增加CO2向上运移的阻力，CO2将不得不沿高渗透层水平运移到储层的中部和底部。显然，这增加了CO2与咸水的接触时间和接触面积，促进了CO2溶解，进而促进了矿物的溶解和新矿物沉淀[14-15]。kv降低又会抑制对流作用，主要是因为：储层上部的CO2径向运移距离减小，减小了对流发生的区域[16-17]；垂向渗透速率减小[16]；更多的CO2运移到储层的中部和底部后，下部咸水的密度也会增加，使上下密度差减小[18-20]。由此可见，改变kv对CO2的水平运移和系统的对流作用会产生相反的效果(促进或抑制)，这使各向异性对封存机制的影响变得复杂。但目前在研究各向异性的影响时，大都只分析了各项异性对其中一种作用(对流作用或CO2的水平运移)的影响，并且分析的重点是对CO2溶解的影响，对矿物封存的影响研究不够。在Zhang等[21]的研究中，当降低kv时，CO2水平运移增强，促进了溶解和矿物封存，但该研究未说明降低kv对对流作用的影响。Jahangiri等[14]、Han等[17]和Ukaegbu等[22]也发现kv/kh(kh指水平渗透率)值越小，CO2溶解越多，但他们均未考虑CO2的矿物封存机制。Audigane等[13]和Lindeberg等[16]的研究均表明:对流作用对CO2溶解的促进作用在系统的流体力学不稳定性足够大时才比较明显；而kv越小，不稳定性增加得越慢[16, 18, 23]。Kihm等[24]的研究则表明，对流作用首先促进溶解封存进而促进矿物封存。综上所述，有必要系统地研究储层各向异性对评估CO2的封存机制的影响。

苏北盆地是含油气沉积盆地，其深部存在体积巨大的咸水层，具有成为CO2地质封存场所的地质条件。本研究利用TOUGHREACT并行版软件建立了二维径向模型，针对苏北盆地盐城组下段(N1y)砂岩层内的长期CO2-水-岩作用，分析了阳离子的反应路径及CO2的矿物封存潜力，评估了该岩层作为碳封存储层的可行性，并探讨网格剖分精度和储层的各向异性对评估矿物封存量的影响。

1 模拟软件简介

TOUGHREACT是基于非等温多相流体模拟软件TOUGH2的框架，将水流、溶质运移和地球化学反应耦合形成的一套可变饱和地质介质中非等温多相流体地球化学反应运移模拟软件[25]，其中的ECO2N模块已被广泛运用于CO2地质封存的研究中[13, 21, 25-26]。TOUGHREACT-MP是TOUGHREACT的并行版本，采用类似于TOUGH2-MP[27]的并行计算方案，引进区域分解算法对单机版进行了粗粒并行化，从而显著提高了计算规模和效率。

2 模型的建立

2.1 区域地质条件和模型概化

苏北盆地内部发育一系列北东向的南陡北缓的断陷和凹陷，是油气聚集的良好场所，也为CO2地质封存提供了良好的储存空间。盐城组为河流环境的沉积产物，其下段(N1y)由3个沉积旋回组成，以粉砂质泥岩与中粗砂岩、细砾岩互层为特征，顶部埋深一般大于1 000 m，孔隙度为23%～32%，渗透率为(0.05～1.10)×10-12m2；上段(N2y)为泥岩段，可作为区域盖层。从埋深、岩层特征、区域盖层封堵性、地质圈闭条件等因素判断，盐城组下段是潜在的CO2封存目标层[28-29]。

在基础模型(模型A)中，将储层概化为径向100 km、垂向40 m的圆柱形二维模型，储层顶板的标高为-1 000 m。径向上不等距递增剖分了111层，最外一层体积设为无限大，视为一类边界。垂向上等距剖分40层，总共剖分4 440个网格。此外，假设储层上下边界均为隔离边界(不参与反应)。注入井位于模拟区中心，井孔位于-1 032～-1 040 m，以4 kg/s的恒定速率注入超临界CO2，持续时间为10 a，总模拟时间为5 000 a。

为了研究网格剖分精度对评估CO2矿物封存量的影响，建立一个具更粗网格的模型(模型B)和一个具更精细网格的模型(模型C)，模型B在径向和垂向上的分层数较模型A均减少一半，模型C则分别增加一倍；为了研究储层的各向异性对评估CO2矿物封存量的影响，在模型A的基础上另建立2个模型，分别将kv降低1个(模型D)和2个(模型E)数量级(表1)。

2.2 水文地质参数值和储层物质组分的确定

表1 模型方案及研究目的

盐城组砂岩储层的矿物成分见文献[32]，主要由石英、钾长石、白云母和钠长石组成，其次是赤铁矿，绿泥石和绿帘石。绿泥石溶解后会产生碳酸盐矿物沉淀所需的Fe2+和Mg2+，可显著促进CO2的矿物封存[21, 24, 33-34]；绿帘石在相关研究中较少涉及，但其金属阳离子含量较高，是碳酸盐矿物沉淀所需二价阳离子的潜在供源。根据原生矿物的类型，参考其他学者的研究成果[13, 21, 33, 35]，也得益于TOUGHREACT并行版的运用，在模型中考虑了几乎所有可能在碳封存过程中形成的碳酸盐矿物和黏土矿物，次生碳酸盐矿物包括方解石、白云石、铁白云石、菱镁矿、菱铁矿和片钠铝石，非碳酸盐次生矿物为高岭石、伊利石、钙/钠蒙脱石和黄铁矿(表3)。

在天然情况下，深层地下含水层中水岩相互反应已达到近似平衡状态。为得到反应溶质运移模型中溶液组分接近平衡的初始状态，在进行CO2注入模拟过程之前，利用储层的矿物组分和盐度为2%的咸水(盐度数据参见文献[31]和[32])建立一个水-岩反应模型，进行1 000 a的模拟反应，将其反应1 000 a后所得的水溶液组分作为反应溶质运移模型的初始水溶液组分。平衡后初始咸水组分见表4。

表2 模型中水文地质参数和热力学参数设置[30-31]

注：Slr为残余水饱和度；Sgr为残余气体饱和度；p0为毛管进气压力；λ为指数。

表3 储层的原生矿物、考虑的次生矿物及其反应动力学参数值[36-37]

注：假设平衡即假设该矿物和溶液可以达到局部平衡状态。

表4 水溶液组分初始浓度和pH值

Table 4 Initial concentrations of water compositions and pH used in the simulations

组分平衡质量摩尔浓度/(mol/kg)组分平衡质量摩尔浓度/(mol/kg)Ca2+1.307e-5HCO-36.902e-2Mg2+2.678e-14SO2-41.000e-10Na+3.215e-1AlO2-1.131e-6K+9.461e-4Cl-2.569e-1Fe2+3.100e-7pH9.05SiO2(aq)2.509e-4

2.3 矿物反应动力学参数

在TOUGHREACT中，矿物反应(溶解或沉淀)速率所采用的模型[38-39]为

(1)

Sg为气相CO2饱和度;b溶解相CO2为溶解相CO2质量摩尔浓度。图1 5 000 a后储层中气相CO2饱和度(a)和溶解相CO2质量摩尔浓度(b)的分布Fig. 1 Distribution of CO2 gas saturation (a) and dissolved CO2 concentration(b)in the aquifer after 5 000 a

其中：r为矿物反应速率(mol/(L·s))；k为速率常数(mol/(m2·s))；A为每千克水的矿物比表面积(cm2/g)；K为平衡常数(无量纲)；Q为相应的离子活度积(无量纲)；θ和η为实验测定的常数，为正值。

对于大多数矿物来说，速率常数k通常是3种机制下的综合[33, 35]：

(2)

其中：k25为25℃时的速率常数(mol/(m2·s))；Ea为活化能(kJ/mol)；R为气体常数(8.31 J/(mol·K))；T为绝对温度(K)；α为组分活度(mol/L)；n为幂指数(常数)；式中上、下标nu, H和OH分别代表中性、酸性和碱性机制。

表3列出了各矿物反应动力学参数的取值。其中除方解石为平衡控制外，其余矿物均为反应动力学控制。

3 模拟结果及分析

3.1 CO2饱和度的分布和pH值的变化

本文为方便起见将超临界CO2流体称为气相。在超临界CO2注入后，由于其自身的密度小于咸水密度，因而在浮力作用下向上运移，遇到盖层阻隔后则沿其底部水平运移。图1显示了气相CO2饱和度和溶解相CO2在储层中的分布。由图1a知，在CO2注入5 000 a后，气相CO2羽的最大半径约为2 600 m。在CO2羽内，CO2溶于咸水中，使咸水密度增加，由此引起的对流作用又促进了CO2的进一步溶解(图1b)[13, 18-19, 23]，使CO2羽本身及其附近区域形成酸性区。图2显示了5 000 a后储层中pH值的分布。由图2可见，酸性区内的pH<5.5，酸性区外pH>10.0。

图2 5 000 a后储层中的pH值分布Fig. 2 Distribution of pH in the aquifer after 5 000 a

3.2 矿物的变化和阳离子反应路径

储层pH值的变化破坏了原有的平衡体系，导致原生矿物溶解和次生矿物沉淀，但由于溶液pH值的差异，酸性区内外矿物的变化差别明显，具体分析如下。

3.2.1 酸性区内矿物的变化和阳离子反应路径

在酸性区内，绿帘石、绿泥石、白云母、钠长石和赤铁矿发生溶解。其中白云母体积分数减少最多(图3a)，绿帘石(图3b)、绿泥石(图3c)和钠长石(图3d)溶解也较明显，赤铁矿溶解很少。绿帘石溶解产生的Fe2+和Ca2+，分别主要被菱铁矿和方解石沉淀消耗(式(3))。绿泥石溶解产生Fe2+和Mg2+，Fe2+也主要被菱铁矿消耗(式(4))，伊利石大量沉淀消耗了Mg2+(式(4)，图3e)，是酸性区内沉淀最多的矿物。钠长石溶解产生的Na+主要被片钠铝石沉淀消耗(式(5))，少部分消耗于钠蒙脱石的沉淀。另外两种原生矿物钾长石和石英在酸性区内发生沉淀，其中钾长石的K+来自白云母(式(6))。主要反应方程式如下：

Ca2FeAl2(SiO4)3(OH)(绿帘石)+3SiO2(石英)+2K++3CO2+0.5H2O→2KAlSi3O8(钾长石)+2CaCO3(方解石)+2H++0.25O2+FeCO3(菱铁矿) ；

(3)

0.25Fe2.5Mg2.5Al2Si3O10(OH)8(绿泥石)+0.6KAl2(AlSi3O10)(OH)2(白云母)+0.95SiO2(石英)+0.75H++0.625CO2→K0.6Mg0.25Al1.8(Al0.5Si3.5O10)(OH)2(伊利石)+0.625FeCO3(菱铁矿)+0.375Mg2++0.975H2O ；

(4)

NaAlSi3O8(钠长石)+CO2+H2O→NaAlCO3(OH)2(片钠铝石)+3SiO2(石英)；

(5)

KAl3Si3O10(OH)2(白云母)+ 6SiO2(石英)+2K+→ 3KAlSi3O8(钾长石)+2H+。

(6)

体积分数变化值为负数表示溶解。图3 5 000 a后矿物的体积分数(a--e)的变化和孔隙度Φ的分布Fig. 3 Changes of minerals in volume fraction (a-e)and distribution of porosity Φ after 5 000 a

3.2.2 酸性区外矿物的变化和阳离子反应路径

在酸性区外，绿帘石、石英和白云母溶解较明显，但由于pH值较高，反应强度小于酸性区内，绿泥石则几乎没有溶解,钠长石、赤铁矿和钾长石发生沉淀。次生矿物中，只有方解石沉淀较明显，而在酸性区内沉淀最多的碳酸盐矿物片钠铝石则没有沉淀，这是由于片钠铝石在高CO2分压下稳定[40]，只有在注入井附近能满足这一条件。根据矿物的变化推知，方解石的Ca2+和赤铁矿的Fe3+均来自绿帘石(式(7))，钾长石的K+来自白云母(式(6))，而由于没有矿物溶解提供Na+，所以钠长石的沉淀消耗了溶液中原有的Na+，反应方程式为

Ca2FeAl2(SiO4)3(OH)(绿帘石)+3SiO2(石英)+2Na++2CO2+0.5H2O→2NaAlSi3O8(钠长石)+2CaCO3(方解石)+0.5Fe2O3(赤铁矿)+2H+。

(7)

在模拟过程中，含镁碳酸盐矿物铁白云石在模拟区域内沉淀极少，菱镁矿和白云石则没有沉淀，原因是伊利石作为酸性区内沉淀最多的矿物，消耗了大量Mg2+，导致含镁碳酸盐矿物沉淀极少。高岭石和钙蒙脱石只在酸性区内有很少量沉淀，黄铁矿则没有沉淀。

3.3 孔隙度和渗透率的变化

图3f显示5 000 a时储层的孔隙度在模拟区域内均降低，在酸性区内降低幅度最大。这是由于沉淀的矿物主要是密度较小的黏土矿物伊利石、碳酸盐矿物(片钠铝石、菱铁矿和方解石)及石英等,而溶解的矿物中，重矿物绿帘石的密度较大，结果导致储层孔隙度降低，最大降低幅度达3.2%。渗透率的变化类似于孔隙度的变化，最大降低幅度达9.3%；kv减小将减弱对流作用[16]，不利于CO2的溶解和矿物封存。

3.4 矿物封存

碳酸盐矿物沉淀能永久固定CO2，图4a、b显示5 000 a时片钠铝石和菱铁矿在酸性区内沉淀均较明显，在酸性区外没有沉淀；图4c、d显示了方解石在1 000和5 000 a时体积分数的变化情况，1 000 a时方解石在酸性区内沉淀很少，只在酸性区外沉淀较明显，但5 000 a时在酸性区内也大量沉淀，并且沉淀量高于酸性区外。这是由于矿物溶解缓冲了pH值，同时绿帘石不断溶解提供丰富的Ca2+，利于方解石沉淀[36, 41]。 图5显示了储层矿物封存量的空间分布及其随时间的变化，表现出逐渐增加的趋势，根据碳酸盐矿物的沉淀情况判断，盐城组下段砂岩层具有较高的矿物封存潜力。储层中所含的少量绿帘石在酸性区内外溶解均比较明显，产生的二价阳离子(Fe2+和Ca2+)主要被菱铁矿和方解石沉淀消耗(式(3)和(7))，说明绿帘石促进了CO2的矿物封存。

a、b、d. t=5 000 a；c. t=1 000 a。图4 片钠铝石(a)、菱铁矿(b)和方解石(c-d)体积分数变化Fig. 4 Changes of dawsonite (a), siderite (b) and calcite (c-d) in volume fraction

a.t=1 000 a；b. t=5 000 a。SMCO2.储层每立方米介质的矿物封存量。图5 不同时刻SMCO2的分布Fig. 5 Distribution of SMCO2 at different time

4 敏感性分析

4.1 网格剖分精度对评估矿物封存量的影响

数值模拟的结果中可能包含离散误差、舍入误差和迭代残差[6]，本次研究拟从网格剖分精度的角度来研究这些误差对评估CO2矿物封存量的影响。为此，在模型A的基础上，建立了如表1所列的粗网格模型B和精细网格模型C。由于模型C计算量太大，因此3个模型都只模拟到1 000 a。

图6显示了模型B和C在1 000 a时矿物封存量的空间分布(相同时刻采用同一图例，下同)。结合图5a可发现，剖分精度越低，模拟所得的矿物封存量越大。1 000 a时，模型A、B和C矿物封存比例分别为4.5%、6.2%和2.7%。剖分最粗的模型B的矿物封存量较大的区域在水平方向上延伸约为1 900 m，但模型A和C均超过了2 000 m。这说明网格剖分越精细，越能精确刻画CO2的运移，这与Juanes等[8]的研究结果类似。值得注意的是，网格剖分精度的不同对离子的反应路径并未产生影响。

4.2 各向异性对评估矿物封存量的影响

为研究储层各向异性对评估CO2矿物封存量的影响，在模型A的基础上，另建立了kv/kh=0.10(模型D)和kv/kh=0.01(模型E)2个模型，每个模型的kh保持一致。图7显示了模型D和E在5 000 a时CO2饱和度的分布。结合图1a可知：kv/kh越小，越多的气相CO2运移到储层的中部和底部；但在储层上部，则是kv/kh值越大CO2运移越远，这与Zhang等[21]的模拟结果类似。在3 000 a之前，kv/kh越小，溶解封存量越高(图8)，说明水平运移对CO2溶解的促进是主要的；但3 000 a时，模型A的溶解封存量超过了模型D，说明系统此时的流体动力学不稳定性已足够大，对流作用强度增大，减小kv对对流作用强度的抑制和减小对流发生的区域的作用开始明显体现。3个模型的溶解封存量在达到最大值后均缓慢减少，这与Kihm等[24]的模拟结果类似；这是由矿物封存量逐渐增加导致的。但3个模型的溶解封存量开始下降的时间不同，具中等kv值的模型D开始最早，而3 000 a后模型A下降最快，这可用对流作用对矿物封存的影响来解释。

图9显示了模型D和E不同时刻的矿物封存量的空间分布，结合图5，可看出对流作用对矿物封存的显著影响。模型A在1 000 a以后矿物封存量较大的区域(位于酸性区及其附近区域)无论在垂向上还是水平方向，都明显大于模型E，这种差异在5 000 a时最明显。图5和图9c中矿物封存量的指状分布也说明对流作用影响的存在，因为对流作用发生时，高密度和低密度咸水的过渡带呈指状分布[18-19]。矿物封存量分布的差异，是由于CO2经过长期的运移和溶解，使系统产生足够大的不稳定性，对流作用已很明显，降低kv抑制了对流作用，进而抑制了CO2和矿物的溶解，不利于矿物封存。此外，5 000 a时，在注入井周围数百米的区域内，模型A和D每立方米介质的矿物封存量要比模型E多 0.5 kg左右(模型D尤其明显)。研究表明，低pH值不利于碳酸盐矿物沉淀[38, 41]。在模型E中，降低kv使更多CO2水平运移开后，CO2的大量溶解使溶液pH降低，虽然pH降低促进了矿物溶解，提供更丰富的阳离子，但低pH环境本身并不利于碳酸盐矿物沉淀。由图8可见，1 000 a后，模型D的矿物封存量超过了模型E，模型A的矿物封存量在3 000 a后快速增加，到5 000 a时超过了模型E(但仍小于模型D)。5 000 a时，模型A、D和E的矿物封存比例分别为34.0%，35.8%和33.4%。

上述结果说明，各向异性(本研究为kh固定，kv降低不同数量级)影响评估矿物封存量。在短期内(本文为1 000 a左右)，对流作用由于系统的流体力学不稳定性还未达到足够大而并不明显，降低kv促进气相CO2的水平运移,从而促进了溶解和矿物封存。1 000 a后，降低kv对对流作用的的抑制开始明显体现。因此在降低kv对水平运移的促进和对对流作用的抑制双重影响下，表现为中等kv/kh值(模型D)相比较大kv/kh值(模型A)和较小kv/kh值(模型E)的各向异性储层的模拟矿物封存量更高，而且在长期封存中，模型A的模拟矿物封存量超过了模型E。对流作用对矿物封存有明显的促进作用。

图6 1 000 a后模型B(a)和模型C(b)的SMCO2分布Fig. 6 Distribution of SMCO2 in model B (a) and C (b) after 1 000 a

图7 5 000 a后模型D(a)和模型E(b)气相CO2饱和度分布Fig. 7 Distribution of CO2 gas saturation in model D (a) and E (b) after 5 000 a

图8 CO2的封存机制及其随时间的变化Fig. 8 Time evolution of the injected CO2 in different trapping mechanisms

a、b. t=1 000 a；c、d. t=5 000 a。图9 模型D(a, c)和模型E(b, d)的SMCO2不同时刻的分布Fig. 9 Distribution of SMCO2 in model D (a, c) and C (b, d) at different time

5 结论

本研究基于TOUGHREACT-MP对苏北盆地盐城组下段砂岩层CO2封存进行了二维CO2-水-岩反应数值模拟，得出了以下结论：

1) 在超临界CO2注入后，储层中的钠长石、绿帘石和绿泥石溶解，促使碳酸盐矿物片钠铝石、方解石和菱铁矿沉淀，促进了CO2的矿物封存。5 000 a后矿物封存比例占34.0%，说明苏北盆地盐城组下段砂岩层具有较高的矿物封存潜力。由于溶解的绿帘石密度较大，而沉淀的矿物密度较小，导致储层孔隙度和渗透率降低。

2) 对不同网格剖分精度的研究结果表明，剖分精度的不同对离子的反应路径没有影响，但粗网格模型计算得到的矿物封存比例偏高，数值误差较大。因此虽然提高剖分精度会增加计算量，但为使结果尽可能精确，模型应达到一定的剖分精度。

3) 储层的各向异性影响对CO2矿物封存量的评估。短期内，降低kv促进了CO2的溶解和矿物封存；但随着时间推移，降低kv对对流作用的抑制会随着流体力学不稳定性的增加而逐渐体现，进而使计算得到的矿物封存量减少。对流作用对矿物封存有明显的促进作用。

[1] Metz B, Davidson A, de Coninck H, et al. IPCC Special Report on Carbon Dioxide Capture and Storage[R]. Geneva: IPCC, 2005.

[2] Dahowski R, Li X, Davidson C, et al. A Preliminary Cost Curve Assessment of Carbon Dioxide Capture and Storage Potential in China[J]. Energy Procedia, 2009, 1(1): 2849-2856.

[3] Krupka K M, Cantrell K J, Mcgrail B P. Ther-modynamic Data for Geochemical Modeling of Carbonate Reactions Associated with CO2Sequestration-Literature Review[R]. Richland: Pacific Northwest National Laboratory, 2010.

[4] Schnaar G, Digiulio D C. Computational Modeling of the Geologic Sequestration of Carbon Dioxide[J]. Vadose Zone Journal, 2009, 8(2): 389-403.

[5] 莫绍星, 李瀛, 龙星皎, 等. 咸水层CO2矿物封存数值模拟研究进展[J]. 地质科技情报, 2013, 32(6): 150-158. Mo Shaoxing, Li Ying, Long Xingjiao, et al. Development of Numerical Simulation for CO2Mineral Sequestration[J]. Geological Science and Technology Information, 2013, 32(6): 150-158.

[6] 薛禹群, 谢春红. 地下水数值模拟[M]. 北京: 科学出版社, 2007. Xue Yuqun, Xie Chunhong. Groundwater Numerical Simulation[M]. Beijing: Science Press, 2007.

[7] Beni A N, Clauser C, Erlström M. System Analysis of Underground CO2Storage by Numerical Modeling for a Case Study in Malmö[J]. American Journal of Science, 2011, 311(4): 335-368.

[8] Juanes R, Spiteri E, Orr F, et al. Impact of Relative Permeability Hysteresis on Geological CO2Storage[J/OL]. Water Resources Research, 2006, 42: W12418[2013-10-23].doi:10.1029/2005WR004806.

[9] Basbug B, Gumrah F. Parametric Study of Carbon Dioxide Sequestration in Deep Saline Aquifers[J]. Energy Sources：Part A: Recovery, Utilization, and Environmental Effects, 2009, 31(3): 255-272.

[10] Green C P, Ennis-King J. Spatial Grid Correction for Short-Term Numerical Simulation Results of Carbon Dioxide Dissolution in Saline Aquifers[J]. Computational Geosciences, 2012, 16(4): 1153-1161.

[11] Ennis-King J, Gibson-Poole C M, Lang S C, et al. Long-Term Numerical Simulation of Geological Storage of CO2in the Petrel Sub-Basin, North West Australia[C]//Proceedings of the Sixth International Conference on Greenhouse Gas Control Technologies. Oxford: Pergamon Press, 2003.

[12] Ennis-King J, Lincoln P. Engineering Aspects of Geological Sequestration of Carbon Dioxide[C]//Proceedings of SPE Asia Pacific Oil and Gas Conference and Exhibition. Melbourne：Society of Petroleum Engineers, 2002.

[13] Audigane P, Gaus I, Czernichowski-Lauriol I, et al. Two-Dimensional Reactive Transport Modeling of CO2Injection in a Saline Aquifer at the Sleipner Site, North Sea[J]. American Journal of Science, 2007, 307(7): 974-1008.

[14] Jahangiri H R, Zhang D X. Effect of Spatial He-terogeneity on Plume Distribution and Dilution During CO2Sequestration[J]. International Journal of Greenhouse Gas Control, 2011, 5(2): 281-293.

[15] Anbar S, Akin S. Development of a Linear Predictive Model for Carbon Dioxide Sequestration in Deep Saline Carbonate Aquifers[J]. Computers & Geosciences, 2011, 37(11): 1802-1815.

[16] Lindeberg E, Bergmo P. The Long-Term Fate of CO2Injected into an Aquifer[J]. Greenhouse Gas Control Technologies, 2003, 1: 489-494.

[17] Han W S, Lee S Y, Lu C, et al. Effects of Permeability on CO2Trapping Mechanisms and Buoyancy-Driven CO2Migration in Saline Formations[J/OL]. Water Resources Research, 2010, 46(7): W07510[2013-10-23]. dio: 10.1029/2009wr007850.

[18] Zhang W, Li Y L, Omambia A N. Reactive Transport Modeling of Effects of Convective Mixing on Long-Term CO2Geological Storage in Deep Saline Formations[J]. International Journal of Greenhouse Gas Control, 2011, 5(2): 241-256.

[19] Pruess K, Zhang K N. Numerical Modeling Studies of the Dissolution-Diffusion-Convection Process During CO2Storage in Saline Aquifers[R]. Berkely: Lawrence Berkeley National Laboratory, 2008.

[20] Ennis-King J, Preston I, Paterson L. Onset of Con-vection in Anisotropic Porous Media Subject to a Rapid Change in Boundary Conditions[J/OL]. Physics of Fluids, 2005, 17:084107[2013-10-23]. http://dx.doi.org/10.1063/1.2033911.

[21] Zhang W, Li Y L, Xu T F, et al. Long-Term Va-riations of CO2Trapped in Different Mechanisms in Deep Saline Formations: A Case Study of the Songliao Basin, China[J]. International Journal of Greenhouse Gas Control, 2009, 3(2): 161-180.

[22] Ukaegbu C, Gundogan O, Mackay E, et al. Si-mulation of CO2Storage in a Heterogeneous Aquifer[J]. Proceedings of the Institution of Mechanical Engineers: Part A: Journal of Power and Energy, 2009, 223(3): 249-267.

[23] Xu X F, Chen S Y, Zhang D X. Convective Stability Analysis of the Long-Term Storage of Carbon Dioxide in Deep Saline Aquifers[J]. Advances in Water Resources, 2006, 29(3): 397-407.

[24] Kihm J H, Kim J M, Wang S, et al. Hydro-geochemical Numerical Simulation of Impacts of Mineralogical Compositions and Convective Fluid Flow on Trapping Mechanisms and Efficiency of Carbon Dioxide Injected into Deep Saline Sandstone Aquifers[J/OL]. Journal of Geophysical Research, 2012, 117: B06204[2013-10-23].doi:10.1029/2011JB008906.

[25] Xu T F, Sonnenthal E, Spycher N, et al. TOU-GHREACT User’s Guide: A Simulation Program for Non-Isothermal Multiphase Reactive Geochemical Transport in Variably Saturated Geologic Media, V1. 2.1[R]. Berkely: Lawrence Berkeley National Laboratory, 2008.

[26] 许天福, 金光荣, 岳高凡, 等. 地下多组分反应溶质运移数值模拟:地质资源和环境研究的新方法[J]. 吉林大学学报: 地球科学版, 2012, 42(5): 1410-1425. Xu Tianfu, Jin Guangrong, Yue Gaofan, et al. Surface Reactive Transport Modeling: A New Research Approach for Geo-Resources and Environments[J]. Journal of Jilin University: Earth Science Edition, 2012, 42(5): 1410-1425.

[27] Zhang K N, Wu Y S, Pruess K. User’s Guide for TOUGH2-MP-A Massively Parallel Version of the TOUGH2 Code[R]. Berkely: Lawrence Berkeley National Laboratory, 2008.

[28] 祝厚勤, 朱煜, 郑开富. 苏北盆地盐城组天然气藏成藏条件及控制因素探讨[J]. 海洋地质动态, 2003,19(9): 22-26. Zhu Houqin, Zhu Yu, Zheng Kaifu. Pooling Conditions and Controlling Factors of Gas Pools in Yancheng Formation of Subei Basin[J]. Marine Geology Frontiers, 2003, 19(9): 22-26.

[29] 翟光明, 王慎言, 邱中建, 等. 中国石油地质志: 卷八[M]. 北京: 石油工业出版社, 1992. Zhai Guangming, Wang Shenyan, Qiu Zhongjian, et al. Petroleum Geology of China: Vol 8[M]. Beijing: Petroleum Industry Press, 1992.

[30] Pruess K. ECO2N: A TOUGH2 Fluid Property Module for Mixtures of Water, NaCl, and CO2[R]. Berkely: Lawrence Berkeley National Laboratory, 2005.

[31] 郑菲, 施小清, 吴吉春, 等. 苏北盆地盐城组咸水层CO2地质封存泄漏风险的全局敏感性分析[J]. 高校地质学报, 2012, 18(2): 232-238. Zheng Fei, Shi Xiaoqing, Wu Jichun. et al. Global Sensitivity Analysis of Leakage Risk for CO2Geological Sequestration in the Saline Aquifer of Yancheng Formation in Subei Basin[J]. Geological Journal of China Universities, 2012, 18(2): 232-238.

[32] 江苏省地质局石油普查队苏北专题队. 苏北中新生界岩矿特征及其在地层对比上的意见[R]. 南京: 江苏省地质资料馆, 1961. Subei Oil Survey Team in Geological Bureau of Jiangsu Province. The Petrological and Mineralogical Characteristics of the Mesozoic and Cenozotic Stratum and Views on the Stratigraphic Correlation in North Jiangsu Province[R]. Nanjing: Geological archive of Jiangsu Province, 1961.

[33] Gundogan O, Mackay E, Todd A. Comparison of Numerical Codes for Geochemical Modelling of CO2Storage in Target Sandstone Reservoirs[J]. Chemical Engineering Research and Design, 2011, 89(9): 1805-1816.

[34] Xu T F, Apps J A, Pruess K, et al. Numerical Modeling of Injection and Mineral Trapping of CO2with H2S and SO2in a Sandstone Formation[J]. Chemical Geology, 2007, 242(3): 319-346.

[35] Xu T F, Apps J A, Pruess K. Mineral Sequestration of Carbon Dioxide in a Sandstone-Shale System[J]. Chemical Geology, 2005, 217(3): 295-318.

[36] Lasaga A C, Soler J M, Ganor J, et al. Chemical Weathering Rate Laws and Global Geochemical Cycles[J]. Geochimica et Cosmochimica Acta, 1994, 58(10): 2361-2386.

[37] Steefel C I, Lasaga A C. A Coupled Model for Tran-sport of Multiple Chemical Species and Kinetic Precipitation/Dissolution Reactions with Application to Reactive Flow in Single Phase Hydrothermal Systems[J]. American Journal of Science, 1994, 294(5): 529-592.

[38] Palandri J L, Kharaka Y K. A Compilation of Rate Parameters of Water-Mineral Interaction Kinetics for Application to Geochemical Modeling[R]. Denver: US Geological Survey, 2004.

[39] Xu T F, Sonnenthal E, Spycher N, et al. TOU-GHREACT：A Simulation Program for Non-Isothermal Multiphase Reactive Geochemical Transport in Variably Saturated Geologic Media: Applications to Geothermal Injectivity and CO2Geological Sequestration[J]. Computers & Geosciences, 2006, 32(2): 145-165.

[40] Worden R H. Dawsonite Cement in the Triassic Lam Formation, Shabwa Basin, Yemen: A Natural Analogue for a Potential Mineral Product of Subsurface CO2Storage for Greenhouse Gas Reduction[J]. Marine and Petroleum Geology, 2006, 23(1): 61-77.

[41] Druckenmiller M L, Maroto-Valer M M. Carbon Se-questration Using Brine of Adjusted pH to Form Mineral Carbonates[J]. Fuel Processing Technology, 2005, 86(14): 1599-1614.

Numerical Modeling of CO2Sequestration in the Saline Aquifer of Yancheng Formation in Subei Basin Using TOUGHREACT-MP

Mo Shaoxing1,Long Xingjiao1,Li Ying1,Zheng Fei1,Shi Xiaoqing1,Zhang Keni2,Zhao Liang1

1.KeyLaboratoryofSurficialGeochemistry,MinistryofEducation/SchoolofEarthSciencesandEngineering,NanjingUniversity,Nanjing210023,China2.CollegeofWaterSciences,BeijingNormalUniversity,Beijing100875,China

The two-dimensional radial model of CO2-water-rock reaction in carbon sequestration of Yancheng Group under section sandstone layer in Subei basin is built based on TOUGHREACT parallel version. The mineral sequestration potential in reservoir is evaluated. The influence of grid subdivision precision and anisotropy of reservoir on mineral sequestration simulations is analyzed. The simulation results show that dawsonite, calcite and siderite, significantly precipitated in the process of CO2sequestration. The total mass of CO2sequestrated by mineral trapping is as high as 34.0% after 5 000 a. The grid resolution has little impact on the reaction path of minerals sequestration, however, the total amount of CO2mineral sequestration with coarser grid is overestimated comparing to that with finer grid. The discretion of the vertical permeability (kv) would lead to a more rapid horizontal migration of injected CO2to the middle and bottom, which results in an incretion in solubility trapping and mineral trapping for a short period of time. However, for a long term, with the increasing instability of the system due to CO2dissolution, the vertical convective mixing becomes significant, it turns out to be a significant inhibition of convection by decreasingkv. It shows that the mineral sequestration calculated by middlekvvalue model is higher than the largerkvand smallerkvvalue model after 1 000 a.

mineral sequestration of CO2; grid subdivision precision; anisotropy; Subei basin

10.13278/j.cnki.jjuese.201405208.

2013-12-23

江苏省自然科学基金项目(BK2012313)；国家大学生创新训练计划项目(G1210284016)；江苏省产学研联合创新资金计划项目(BY2010136)

莫绍星(1992--)，男，布依族，主要从事CO2地质封存技术研究，E-mail:njujinchun@163.com

施小清(1979--)，男，副教授，主要从事反应溶质运移模拟研究，E-mail:shixq@nju.edu.cn。

10.13278/j.cnki.jjuese.201405208

P641.69

A

莫绍星, 龙星皎, 李瀛, 等. 基于TOUGHREACT-MP的苏北盆地盐城组咸水层CO2矿物封存数值模拟.吉林大学学报:地球科学版,2014,44(5):1647-1658.

Mo Shaoxing, Long Xingjiao, Li Ying, et al. Numerical Modeling of CO2Sequestration in the Saline Aquifer of Yancheng Formation in Subei Basin Using TOUGHREACT-MP.Journal of Jilin University:Earth Science Edition,2014,44(5):1647-1658.doi:10.13278/j.cnki.jjuese.201405208.