胡婷¹　芮振华^1,2

1.油气资源与工程全国重点实验室·中国石油大学（北京）

2.中国石油大学（北京）碳中和未来技术学院

0　引言

碳达峰、碳中和的时代背景下，二氧化碳捕集、利用与封存技术（CCUS）作为一种大规模的温室气体减排手段，成为助力“双碳”目标实现的重要技术。地下广袤的储集空间提供了巨大的CO₂地质利用与封存潜力。刻画地下复杂地质条件，并表征CO₂注入后的物理化学过程等相关技术水平的提高至关重要^[1]。注入油藏的CO₂一部分与原油发生相间传质作用导致原油黏度降低，流动性增强，以提高原油采收率；另一部分CO₂溶解于地层水，甚至与岩石反应形成碳酸盐矿物以实现安全封存。CO₂提高采收率与封存的本质是CO₂迁移与相态转化的过程。因此，准确理解CO₂注入储层后，CO₂的迁移转化规律及分布相态演化，对CCUS工程的地质封存储量评估以及安全性评价至关重要^[2]。

数值模拟技术是实现矿场碳封存潜力评价的主要方法，目前被广泛应用的商业化碳封存模拟软件是加拿大发布的CMG软件[3]和美国发布的Eclipse软件。此外，还有TOUGH、UNCONG、CO₂-SCREEN等其他模拟软件。然而，由于没有考虑地层化学反应的非均质性、未配套完善的地球化学反应参数数据库而忽略可能形成的次生矿物，或只能实现用户自行定义的简单水—岩反应等原因，现有的成熟商业软件不能构建符合实际地层条件的地球化学系统。实际地层中矿物种类复杂，除碳酸盐矿物外，还有长石等硅酸盐矿物和黏土矿物等铝硅酸盐矿物，而硅酸盐、铝硅酸盐矿物的溶解正是为矿物态CO₂提供阳离子的主要途径^[4]。矿物的溶解、沉淀反应大部分为反应动力学过程而非瞬间平衡过程，如不同酸碱度（pH值）的反应速率常数不同，影响矿物反应速率的反应表面积随着矿物体积、润湿性及离子半径等参数变化而改变，矿物溶解及沉淀反应活化能参数的差异等，都构成矿物的复杂反应动力学过程^[5]。

目前的数值模拟软件中，用户自定义的简单的矿化反应过程对于CO₂—水—油—岩系统复杂的地球化学反应过程的量化表征是不能完全实现的，对于地层的碳封存潜力，尤其是固相碳储潜力的量化研究尚不完备。为此，通过程序开发手段，将闪蒸计算加入开源的反应溶质运移模拟软件中，开发配套地球化学数据库的热—流—化耦合油藏CCUS数值模拟程序，利用改进后的软件建立松辽盆地大情字井油田H59区块的三维地质模型，厘清了CO₂—水—油—岩系统的相互作用机理，明确地球化学反应过程对渗流的影响机制，通过历史拟合注采过程校正模型地层参数，量化表征不同注入阶段及注采后CO₂迁移与相态转化的时空演化过程，以期解决CO₂高效地质利用与封存的关键问题。

1 　基于热—流—化耦合模拟软件的开发

笔者基于热—流—化耦合TOUGHREACT数值模拟软件开展研究^[6]，该软件是在原有的多相流体流动及热传导模拟软件TOUGH^[7]的基础框架上开发而成，软件中的地球化学反应模块与流体流动模块为顺序耦合^[8]，可考虑不同温度、压力、液体饱和度和离子强度条件下的地下热物理化学过程^[9]。笔者在TOUGHREACT软件的基础上，增加了含多种烃类组分的油相，并将多种烃类组分的状态方程计算和含油相的多相渗流计算耦合到该软件的多相流体流动运算模块中^[10]。

1.1 　相态数及主变量变化

对于储层中可能存在的相态组合如表1所示。在该热—流—化耦合数值模拟软件的基础上进行的模块开发，尽可能地沿用了原有的主变量及其所在位置^[11]。若增加NHC个烃组分，加上水组分，那么整个系统的组分个数N_k=N_HC+1。N_HC个烃组分，需用（N_HC-1）个烃组分摩尔分数进行状态刻画。保留温度（T）、压力（p）参数及其位置，根据系统相态存在的情况主变量设置如表1所示。其表达式为：

