高冉 吕成远 周凯 伦增珉 周冰

1.中国石化石油勘探开发研究院博士后工作站；2.中国石化石油勘探开发研究院；3.中国矿业大学地球科学与测绘工程学院

到2017年，全球大气CO2含量已经达到4.05%。CO2的温室效应导致全球气候变暖、海平面上升、极端气候增加等一系列环境危害，影响着人类的生活［1］。国内外研究表明，CO2捕集利用与埋存技术可以有效地缓解温室效应［2］。当前可行的埋存方式主要为地下埋存、海洋埋存、森林和陆地生态埋存。CO2的地下埋存技术较成熟，埋存潜力大［3-4］。地下埋存主要选择枯竭油气藏(CO2-EOR)、深部的咸水层以及不可开采的煤层(CO2-ECBM)等。CO2-EOR不仅减少温室气体CO2的排放，同时还大幅度提高了原油的采收率。CO2-EOR，即CO2驱的CO2埋存潜力计算对油藏开发方案设计和安全封存意义重大［5-6］。

CO2驱埋存潜力计算方法主要分为静态法和动态法。静态法计算得到的是理论埋存潜力，以美国能源部(USDOE)、碳领导人论坛(CSLF)、Bachu和Brennan等人的研究为代表［7］。基于物质平衡方程，对CO2埋存潜力进行粗略计算。通过静态法得到的理论埋存潜力是最高上限的埋存潜力。

动态法计算得到的是有效埋存潜力。现有的动态法主要有3类：物质平衡动态法［8］、类比法动态法［9］；蒙特卡罗动态法［10］。基物质平衡的动态方法考虑了浮力、重力超覆、流度比和非均质性等因素，但是没有方法能够准确定量计算这些因素对埋存潜力的影响。类比法通过与已开发的相似油藏的对比来获得CO2利用系数，从而计算CO2埋存潜力。但是，该方法不能考虑油藏的生产对埋存潜力的影响。蒙特卡罗方法可以计算出相对准确的埋存潜力，但是蒙特卡洛方法依赖于大量的油田数据，需要耗费大量的计算时间。对于新开采的油藏，由于缺乏足够的数据，不太适用。提出了一种基于组分闪蒸运算的CO2驱动态埋存潜力计算方法。方法基于组分闪蒸运算，考虑了油田的生产实际和CO2驱的埋存机理，可以计算得到溶解CO2、束缚CO2、自由CO2和总的CO2埋存潜力。同时，选取了2类油藏进行了CO2驱动态埋存潜力计算。

1 计算机理

1.1 CO2埋存机理

CO2的地下埋存机理［11］主要有构造地层静态埋存、束缚气埋存、溶解埋存和矿化埋存。不同的埋存机理具有不同的作用时间尺度，如图1所示。

图1 各种埋存机理随时间变化的贡献比例和安全性示意图Fig.1 Schematic contribution ratio and safety of different storage mechanisms over the time

不同的埋存机理产生不同类型的CO2，如图2所示［12］。构造埋存指CO2封存在构造圈中，存在于构造中的可以自由流动的CO2叫做自由CO2。束缚埋存是由气液相界面张力不同导致CO2被长期滞留在岩石孔隙中，存在束缚空间中的CO2叫做束缚CO2。溶解埋存即CO2在运移的过程中溶解在地层水和地层油中，由于溶解作用而埋存的CO2叫做溶解CO2。矿化埋存是地层中的CO2与地层水发生化学反应产生碳酸类矿物沉淀从而长久的滞留在地下的过程。矿化埋存作用时间为100~10 000年。影响矿化埋存的主要因素为地层岩石的矿物成分、流体类型和化学反应类型。由于矿化埋存的时间尺度比较大，基于组分闪蒸运算的CO2驱动态埋存潜力计算方法未对矿化于岩石中的CO2进行计算。

图2 油藏埋存机理示意图Fig.2 Sketch of storage mechanism

1.2 闪蒸计算方法

闪蒸计算包括稳定性测试和相分离2个过程。在相分离过程中，需要求解一个相平衡方程(也叫质量守恒方程)和一个逸度平衡方程。一般采用分步连续迭代方法对这2个方程进行有序求解：首先利用初始的相平衡常数求解相平衡方程；然后把用相平衡方程求解出来的相平衡常数带入到状态方程和逸度方程中进行求解，从而更新出新的相平衡常数。用新的相平衡常数验证初始相平衡常数，从而不断迭代计算，得到正确的解，即确定混合物在指定温度、压力条件下组分在各相中的摩尔百分含量［13］。

油-气两相闪蒸计算方法已经非常成熟，因而可以用商业油藏数值模拟软件如Eclipse，CMG等进行该部分的求解计算［14］，也可以自主编写程序进行闪蒸求解。基于组分闪蒸运算的CO2驱动态埋存潜力计算方法针对这2种求解方法提出了2套CO2驱埋存潜力计算公式。

2 计算方法

