数值模拟基于已有THMC多场耦合模拟器FLAC3D-TOUGHREACT, 该模拟软件可表征多相流体在含有离散裂缝网络多孔介质中流动规律, 同时软件可动态预测岩石流体特性与力学特性的参数演化, 包括渗透率、孔隙度、力学模量的数值演化, 目前已有学者对该模拟软件进行了验证^{[17, 18]}。该模拟器主要运用孔弹性力学、渗流力学、地球化学等方程描述改变岩石力学和流体性质的机理过程, 下面将介绍相关控制平衡方程及其对应的本构方程。

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)}$为裂缝⁽¹⁾和基质⁽²⁾的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]。

以罗家寨气田嘉陵江组为背景, 建立热-流-固（thermal-hydraulic-mechanical, THM）耦合三维数值模型。此三维地质模型按照实际工程中CO₂注入量进行设计, 模型范围以注入井为中心, 为1 500 m × 1 500 m × 700 m, 如图1所示。该三维数值模型的三个不同区域分别代表上覆盖层（300 m厚）、目标储层（100 m厚）和下伏盖层（300 m厚）, 其中储层选取为含气层位相对较多的嘉二²亚段和嘉二¹亚段, 上覆盖层的上边界距离地面3 110 m。

图1

Fig. 1

	Figure Option ViewDownloadNew Window
	图1 三维数值模型的概念模型图Fig. 1 Conceptual model diagram of three-dimensional numerical model

基于图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）所示。

图2

Fig. 2

	Figure Option ViewDownloadNew Window
	图2 （a）三维数值模型1网格划分图; （b）模型2网格划分图Fig. 2 (a) 3D numerical Model 1 grid division diagram; (b) Model 2 grid division diagram

为确定网格划分对数值模拟的影响, 在相同储层参数和外部条件下, 对三维数值模型1和模型2的模拟结果进行了比较。如图3所示, 模型1和模型2的模拟结果显示, 两者注入井的井筒压力和CO₂饱和度的变化趋势相似, 且最终的井筒压力和CO₂饱和度基本趋于一致。这说明网格划分方式对模拟结果的影响不大。因此为了缩短模拟计算时间, 选取模型1作为本次数值模拟的三维耦合数值模型。

图3

Fig. 3

	Figure Option ViewDownloadNew Window
	图3 模型1与模型2的注入井井筒压力和CO2饱和度模拟结果对比Fig. 3 Comparison of wellbore pressure and CO2 saturation simulation results of injection wells in Model 1 and Model 2

三维地质模型1中所使用的物理参数主要根据罗家寨气田的测井数据和相关文献^{[20, 21, 22, 23]}确定, 见表2。地层的初始孔隙压力分布为9.8 × 10³ Pa/m, 最高地层温度为83β ℃, 温度梯度为0.024 6β ℃/m。下伏围岩层初始主应力状态在x、y、z三个方向分别为99 MPa、51 MPa和81 MPa。

表2

Table 2

表2(Table 2)

表2 三维数值模型1的模拟参数表 Table 2 Simulation parameters of 3D numerical Model 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-16m2) 0.1 盖层孔隙度/% 1 储层渗透率/(× 10-16m2) 8.65 储层孔隙度/% 3.09 井筒初始压力/MPa 34.1 CO2注入速率/(kg/s) 0.1

表2 三维数值模型1的模拟参数表 Table 2 Simulation parameters of 3D numerical Model 1

从地层压力场分布变化特征可以看出, 随着CO₂的注入, 地层原有的压力平衡被打破, 注入井区域附近的压力快速上升。在y平面上分布的储层压力均随着CO₂的注入而升高, 如图4（a1 ~ d1）所示。但就z平面上的压力分布而言, 上覆盖层的压力因CO₂的浮力作用逐渐上升, 而下伏盖层的压力则是先下降后升高, 如图4（a2 ~ d2）所示。随着注入时间的推移, 高压场的影响范围以注入井为核心向四周逐渐扩大, 直至CO₂封存结束整个地层压力也随之增大。

图4

Fig. 4

	Figure Option ViewDownloadNew Window
	图4 CO2地质封存期间地层压力空间分布演变Fig. 4 Spatial distribution evolution of formation pressure during CO2 geological storage

从图5中可以观察到, 注入井的井筒压力先是迅速上升, 这可归因于储层物性较差使得CO₂注入初期主要集中在注入井附近区域。随着CO₂的持续注入导致储层的孔渗性增大, 使得CO₂以注入井为中心在压差驱动下逐渐向四周扩散。因此, 注入井井筒压力虽在A点（t = 92 d）之后仍直线上升至约39 MPa, 但其上升速率（0.001 MPa/d）相较于A点之前（0.007 MPa/d）却呈现下降趋势。

图5

Fig. 5

	Figure Option ViewDownloadNew Window
	图5 注入井井筒压力演变曲线Fig. 5 Evolution curve of the wellbore pressure in injection well

由于中间储层的初始孔隙度和渗透率明显高于上下盖层, 因此CO₂主要在中间储层中运移。在进行碳封存的过程中, CO₂既有气相又有水相, 且气相与水相处于平衡状态, 如方程所示：

$C{{O}_{2}}(sc)\rightleftarrows C{{O}_{2}}(aq)$ (5)

$C{{O}_{2}}(aq)+{{H}_{2}}O\to HCO_{3}^{-}+{{H}^{+}}$ (6)