式中z_i定义为某一相中烃组分i在总烃组分中的摩尔分数；m_HC,i 表示某一相中烃组分i的摩尔质量，mol；m_HC,k表示某一相中第k个烃组分的摩尔质量，mol；N_HC表示总烃组分数，个。

表1    相态组合及主变量设置表

保留原软件中第一个相态位置（即序号为N_PH＝1）为气相指针位置、第二个相态位置（即序号为N_PH＝2）为水相指针位置，设定油相的热物性参数数组序号为N_PH＝3，代替原有NaCl析出形成的固相序号位置。因此，更改序列后的程序将原本N_PH＝3序列中的固相饱和度与孔隙度的变化相联结。各相态满足的关系式的变化如下式所示：

式中V_g、V_w、V_o分别表示地层中气相、液相、油相态的体积，m³；ϕ表示地层的孔隙度；V_m表示地层的体积，m³。

1.2 　状态方程耦合多相渗流计算

针对系统中含有多组分烃的特征，利用PengRobinson（PR）立方方程^[12-13]对多组分在气相或油相中的逸度进行计算，以此确定组分的相平衡状态。该方程利用单组分方程结合二元相互作用系数计算多组分混合物，考虑了多组分间的相互作用，保证预测结果更准确，方程有2种表达形式，如式（3）、（4）所示。二元交互作用系数aij通常由实验室测得的数据建立数据库确定，其他如分子量、标态性质、偏心因子以及理想气体的比热系数等参数取自于已发布的属性数据库^[14]。

1.3 　耦合地球化学过程

在化学求解模块中，所需求解的过程包括：①气体组分CO₂溶解于水相形成HCO₃^–的过程；②水相中溶质迁移过程；③矿物溶解及沉淀过程。因此，在完成流动及热流模块的计算后，分别调用序列矩阵N_PH=1气相中CO₂分压及N_PH＝1气相、N_PH＝2水相的流体流动速率及上一个时间步的地球化学系统的变化参数进行上述过程的计算。

如图1所示，首先将考虑油相的流动及热流计算模块求解的流体流动速率导入溶质运移模块，通过溶质运移计算确定各相态中的组分分布^[11]；之后，将气相中CO₂分压及水相中的各组分浓度分布导入化学计算模块进行单网格的动力学平衡计算及收敛判断；收敛平衡后，将更新的孔隙度（ϕ）、绝对渗透率（K）及CO₂质量的变化量（）再传回流动及热流计算模块，至此完成二次开发软件中的流动及热流计算模块、化学结算模块的顺序耦合过程，实现对改进的流动及热流计算模块与原有的化学计算模块的搭接^[18]。

图1    CCUS-EOR数值模拟程序工作流程图

1.4 　CO₂分布相态的计算方法

注入的CO₂最终变为4种不同的分布形态，即气相、溶解于水相、溶解于油相、矿物固相。由于气相的密度随压力变化大，因此，气相CO₂封存量用质量守恒方法进行计算^[19]，各相态中CO₂分布量在改进的软件中将通过以下公式进行计算：

2 　数值模拟及历史拟合

2.1　模型构建

以中国CCUS工业化应用的示范场地——松辽盆地大情字井油田H59区块为例，通过收集场地地层参数特征及井流监测数据，构建研究区的三维数值模型^[20]。目标层上白垩统嫩江组黑帝庙油层为轻质油，油密度为0.765 kg/m³。前人已开展的细管实验表明，该油层的最小混相压力为22.3 MPa，地层破裂压力为49.0 MPa^[21-22]，因此，CO₂注入后可保证地层安全的条件下与原油发生混相^[23]。H59区块采用正方形反九点井网进行开采，其目标储层深度介于2 230～2 280 m，岩性以长石岩屑砂岩为主，平均孔隙度为12.7%，平均渗透率为3.0 mD^[24]。根据收集的井流数据资料，选取H59区块边缘位置的一个注采单元中的注入井H59-1井及开采井H59-3-3井建立注采模型，复现目标注采井13个月的历史注采过程，并分析注采期内注入碳的迁移、转化及分布相态特征，定量刻画注采结束后第10年、第100年时储层内注入CO₂的分布相态演化。

2.1.1　概念模型和边界条件