基于组分闪蒸计算和CO2埋存机理，耦合CO2在水中的溶解，从而计算自由CO2、溶解CO2、束缚CO2和总的CO2埋存潜力。

2.1 计算溶解CO2质量

地层有油、气和水相，CO2既溶于油相，又溶于水相，因此溶解CO2包括CO2在油相中的溶解和CO2在水相中的溶解2个部分。

2.1.1 CO2在油相中的溶解质量

(1)CO2在油相中的溶解由3步计算得到，由式(1)计算油相的摩尔数为

(2)由组分闪蒸计算得到CO2在油相中的摩尔百分数x；

(3)由式(2)计算CO2在油相中的溶解质量为

式中，no为油相的摩尔数，mol；V为油相总体积，m3；so为含油饱和度；Vmo为油相摩尔体积(闪蒸计算得到)，m3/mol；So为CO2在油相中的溶解质量，kg。

若使用数值模拟软件Eclipse进行组分闪蒸计算，由式(3)计算CO2在油相中的溶解质量为

式中，Vi，soi和xi分别为每个网格的体积、含油饱和度、CO2在油相中的摩尔百分数；Boi为每个网格的原油摩尔密度，mol/m3。Boi，soi和xi分别对应Eclipse软件自带的关键字BOIL，SOIL和XMF1。进行组分模拟计算时，选择输出每一步的这几个关键字的计算结果，根据式(3)编写程序得到溶解在油中的CO2的质量。

2.1.2 CO2在水相中的溶解质量

CO2在水中的溶解由溶解度公式计算得到［15］

若使用数值模拟软件Eclipse进行组分闪蒸计算，计算CO2在水相中的溶解质量为

式中，Sw为CO2在水相中的溶解质量，kg；mw为水的质量，kg；c为CO2在水中的溶解度，g/100 g；swi为每个网格的含水饱和度，对应Eclipse软件中的关键字SWAT；ρi为CO2在压力pi下的密度，kg/m3；ci为每个网格的CO2在水中的溶解度，g/100 g，对应Eclipse软件中的关键字RSW。总的溶解CO2质量S为CO2溶解在水中的质量Sw与CO2溶解在油中的质量So之和。

2.2 计算自由CO2和束缚CO2的质量

自由CO2和束缚CO2的质量分2种情况计算。

(1)当含气饱和度大于束缚气饱和度，即sG＜sg，自由CO2和束缚CO2同时存在，求解自由CO2和束缚CO2的质量为

求解自由气摩尔数和束缚气摩尔数为

式中，sG为束缚气饱和度；sg为含气饱和度；F为自由CO2质量，kg；y为CO2在气相中的摩尔百分含量，由闪蒸计算得到；ng为气体的摩尔数，mol；nG为束缚气摩尔数，mol；T为束缚CO2质量，kg；Vmg为气相摩尔体积，m3/mol，由闪蒸计算得到。

(2)当束缚气饱和度大于等于含气饱和度，即sG≥sg，则只有束缚CO2，没有自由CO2，即F=0。求解公式为

同样，可以用数值模拟软件Eclipse进行组分模拟闪蒸计算，再计算自由CO2和束缚CO2的质量。若sG＜sg，既有自由CO2也有束缚CO2，求解公式为

若sG≥sg，只有束缚CO2，即F=0，求解公式为

总的CO2驱动态埋存潜力为

式中，Bgi为每个网格的气体摩尔密度，mol/m3，对应Eclipse中的BGAS关键字。sGi、sgi、yi分别表示每个网格的束缚气饱和度、每个网格的含气饱和度和每个网格的CO2在气相中的摩尔百分数，分别对应Eclipse软件中的关键字SGI、SGAS、YMF1；M为总的CO2驱动态埋存潜力，kg。

3 应用示例

3.1 常规砂岩储层CO2驱动态埋存潜力计算

依据大庆榆树林油田扶余油层的数据［16］，应用CO2动态埋存潜力计算方法，研究CO2的驱油埋存过程。大庆榆树林扶余油层含油面积2.36 km2，地质储量217.8×104t，平均孔隙度10%~30%，渗透率(10~130)×10−3μm2，油藏埋深1 000~2 300 m，原始地层压力22 MPa，地层原始温度108℃，地层原油密度为0.815 g/cm3，地层原油黏度3.6 mPa · s。榆树林油田地层水总矿化度平均为4 767 mg/L，pH值7.9，地层水为NaHCO3型。

采用机理模型均质网格，模拟区域大小为190 m×110 m×70 m，孔隙度为0.3，渗透率为30×10−3μm2。共有8个组分，各组分的初始摩尔百分含量：CO2为0、C1为3.00%、C2为4.69%、C3为5.91%、C4-6为4.67%、C7+1为9.75%、C7+2为8.51%、C7+3为63.48%，采用PR状态方程进行求解［17］。原始地层压力为22 MPa，采用五点井网四注一采。注入井定井底压力30 MPa注入，设置注入50 a；生产井定井底压力10 MPa生产，生产了25 a。一共模拟了50 a用以观察CO2的埋存过程。在第25 a，生产井关井导致地层压力急剧升高，CO2注入速度逐渐降低；在第26 a，注入井达到定井底压力，从而转为关井。