式中：sc代表超临界状态; aq代表水相。上述反应式表明CO₂除了被束缚于孔隙空间内以实现构造封存（图6）, 还可通过生成碳酸实现CO₂的溶解封存（图7）。

图6

Fig. 6

	Figure Option ViewDownloadNew Window
	图6 地质封存期间CO2气体摩尔分数空间分布演化Fig. 6 Spatial distribution evolution of CO2 gas mole fraction during geological storage

图7

Fig. 7

	Figure Option ViewDownloadNew Window
	图7 地质封存期间水相中CO2浓度空间演化Fig. 7 Spatial distribution evolution of aqueous CO2 concentration during geological storage

由于CO₂的运移主要是在压力梯度、浓度差和密度差即浮力的作用下进行的。因此, 随着CO₂的注入, CO₂在压力梯度和浓度差的联合作用下在y平面上[图6（a1 ~ d1）]向四周扩散, 进行水平方向上的运移。此外, CO₂和地层水之间的密度差使得CO₂在浮力的作用下在z平面发生垂向运移。由于注入初期静态封存就开始作用, CO₂在地层中运移时间越久, CO₂持续溶解于地层水和被束缚于孔隙中, 因此注入井附近区域的CO₂气体浓度和液相CO₂浓度均持续增加。从图6中还可以观察到, 高浓度液相CO₂更偏向于存在于储层下方（z平面）, 这是由于饱和CO₂的地层水密度相较于未饱和地层水更大。

从图6（c2）和图7（c2）均可观察到, 在CO₂羽流占据中间储层大部分空间后, 气相以及液相CO₂开始慢慢扩散至上下盖层中。但由于上下盖层的低孔渗性, 直至CO₂注入10年后仍未突破盖层, 且侵入盖层的高度很小, 见图6（d2）和图7（d2）。这说明所选层位具有良好的连续性、完整性和封闭性, 足够承受羽状流产生的浮力, 可保证封存于地层中的大部分CO₂处于盖层下方。

由于暂不考虑CO₂矿化封存量, 因此最终所获得的总CO₂地质封存潜力可认为是构造封存和溶解封存两种机理下的地质封存潜力之和。

（1）构造封存。在构造地层圈闭和水动力圈闭作用下, 地层构造捕获的CO₂封存潜力计算公式如下：

${{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}}}}$为构造封存的CO₂潜力, kg; $\phi $为储层的孔隙度, %; ${{\rho }_{C{{O}_{2}}}}$为储层中CO₂的密度, kg/m³; ${{S}_{w}}$为残余水饱和度, 即CO₂排水后岩石孔隙中的水的体积与孔隙体积的百分比, %; V为网格体积, m³。

（2）溶解封存。CO₂溶解到地层水中的封存量可视为原始地层水达到CO₂饱和时, 所能溶解的CO₂量, 计算公式如下：

${{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}}}}$为CO₂在地层水中的溶解度, mol/kg; ${{M}_{C{{O}_{2}}}}$为CO₂的摩尔质量, 0.044 kg/mol; ${{\rho }_{\text{w}}}$为初始地层水的密度, kg/m³。

基于上述构造与溶解封存量的计算模型, 通过数值模拟获得了CO₂在地层中的构造封存量和溶解封存量, 分别为1.510 × 10⁷ kg和1.701 × 10⁷ kg, 即最终获得了3.211 × 10⁷ kg的总封存量。

如图8所示, 在CO₂地质封存期间, 中间储层的孔隙度随着CO₂的注入逐渐增大。此外, 从图8（a1 ~ d1）（y平面）可以观察到注入井区域附近的孔隙度上升幅度更为明显, 而孔隙度的增大可加快地层流体的流动, 进而提高CO₂被束缚封存的效率, 这也是注入井区域附近的CO₂浓度更高的原因。同样, CO₂的注入可提高储层渗透率（图9）, 进而提升流体的渗流能力, 且在注入井附近区域的渗透率增大速率相对更为明显。这可归因于CO₂注入期间该区域内的高孔隙压力。而上下盖层由于良好的连续性和封闭性, 使得CO₂羽流主要在储层中运移, 进而导致盖层的孔渗性在整个碳封存期间未发生明显变化。

图8

Fig. 8

	Figure Option ViewDownloadNew Window
	图8 CO2地质封存期间储层孔隙度演化Fig. 8 Reservoir porosity evolution during CO2 geological storage

图9

Fig. 9

	Figure Option ViewDownloadNew Window
	图9 CO2地质封存期间储层渗透率演变Fig. 9 Reservoir permeability evolution during CO2 geological storage

如图10所示, 地层平均渗透率与孔隙度在整个CO₂地质封存期间均逐渐增大。这是由于CO₂溶于水后呈弱酸性, 地层中的碳酸盐矿物与酸性溶液发生反应, 易使岩石产生大量的次生孔隙甚至微裂缝。且在动力学作用下, 溶解于水中的CO₂易沿产生的孔隙或微裂缝路径流动, 也会导致孔隙度和渗透率进一步增大。相较于孔隙度在整个碳封存期间匀速增大, 渗透率的增大主要集中于碳封存前期（1年内）。这可能是由于前期CO₂因储层孔渗性较差, CO₂主要集中在井筒附近, 导致孔隙压力快速升高（图5）, 进而使得地层产生了较大的微裂缝, 以及提高了孔隙间的连通性。

图10

Fig. 10

	Figure Option ViewDownloadNew Window
	图10 CO2地质封存期间地层平均渗透率与孔隙度的演变Fig. 10 Evolution of average permeability and porosity during CO2 geological storage