目标注入井H59-1井及开采井H59-3-3井间距200 m，建立的单注单采三维场地数值模型不考虑H59-1井注采单元内其他开采井的影响。设置侧向边界位置距离注采井5 000 m及3 000 m，边界网格设置体积无穷大，并设定为第一类边界条件，降低边界效应的影响。大情字井油田H59区块黑帝庙油层单层厚度介于11.2～18.2 m^[25]，模型垂向厚度设置为12 m，井周网格相较于其他区域进行了加密剖分，加密部分网格数量为10×10×4共400个网格，每个网格在x、y、z方向上的距离分别为20 m、20 m、3 m（图2）。概念模型设置在H59-1井储层底部注入，在H59-3-3井储层顶部开采，注入段和采出段储层厚度均为3 m。

图2    大情字井油田H59区块概念模型图

2.1.2　初始物性条件

根据大情字井H59区块的地质资料，取各参数平均值为模型的初始物性条件，如表2所示。

表2    H59区块黑帝庙油层初始物性条件数据表^[26]

初始原油拟组分依据前人研究结果设置，具体参数如表3所示。

表3    H59区块黑帝庙油层初始原油拟组分数据表^[10]

相对渗透率参数^[19]设定如图3所示，依据水—油及气—油相对渗透率参数，利用STONE II模型计算气—水—油三相系统中的流体相对渗透率。STONE II模型中相关参数设置如表4所示。

图3    H59区块黑帝庙油层储层流体相对渗流率曲线图

（资料来源：据本文参考文献[19]）

表4    H59区块黑帝庙油层STONE II模型参数设置表

2.1.3　初始化学条件

由于现场资料有限，模型中原始地层水的化学条件设定参考雷浩^[26]研究中的大情字井油田H59区块的地层水化学组成进行设定，如表5所示。初始矿物组成设定为所收集资料中矿物含量的平均值^[27]，矿物组成及矿物的反应动力学参数如表6所示。

表5    H59区块黑帝庙油层初始地层水离子组成及摩尔浓度数据表^[26]

表6    H59区块黑帝庙油层初始矿物组成及其反应动力学参数数据表

2.2 H59区块注采史拟合与预测

已开发软件的可靠性和准确性已通过严格的对比验证，并进行场地应用^[10]。针对目标区地质模型的参数设置进行历史拟合，定量分析目标区CO₂的迁移与分布相态规律。

2009年5月至2010年5月，对注入井H59-1井井底压力进行了连续监测，井流数据如图4所示，图中显示了监测期内H59-1井注入水及CO₂的注入速率。H59-1井分3个阶段进行注入，分别为注水阶段、注CO₂阶段和气水同注阶段，监测数据显示井底压力平均值为33.2 MPa。运用数值模型对监测的注入井井底压力进行历史拟合，并计算出开采井井底压力为14.5 MPa。根据井流数据计算出向储层中注入水和CO₂的累计量分别为0.654×10⁴ t、0.942×10⁴ t（图5），依据井流数据复现H59-1井的流体注入过程，分析3种典型注入方式下的流体运移、转化规律。

图4    H59-1井生产曲线图

（资料来源：据本文参考文献[22]）

图5    H59-1井的实际注入曲线图

H59区块H59-1井注采单元包括注入井H59-1井和开采井H59-3-3井，开采井H59-3-3井位于注入井H59-1井西南200 m处。数值模拟计算的生产数据如图6所示，3个不同注入阶段的日产油量、日产气量均表现为增大的趋势，在气水同注时日产油量、日产气量达到最大，而日产水量基本不变。这主要是由于注气阶段注入CO₂在注采井间大面积溶解于油相，改善了油相的流动性，有利于后续流体的注入及运移^[28]。此外，注采期3个阶段共采出已注入的CO₂量为0.045×10⁴ t，主要来源于溶解在油相的CO₂，最终保留在储层中的注入的CO₂量为0.897×10⁴ t。在注气阶段转变为气水同注阶段时，模型计算得到的CO₂采出量迅速下降的现象与实测值变化情况一致，这说明软件具有较好的场地应用效果。此外，在整体上计算值较实测值均匀分布，这主要是由于实际条件下生产数据受多种因素的影响，如地层非均质性、井间干扰等^[29-30]。

图6    开采井H59-3-3井的模型计算产量曲线图

（资料来源：实测数据来源于本文参考文献[22]，有修改）

3 　注采期碳迁移和分布相态特征

3.1  CO₂迁移及空间分布特征