采用基于组分闪蒸运算的CO2驱动态埋存潜力计算方法进行计算，预测CO2动态埋存潜力为27×104t。随着时间的增加，自由CO2、束缚CO2和总的CO2埋存质量逐渐增加，并且保持同样的变化趋势，如图3所示。第25 a，由于采油井的关闭，CO2的持续注入，导致自由CO2、束缚CO2和总的CO2埋存质量发生大幅增加。前25 a，随着采油井的采出，溶解CO2逐年减少；第25 a关井后，溶解CO2开始增加。溶解CO2随时间的质量变化曲线如图4所示。

图3 总的CO2埋存质量随时间变化Fig.3 Variation of total CO2 storage mass over the time

图4 溶解CO2质量随时间变化Fig.4 Variation of dissolved CO2 mass over the time

25~26 a中间出现小波折的原因：在前25 a开采过程中，注入到地下的CO2有很大一部分会跟随原油被开采出来。第25 a生产井关井，注入的CO2全部埋存在地下，从而束缚CO2、溶解CO2、自由CO2和总的CO2质量急速增加。由于生产井关井导致的CO2的大量注入和地层压力的急速升高，在第26 a，注入井达到定井底压力，从而转为关井停止注入CO2，总的CO2质量不再变化。各种类型的CO2埋存潜力占比随时间的变化数据见表1。

表1 不同埋存机理的CO2占比Table 1 CO2 reserve ratio of different storage mechanisms

3.2 低渗油藏水驱转CO2驱动态埋存潜力计算

中原油田卫42-43块构造为深层特低渗透油藏，构造面积约10.5 km2；油藏温度114~128℃，原始地层压力34.5~39.29 MPa；油藏含油面积5.69 km2，有效井段60 m，石油地质储量270.55×104t，平均孔隙度12.7%，平均渗透率2.25×10−3μm2；地面原油密度为0.853 4 g/cm3，地面原油黏度为8.84 mPa · s；地层水总矿化度26.9×104mg/L，氯离子含量为16.1×104mg/L，水型为CaCl2型。

采用角点网格，建立非均质网格模型，共有网格数为52×262×13=177 112，如图5所示。地层平均孔隙度为0.14，渗透率为1.9×10−3μm2。通过相态模拟得到7个拟组分，各组分的初始摩尔百分含量：N2为1.97%、CO2为0、C1为30.98%、C2为14.66%、C5+为3.77%、C10+为6.42%、C19+为42.20%，采用PR状态方程进行求解。原始地层压力为35 MPa，2口注入井4口采油井。注入井定井底压力为50 MPa注入，注入了5 a，一共注入了79.2×104t CO2；生产井定井底压力为20 MPa生产，生产了5 a；后期埋存了40 a，总的模拟时间为45 a。

采用基于组分闪蒸运算的CO2驱动态埋存潜力计算方法进行计算，预测动态埋存储量为50×104t。生产阶段(即模拟的前5 a)，自由CO2和溶解CO2都在减少，低渗透油藏经历过水驱开发，CO2的注入能力变差，注入气的速度小于产出气的速度，所以自由CO2和溶解CO2逐年减少。但是注入初期，注入的CO2容易进入到地层孔隙中，从而束缚CO2逐渐增加；后期随着注入气的速度小于产出气的速度，束缚CO2也呈现逐年减少的趋势。第5 a关井埋藏后，自由CO2减少，束缚CO2和溶解CO2增加。随着埋存时间的增加，自由CO2转化为溶解CO2和束缚CO2。第45 a时各种埋存机理的贡献值见表2，可以看出低渗油藏采收率低，地下残余油多，溶解在油中的CO2多。油藏经历过水驱开发，地层水含量高，不能忽略CO2在地层水中的溶解。

图5 研究区域实际模型非均质网格示意图(模拟45 a时的含油饱和度场分布)Fig.5 Schematic heterogeneous grid in the actual model of the study area (simulating the distribution of the oil saturation in the 45th year)

表2 第45 a不同埋存机理的CO2占比Table 2 CO2 reserve ratio of different storage mechanisms in the 45th year

4 结论

(1)基于组分闪蒸运算的CO2驱动态埋存潜力计算方法考虑了油田的生产实际和CO2驱的埋存机理，可以计算溶解CO2、束缚CO2、自由CO2和总的CO2埋存潜力。

(2)随着埋存时间的增加，自由CO2转变成束缚CO2和溶解CO2，而油藏经历过水驱开发，地层水含量高，不能忽略CO2在地层水中的溶解。

(3)基于组分闪蒸运算的CO2驱动态埋存潜力计算方法可以计算不同种类油藏不同开发方式的动态埋存潜力，方法简单实用。