地层各向异性对CO2在含水层中迁移转化的影响

李恋， 吴亚红*， 种珊， 温庆志， 林千果

(1.中国石油大学(北京)石油工程教育部重点实验室， 北京 102249； 2.中国矿业大学(北京)煤炭资源与 安全开采国家重点实验室， 北京 100083; 3.北京大学工学院工程科学与新兴技术高精尖创新中心， 北京 100871； 4.华北电力大学苏州研究院二氧化碳捕集实验室， 苏州 215000)

CO2捕集、利用与封存(carbon capture, utilization and storage, CCUS)是一种既能有效利用CO2提高油藏采收率，又能实现CO2地质封存的新技术[1-3]。CCUS技术对全球长期应对气候变化和推进落实好“碳中和”具有重要意义[4]。随着CO2驱油项目在世界各地的广泛开展，其技术的安全性和潜在的泄漏风险越来越引起学者们的广泛关注[5]。由于地层的不连续性，封存的CO2可能沿断层、裂缝、开采井和废弃井等途径泄漏到浅层含水层[6]。CO2一旦泄漏到浅层含水层，会打破浅层含水层原来的水化学平衡条件，对浅层含水层造成一系列的影响[7-8]。若CO2突破含水层向上迁移扩散至大气中，不仅会降低CO2驱油封存的效果，而且会对该区域自然环境和人文生活造成一定的影响[9]。

目前研究CO2泄漏对浅层含水层水质的影响主要通过现场试验、室内CO2注入实验以及数值模拟方式。Do等[10]通过收集地下水化学成分和同位素数据，分析了在三种地下水化学类型中CO2长期泄漏后对地下水质量的影响，尤其是微量金属元素的富集。Delkhahi等[11]通过水样分析，探讨了水井中CO2泄漏对含水层中主、次离子的地球化学变化以及硫酸盐的还原作用与同位素δ34SSO4和δ18OSO4之间的关系。Szabó等[12]采用PHREEQC软件研究了CO2泄漏对浅层饮用含水层中离子浓度的影响并估算了地下水的电导率和岩石的孔隙度。Li等[13]通过建立三维反应流模拟模型，评估了储层温度、注水井与渗漏层的距离对CO2泄漏的影响。Fahrner等[14]利用PHREEQC软件模拟了方解石含量、CO2泄漏速率、含水层深度、初始盐度和CO2泄漏时间对地下水电导率的影响。Petit等[15]通过数值模拟研究了CO2泄漏对淡水含水层的电导率、温度和溶解态CO2浓度的影响。

虽然CO2泄漏对浅层含水层的影响已经开展了较为丰富的研究，然而，当前的现场试验和室内实验研究由于条件的限制，难以掌握CO2泄漏在含水层中迁移扩散的空间分布特征。数值模拟研究中主要探讨CO2泄漏速率、温度、初始盐度、含水层深度和泄漏时间对含水层的pH、电导率和微量金属元素迁移的影响。一方面模型中尚未考虑地层各向异性的影响，导致CO2泄漏后突破含水层的时间难以掌握，另一方面缺少对CO2在含水层中迁移转化过程的分析以及地下水化学成分变化特征的研究。

因此，为了解决上述问题，现以延长油田杏子川化子坪试验区白垩系洛河组砂岩含水层为例，采用数值模拟方法建立CO2在浅层含水层中的反应性迁移转化模型，研究地层各向异性对CO2在浅层含水层中迁移转化特征的影响，以期为CCUS项目CO2安全监测提供一定的参考依据。

1 模型设置

1.1 概念模型和初始边界条件

采用非等温-多相流体反应运移模拟软件TOUGHREACT-OMP/ECO2N[16]进行模拟，该软件是在非等温多相流体模拟软件TOUGH2的基础上引入溶质运移和化学反应而开发的一款地下多组分反应运移模拟程序，主要用于模拟多相流体在多孔介质或裂隙中迁移运动及地球化学反应[17]。

以延长油田杏子川化子坪白垩系洛河组砂岩含水层为例，地层概化为均质各向异性、侧向无限延伸的地层。水平方向渗透率为0.5×10-13m2，垂向渗透率为0.5×10-14m2，孔隙度为0.2，岩石热传导系数为2.51 W/(m·℃)，孔隙压缩系数为4.5×10-10Pa-1，扩散系数为1.0×10-9m2/s，岩石比热为920 J/(kg·℃)[18]。含水层的埋深为10 m，垂向厚度为60 m，平均剖分为6个网格，水平方向的长度为500 m，平均剖分为50个网格，侧向边界为一类边界条件。该区地面年平均地温为10.8 ℃，地温梯度约3.1 ℃/100 m，静水压强梯度为10 bar/100 m(1 bar=100 kPa)，底板温度为13 ℃，压力为8 bar。CO2泄漏点源在X=250 m处，CO2泄漏速率为0.000 1 kg/s。CO2泄漏时间为100年。概念模型如图1所示。

