0 引言
有效的CO2地质封存主要取决于地层结构长期安全封存的能力以及孔隙空间和覆盖岩层的可用性[7]。目前,中国CCUS产业发展主要分为探索研究、先导试验和工业化发展三个阶段,如表1所示[8]。但由于CO2注入过程所引起的地质力学效应可能会存在不同,且储层物性不同也会对CO2封存效果产生不同影响,因此在进行矿场试验前对潜在封存点进行系统评估是非常必要的[9,10]。罗家寨气田于2018年申报探明面积94.68 km2和地质储量776.85 × 108 m3,具有较好的储集物性且在纵向上构成了多套生、储、盖组合[11]。四川盆地油气勘探实践表明,嘉陵江组沉积厚度大、发育多套储盖组合且保存条件较好。此外,嘉陵江组具有埋藏浅、投资少、产气层段多、储集类型和圈闭类型多、保存条件较好的成藏特征。受印支期运动影响,古隆起核部嘉三-嘉五段遭受了强烈的剥蚀,早期发现的气藏主要以分布在蜀南地区以嘉一-嘉三段气藏为主。前人对泸州古隆起周缘嘉陵江组开展了油气地质研究,认为嘉二、嘉四段发育优质孔隙型碳酸盐岩储层。其中嘉二段主要为白云岩、膏质云岩及膏岩等,且在横向上分布稳定,各井岩性组合及厚度与邻区均可对比[12]。此外,目前鲜有对罗家寨气田嘉陵江组进行CO2地质封存的研究。
CO2地质封存是一个涉及热-流-固-化（thermal- hydraulic-mechanical-chemical, THMC）多场耦合问题的过程。以罗家寨气田为研究对象,主要针对CO2注入嘉陵江组后的多相流和渗流力学问题,通过FLAC3D-TOUGHREACT模拟器评估枯竭油气藏CO2地质封存的动态过程。通过对CO2注入10年间的地层压力、储层物性变化以及CO2运移分布等进行分析研究,为罗家寨气田CO2地质封存现场先导试验提供理论依据。
1 储层特征及模型
1.1 控制方程和本构模型
1.1.1 动量守恒
应力-应变本构关系是基于考虑热量传递引起的热膨胀和单/双孔隙介质的孔弹性响应,将双孔隙介质（裂缝为1,基质为2）表达为[19]：
${{\xi }_{ij}}=2G{{\varepsilon }_{ij}}+\frac{2Gv}{1-2v}{{\varepsilon }_{kk}}{{\delta }_{ij}}-\left( \alpha _{p}^{1}p+\alpha _{p}^{2}p \right){{\delta }_{ij}}-{{\alpha }_{T}}T{{\delta }_{ij}}$ (1)
式中：${{\xi }_{ij}}$为柯西应力张量的转置散度;G为剪切模量;v为泊松比;$\alpha_{p}^{(i)}$为裂缝(1)和基质(2)的Biot系数;p为孔隙流体压力;T为岩石温度;$\alpha_{T}$为热膨胀系数;$\delta_{ij}$为线性化应变,定义为位移梯度$u{i, j}$的对称部分。应力-应变关系${{\varepsilon }_{ij}}$定义为：
${{\varepsilon }_{ij}}=\frac{1}{2}\left( {{u}_{i,j}}+{{u}_{j,i}} \right)$ (2)
1.1.2 孔渗性演化
由于不考虑化学反应,因此孔隙度主要随力学变形而变化,其关系如下：
$\varphi ={{\varphi }_{\text{h}}}+\left( {{\varphi }_{0}}-{{\varphi }_{\text{h}}} \right){{\text{e}}^{-\lambda \tau '}}$ (3)
式中：${{\varphi }_{0}}$为零有效应力下的孔隙度;${{\varphi }_{h}}$为高有效应力下的孔隙度;$\tau '$为有效平均应力;$\lambda $为岩石可压缩性的常数参数,Pa-1。裂缝和基质渗透率的变化被视为与孔隙度不同的过程。水-力学响应通过有效正应力的变化来改变裂缝孔径：
${{b}_{s}}={{b}_{r}}+\left( {{b}_{\max }}-{{b}_{r}} \right)\exp \left[ -\eta \left( \sigma ' -{{\sigma }_{0}}' \right) \right]$ (4)
式中：${{b}_{s}}$为裂隙孔径,m;${{\operatorname{b}}_{\text{r}}}$为残余孔径,m;${{b}_{\max }}$为零应力下最大孔径,m;$\eta $为非线性断裂刚度,1/MPa;$\sigma '$为有效应力;$\sigma _{0}^{'}$是零变形发生时的有效应力。
除上述介绍的相关控制平衡方程以及本构方程,其他相关本构方程和数学模型可参见文献[19]。
1.2 模型建立
以罗家寨气田嘉陵江组为背景,建立热-流-固（thermal-hydraulic-mechanical, THM）耦合三维数值模型。此三维地质模型按照实际工程中CO2注入量进行设计,模型范围以注入井为中心,为1 500 m × 1 500 m × 700 m,如图1 所示。该三维数值模型的三个不同区域分别代表上覆盖层（300 m厚）、目标储层（100 m厚）和下伏盖层（300 m厚）,其中储层选取为含气层位相对较多的嘉二2亚段和嘉二1亚段,上覆盖层的上边界距离地面3 110 m。
Fig. 1 Conceptual model diagram of three-dimensional numerical model图1 三维数值模型的概念模型图 |
基于图1 的概念模型图,分别建立了三维数值模型1和模型2。模型1的网格划分为15 × 15 × 28,盖层与储层单个网格的模型尺寸为100 m × 100 m × 25 m,如图2（a） 所示。模型2则是基于模型1对储层进一步细化,上下盖层的网格划分与模型1保持一致,即盖层单个网格模型尺寸为100 m × 100 m × 25 m,储层的单个网格则为100 m × 100 m × 10 m,如图2（b） 所示。
Fig. 2 (a) 3D numerical Model 1 grid division diagram; (b) Model 2 grid division diagram图2 （a）三维数值模型1网格划分图;（b）模型2网格划分图 |
1.3 模型选取
为确定网格划分对数值模拟的影响,在相同储层参数和外部条件下,对三维数值模型1和模型2的模拟结果进行了比较。如图3 所示,模型1和模型2的模拟结果显示,两者注入井的井筒压力和CO2饱和度的变化趋势相似,且最终的井筒压力和CO2饱和度基本趋于一致。这说明网格划分方式对模拟结果的影响不大。因此为了缩短模拟计算时间,选取模型1作为本次数值模拟的三维耦合数值模型。
Fig. 3 Comparison of wellbore pressure and CO2 saturation simulation results of injection wells in Model 1 and Model 2图3 模型1与模型2的注入井井筒压力和CO2饱和度模拟结果对比 |
1.4 模型参数
Table 2 Simulation parameters of 3D numerical Model 1表2 三维数值模型1的模拟参数表 |
| 参数 | 数值 |
|---|---|
| 盖层体积模量/GPa | 37.1 |
| 盖层剪切模量/GPa | 19.1 |
| 盖层抗拉强度/MPa | 12.5 |
| 盖层内聚力/MPa | 4 |
| 储层体积模量/GPa | 66.7 |
| 储层剪切模量/GPa | 25.0 |
| 储层抗拉强度/MPa | 20 |
| 储层内聚力/MPa | 8 |
| 泊松比 | 0.28 |
| z方向垂直应力梯度/(Pa/m) | 30 000 |
| y方向水平应力梯度/(Pa/m) | 20 000 |
| x方向水平应力梯度/(Pa/m) | 50 000 |
| 盖层渗透率/(×10-16 m2) | 0.1 |
| 盖层孔隙度/% | 1 |
| 储层渗透率/(×10-16 m2) | 8.65 |
| 储层孔隙度/% | 3.09 |
| 井筒初始压力/MPa | 34.1 |
| CO2注入速率/(kg/s) | 0.1 |
2 结果与讨论
本次模拟研究在中心网格（750, 750, 350）处设置了一口垂直注入井,向储层连续注入10年的CO2。为了更直观地观察CO2在储层中的运移规律及封存机理,在注入井所在位置沿（x, y）和（x, z）方向分别截取剖面（以下简称y平面和z平面）,通过不同方向上二维平面的变化特征来反映三维地质体的空间分布特征。
2.1 地层及注入井压力变化
从地层压力场分布变化特征可以看出,随着CO2的注入,地层原有的压力平衡被打破,注入井区域附近的压力快速上升。在y平面上分布的储层压力均随着CO2的注入而升高,如图4（a1 ~ d1） 所示。但就z平面上的压力分布而言,上覆盖层的压力因CO2的浮力作用逐渐上升,而下伏盖层的压力则是先下降后升高,如图4（a2 ~ d2） 所示。随着注入时间的推移,高压场的影响范围以注入井为核心向四周逐渐扩大,直至CO2封存结束整个地层压力也随之增大。
Fig. 4 Spatial distribution evolution of formation pressure during CO2 geological storage图4 CO2地质封存期间地层压力空间分布演变 |
从图5 中可以观察到,注入井的井筒压力先是迅速上升,这可归因于储层物性较差使得CO2注入初期主要集中在注入井附近区域。随着CO2的持续注入导致储层的孔渗性增大,使得CO2以注入井为中心在压差驱动下逐渐向四周扩散。因此,注入井井筒压力虽在A点（t = 92 d）之后仍直线上升至约39 MPa,但其上升速率（0.001 MPa/d）相较于A点之前（0.007 MPa/d）却呈现下降趋势。
Fig. 5 Evolution curve of the wellbore pressure in injection well图5 注入井井筒压力演变曲线 |
2.2 CO2运移分布及封存总量
由于中间储层的初始孔隙度和渗透率明显高于上下盖层,因此CO2主要在中间储层中运移。在进行碳封存的过程中,CO2既有气相又有水相,且气相与水相处于平衡状态,如方程所示：
$C{{O}_{2}}(sc)\rightleftarrows C{{O}_{2}}(aq)$ (5)
$C{{O}_{2}}(aq)+{{H}_{2}}O\to HCO_{3}^{-}+{{H}^{+}}$ (6)
式中：sc代表超临界状态;aq代表水相。上述反应式表明CO2除了被束缚于孔隙空间内以实现构造封存（图6 ）,还可通过生成碳酸实现CO2的溶解封存（图7 ）。
Fig. 6 Spatial distribution evolution of CO2 gas mole fraction during geological storage图6 地质封存期间CO2气体摩尔分数空间分布演化 |
Fig. 7 Spatial distribution evolution of aqueous CO2 concentration during geological storage图7 地质封存期间水相中CO2浓度空间演化 |
由于CO2的运移主要是在压力梯度、浓度差和密度差即浮力的作用下进行的。因此,随着CO2的注入,CO2在压力梯度和浓度差的联合作用下在y平面上[图6（a1 ~ d1） ]向四周扩散,进行水平方向上的运移。此外,CO2和地层水之间的密度差使得CO2在浮力的作用下在z平面发生垂向运移。由于注入初期静态封存就开始作用,CO2在地层中运移时间越久,CO2持续溶解于地层水和被束缚于孔隙中,因此注入井附近区域的CO2气体浓度和液相CO2浓度均持续增加。从图6 中还可以观察到,高浓度液相CO2更偏向于存在于储层下方（z平面）,这是由于饱和CO2的地层水密度相较于未饱和地层水更大。
从图6（c2） 和图7（c2） 均可观察到,在CO2羽流占据中间储层大部分空间后,气相以及液相CO2开始慢慢扩散至上下盖层中。但由于上下盖层的低孔渗性,直至CO2注入10年后仍未突破盖层,且侵入盖层的高度很小,见图6（d2） 和图7（d2） 。这说明所选层位具有良好的连续性、完整性和封闭性,足够承受羽状流产生的浮力,可保证封存于地层中的大部分CO2处于盖层下方。
由于暂不考虑CO2矿化封存量,因此最终所获得的总CO2地质封存潜力可认为是构造封存和溶解封存两种机理下的地质封存潜力之和。
（1）构造封存。在构造地层圈闭和水动力圈闭作用下,地层构造捕获的CO2封存潜力计算公式如下：
${{m}_{C{{O}_{2ts}}}}=\int{\int\limits_{V}{\int{\phi \times {{\rho }_{C{{O}_{2}}}}\times \left( 1-{{S}_{w}} \right)\text{d}V}}}$ (7)
式中：${{m}_{C{{O}_{2ts}}}}$为构造封存的CO2潜力,kg;$\phi $为储层的孔隙度,%;${{\rho }_{C{{O}_{2}}}}$为储层中CO2的密度,kg/m3;${{S}_{w}}$为残余水饱和度,即CO2排水后岩石孔隙中的水的体积与孔隙体积的百分比,%;V为网格体积,m3。
（2）溶解封存。CO2溶解到地层水中的封存量可视为原始地层水达到CO2饱和时,所能溶解的CO2量,计算公式如下：
${{m}_{C{{O}_{2td}}}}=\int{\int\limits_{V}{\int{\phi \times {{R}_{C{{O}_{2}}}}\times {{S}_{\text{w}}}\times {{M}_{C{{O}_{2}}}}\times {{\rho }_{\text{w}}}\text{d}V}}}$ (8)
式中：${{m}_{C{{O}_{2td}}}}$为溶解捕获地质封存潜力,kg;${{R}_{C{{O}_{2}}}}$为CO2在地层水中的溶解度,mol/kg;${{M}_{C{{O}_{2}}}}$为CO2的摩尔质量,0.044 kg/mol;${{\rho }_{\text{w}}}$为初始地层水的密度,kg/m3。
基于上述构造与溶解封存量的计算模型,通过数值模拟获得了CO2在地层中的构造封存量和溶解封存量,分别为1.510 × 107 kg和1.701 × 107 kg,即最终获得了3.211 × 107 kg的总封存量。
2.3 储盖层孔渗性演化
如图8 所示,在CO2地质封存期间,中间储层的孔隙度随着CO2的注入逐渐增大。此外,从图8（a1 ~ d1） （y平面）可以观察到注入井区域附近的孔隙度上升幅度更为明显,而孔隙度的增大可加快地层流体的流动,进而提高CO2被束缚封存的效率,这也是注入井区域附近的CO2浓度更高的原因。同样,CO2的注入可提高储层渗透率（图9 ）,进而提升流体的渗流能力,且在注入井附近区域的渗透率增大速率相对更为明显。这可归因于CO2注入期间该区域内的高孔隙压力。而上下盖层由于良好的连续性和封闭性,使得CO2羽流主要在储层中运移,进而导致盖层的孔渗性在整个碳封存期间未发生明显变化。
Fig. 8 Reservoir porosity evolution during CO2 geological storage图8 CO2地质封存期间储层孔隙度演化 |
Fig. 9 Reservoir permeability evolution during CO2 geological storage图9 CO2地质封存期间储层渗透率演变 |
如图10 所示,地层平均渗透率与孔隙度在整个CO2地质封存期间均逐渐增大。这是由于CO2溶于水后呈弱酸性,地层中的碳酸盐矿物与酸性溶液发生反应,易使岩石产生大量的次生孔隙甚至微裂缝。且在动力学作用下,溶解于水中的CO2易沿产生的孔隙或微裂缝路径流动,也会导致孔隙度和渗透率进一步增大。相较于孔隙度在整个碳封存期间匀速增大,渗透率的增大主要集中于碳封存前期（1年内）。这可能是由于前期CO2因储层孔渗性较差,CO2主要集中在井筒附近,导致孔隙压力快速升高（图5 ）,进而使得地层产生了较大的微裂缝,以及提高了孔隙间的连通性。
Fig. 10 Evolution of average permeability and porosity during CO2 geological storage图10 CO2地质封存期间地层平均渗透率与孔隙度的演变 |
3 结论
基于罗家寨气田建立地质模型,针对CO2注入地层后的多相流和渗流力学问题,利用FLAC3D-TOUGHREACT软件模拟评估了注入CO2后气藏的地层压力、储层物性变化和CO2运移规律等,得到以下结论：
（1）在CO2注入的10年间,因CO2注入导致的高压场以注入井为中心向四周逐渐扩大,直至封存结束整个地层压力均呈现增大的状态。
（2）CO2主要以气相构造封存于孔隙中和通过生成碳酸溶解封存于地层水中。由于盖层的低渗性,CO2羽流主要在中间储层运移,未发生突破盖层的情况。此外,通过数值模拟研究发现,在注入速率为0.1 kg/s的情况下,该地层的CO2总封存量为3.211 × 107 kg。
（3）储层的孔隙度和渗透率随着CO2的注入逐渐增大,而盖层由于良好的连续性和封闭性,其物性未发生明显变化。此外,由于高压区主要集中在注入井附近,导致该区域孔渗性增长速率相对更快。