对注水阶段、注气阶段及气水同注阶段地层含气饱和度及含水饱和度的空间分布变化情况进行分析发现，在注采过程中，注气阶段及气水同注阶段注入的CO₂只在注入井端小范围内呈现气相，其余位置由于混相而以油相形式存在；含水饱和度的变化显示，气水同注阶段注入水的运移范围明显大于注水阶段，说明注CO₂能进一步增强后续注水的波及效率（图7）。

图7    数值模拟注采期注采井间流体空间分布剖面图

注水前、注水阶段、注气阶段及气水同注阶段含油饱和度的空间分布如图8所示，储层内含油饱和度在注水后显著降低，由于在顶部开采，含油饱和度在不同深度上均呈现斜向上的蝶形分布。在注气阶段及气水同注阶段，由于CO₂与油相发生混相，靠近油相一侧注入前缘的含油饱和度大面积增加。此外，溶解CO₂的油相密度降低而向储层顶部扩散，驱动储层顶部油相向下运移并与CO₂相互作用，不断的传质及扩散过程导致整个储层厚度上的含油饱和度趋于一致^[31]。

图8    数值模拟注采期注采井间含油饱和度剖面图

注水阶段、注气阶段及气水同注阶段油相中CO₂摩尔分数的空间分布如图9所示，注气阶段气相CO₂最远迁移至距注入井160 m处。气水同注阶段CO₂已运移至开采井位置，由于气、水的渗流速度存在差异，图中局部低值显示注入气与注入水的渗流距离出现明显差异现象，局部低值区域则是水段塞分布区。油相中CO₂摩尔分数以注入点为中心，呈现弧状分布，这说明CO₂溶解油相后浮力作用减弱，油相中CO₂组分在垂向上趋于相对均匀分布^[32]。

图9    数值模拟注采期注采井间油相中CO₂的摩尔分数剖面图

3.2　注采期CO₂分布相态特征

确定地层中注入CO₂的分布相态是明确CO₂封存稳定性的关键，也是评估地层碳封存能力的前提。利用更新后的软件对2009年5月至2010年5月间注入CO₂在气相、溶解水相、溶解油相和矿物固相4种分布形态进行了定量表征（图10）。由图中可知，注采期内矿物以溶解为主，CO₂不存在以固相封存的相态。注气阶段CO₂主要溶解于油相，占比70.8%；其次为气相，占比16.3%；溶解水相的CO₂占比12.9%。在气水同注阶段，气相CO₂含量出现明显拐点，该注入方式导致CO₂溶解水相的量明显增加，同时气相中CO₂含量降低。考虑到溶解于水的CO₂能为矿物固相封存提供前提条件，因此，气水同注是提高CO₂提采与安全封存双赢的有效方式^[33]。气水同注阶段整体上占比最多的仍为溶解于油相的CO₂，占比72.2%；其次为溶解于水相的CO₂，占比17.8%；溶解气相的CO₂占比10.0%。

图10    数值模拟注入井H59-1井及开采井H59-3-3井注采期储层内CO₂分布相态图

4 　注采后碳迁移演化特征

4.1 CO₂分布相态演化

5　结论

1）研究区油藏CO₂混相驱条件下，CO₂在注入井端的小范围内呈现气相，在接触到油相前缘后，CO₂会迅速与原油发生混相且受浮力作用减弱，油相中CO₂在垂向上趋于均匀分布，并向开采井端均匀推进。

2）气水同注的方式导致CO₂溶解水相的含量明显增加，为CO₂的稳定固相封存提供了条件，同时还保证了大部分注入的CO₂仍溶解于油相，这表明气水同注的方式是实现CO₂提采与封存双赢的有效方式。

3）注采结束时，注入的CO₂主要溶解于油相，占比72.2%；其次溶解于水相，占比17.8%；气相占比10.0%，矿物固相可忽略不计；注采结束10年时，CO₂分布相态特征表现为溶解水相CO₂逐渐转变为矿物固相，溶解在油相中的CO₂存在逐渐转变为游离态气相的趋势，模型计算的单井注入量中矿物封存量占总量的0.01%，固相态封存量呈现稳定增长，其他相态含量均有所降低；注采结束100年时，地层中主要的沉淀矿物有方解石、钠蒙脱石及菱铁矿，主要的溶解矿物有钾长石、钙蒙脱石、铁白云石及绿泥石。

4）大情字井油田H59区块固碳矿物主要为方解石及菱铁矿，主要矿化过程为铁白云石及绿泥石溶解产生铁离子、钙离子与碳酸根结合形成固碳矿物方解石和菱铁矿。