图1 概念模型示意图Fig.1 Schematic diagram of conceptual model

模型中相对渗透率和毛管压力模型如表1所示，具体参数参考Xu等[19]。

1.2 地层矿物组层

地层的初始矿物组成和可能生成的次生矿物及其体积分数如表2所示。

1.3 浅层含水层的化学组分

初始水化学组分为杏子川化子坪试验区地层水离子浓度的检测数据，由于脱离了原位环境，因此，需要对含水层的化学组分与表2中的地层矿物进行化学平衡计算，浅层含水层化学组分浓度如表3所示。

根据研究区域的地质资料，该地区局部为粉砂岩和砂砾岩，地层渗透率存在一定差异，地层各向异性敏感性分析模拟方案如表4所示。

表1 相对渗透率和毛管压力模型Table 1 Models of relative permeability and capillary pressure

表2 地层矿物组成Table 2 The mineral composition of the formation

表3 浅层含水层化学组分的浓度Table 3 Concentration of chemical components in shallow aquifers

表4 地层各向异性敏感性分析模拟方案Table 4 Simulation program for analyzing sensitivity of parameters

2 模拟结果与讨论

2.1 CO2浓度的分布特征

CO2泄漏到浅层含水层后，气态CO2会在浮力的作用下向上迁移并逐渐向周围扩散，含水层中CO2气体饱和度的分布范围逐渐扩大，且在泄漏点处CO2气体饱和度最高(图2)。到第20年时，气态CO2扩散到浅层含水层的顶部[图2(c)]。另一部分CO2会溶解在含水层中，随着CO2逐渐向周围迁移扩散，含水层中溶解态CO2浓度范围逐渐扩大(图3)。由于压力的不同，含水层底部溶解的CO2浓度相对较大，而顶部溶解的CO2的浓度相对较小。在CO2泄漏点处，溶解态CO2的浓度最高。

2.2 CO2在浅层含水层中的迁移转化特征

图2 基础模拟方案CO2气体饱和度的分布Fig.2 Distribution of CO2 gas saturation in basic simulation scheme

图3 基础模拟方案溶解态CO2的浓度分布Fig.3 Distribution of dissolved CO2 in basic simulation scheme

X为水平长度；Z为垂向深度图4 基础模拟方案浅层含水层pH的分布Fig.4 Distribution of pH values in shallow aquifers in basic simulation scheme

由于浅层含水层的酸性增强，pH显著降低，导致初始矿物发生溶解。由图5(a)可知，在0～5年，由于含水层的酸性明显增强，H+浓度显著增加，方解石的溶解量显著增多。从第5年开始，由于方解石的溶解作用达到稳定状态，方解石的体积分数基本保持稳定。奥长石的体积分数变化量随时间线性增加，方解石和奥长石的溶解显著增加了含水层中Ca2+的浓度。钾长石和绿泥石的体积分数变化量相对较少，而赤铁矿的体积分数变化量接近于0。伊利石在0～10年时溶解，由于钾长石和绿泥石的溶解为伊利石的沉淀提供了大量的K+、Mg2+和Al3+，从第10年开始，伊利石开始沉淀，到第100年时，伊利石的体积分数增加了6.27×10-7。K+浓度增加主要来源于钾长石和伊利石的溶解，由于伊利石从第10年开始沉淀消耗部分K+，因此，从第10年开始K+浓度增加幅度开始变缓[图5(d)]。溶解反应如下：

(1)

(2)

2Mg2.5Fe2.5Al2Si3O10(OH)8(绿泥石)+

6H4SiO4

(3)

(4)

(5)

(6)

奥长石和绿泥石的大量溶解使含水层中Na+和Mg2+的浓度显著增加，进而生成钠蒙脱石沉淀。当含水层中溶解态CO2浓度增加时，由于含水层的酸性明显增强，钾长石会逐渐转化为高岭石[式(8)]。石英、高岭石和钠蒙脱石的沉淀量逐渐增多。伊利石和绿泥石的溶解显著增加了含水层中Mg2+的浓度，Mg2+的浓度增加了1～2个数量级，其浓度变化呈现先增大后减小的趋势。由于钠蒙脱石沉淀消耗大量的Mg2+，因此，Mg2+浓度在10年时出现峰值后显著降低，但仍高于初始值[图5(e)]，说明伊利石和绿泥石溶解后产生的Mg2+含量仍然大于钠蒙脱石沉淀消耗的Mg2+。石英和高岭石的沉淀量相对较少[图5(b)]。沉淀反应为

0.29Na++0.26Mg2++3.97H4SiO4+

(7)

(8)

CO2-水-岩石相互作用会导致次生矿物的生成，方解石和奥长石的溶解为次生矿物钙蒙脱石的生成提供了大量的Ca2+，钙蒙脱石和菱铁矿的体积分数随时间逐渐增大[图5(b)]，沉淀量逐渐增多，其中菱铁矿的体积分数变化量最大。绿泥石的溶解为菱铁矿的沉淀提供了大量的Fe2+，结合图5(a)和图5(c)可知，绿泥石体积分数的变化量比菱铁矿体积分数的变化量高2～3个数量级，到第10年时，菱铁矿体积分数变化量明显增大，因此，Fe2+的浓度在初始阶段显著增加，到第10年时，由于菱铁矿沉淀量显著增加，消耗部分Fe2+，Fe2+浓度增加趋势开始变缓[图5(e)]，次生矿物反应为

(9)

0.145Ca2++0.26Mg2++3.97H4SiO4+

(10)

2.3 地层各向异性敏感性分析

kv/kh(kv保持不变)增大时，CO2在水平方向的迁移扩散明显受到限制，水平方向的扩散范围逐渐缩小(图6)，CO2将不得不沿高渗透层迁移扩散，这增加了气态CO2与地层水在垂向上的接触时间和接触面积，促进了CO2溶解。在CO2泄漏点处，溶解态CO2的浓度最高(图7)，从第10年开始，泄漏点处溶解态CO2浓度逐渐达到饱和状态，溶解态CO2浓度基本保持不变。垂向上溶解的CO2浓度显著增加，导致pH显著降低[(图8)]。当kv/kh从0.01增大到0.5时，溶解态CO2浓度增加了0.35 mol/L[图9(a)]，浅层含水层的pH降低了0.2[图9(b)]，CO2在水平方向上的迁移扩散距离缩小了82 m[图10(a)]，气态CO2突破浅层含水层的时间缩短了约30年[图10(b)]。说明在驱油封存过程中，若CO2泄漏到浅层地下水且地层水平方向的渗透率较小时，CO2会在短时间内突破浅层含水层，从而迁移扩散到含水层的上部区域，这不仅会导致潜在的污染风险，而且会进一步扩散至大气中，导致温室效应加剧。

负值表示矿物溶解；正值表示矿物沉淀图5 矿物体积分数变化量和离子浓度随时间的变化Fig.5 Changes in mineral volume fraction and ion concentration over time

图6 CO2气体饱和度的分布Fig.6 Distribution of CO2 gas saturation

图7 溶解态CO2的浓度分布Fig.7 Distribution of dissolved CO2

图8 浅层含水层pH的分布Fig.8 Distribution of pH values in shallow aquifer

图9 泄漏点处溶解态CO2浓度和pH的变化Fig.9 Changes of dissolved CO2 concentration and pH at the leakage point

图10 CO2迁移扩散距离随时间的变化Fig.10 Changes in the distance of CO2 migration over time

3 结论

以延长油田杏子川化子坪试验区白垩系洛河组砂岩含水层为例，采用数值模拟软件TOUGHREACT研究了地层各向异性对CO2在浅层含水层中迁移转化的影响，得到以下结论。

(1)CO2泄漏到浅层含水层后，浅层含水层的pH显著降低，气态CO2和溶解态CO2浓度逐渐升高，其中在泄漏点处气态CO2和溶解态CO2浓度最高。

(2)溶解矿物为方解石、钾长石、绿泥石、奥长石和赤铁矿；沉淀矿物为钠蒙脱石、高岭石和石英，伊利石从第10年开始沉淀；次生矿物为钙蒙脱石和菱铁矿。

(3)地层垂向渗透率与水平渗透率比值kv/kh(kv保持不变)越大，浅层含水层中溶解态CO2浓度越高，含水层的pH越低，CO2在水平方向上的迁移扩散距离缩小了82 m，CO2突破浅层含水层的时间缩短了约30年。