CO₂地质封存中随机裂隙网络走向对盖层密封性影响

Effect of Random Fracture Network Orientations on Sealing Performance of Caprock in CO₂ Geological Sequestration

摘要

在CO₂地质封存过程中，为探究盖层中裂隙网络走向对CO₂-咸水两相运移的影响规律，进而评价盖层密封性，本文数值模拟选用了6种不同走向的裂隙网络（其角度分别为0°~180°、30°~150°、45°~135°、60°~120°、90°~90°和90°~180°），实现显式裂隙网络下CO₂驱替咸水的两相流研究.研究发现裂隙网络走向会直接改变液相饱和度的赋存情况，从而影响CO₂-咸水驱替难度；当CO₂在盖层中达到相同渗透深度时，驱替时间随裂隙倾角（0°~60°）降低了12.59倍，但倾角增大到60°后不再有明显影响；随着裂隙网络渗透范围的扩大，CO₂在盖层中的渗透量增加.其中，正交裂隙（90°~180°）的渗透量最大.因此，考虑裂隙网络走向对于盖层密封性的评价具有重要意义.

Abstract

In the process of CO₂ geological storage, in order to explore the influence of fracture network orientations on the two-phase transport of CO₂-salt water in the caprock and evaluate the sealing performance of the caprock, six fracture networks with different orientations (with angles of 0°-180°, 30°-150°, 45°-135°, 60°-120°, 90°-90° and 90°-180°, respectively) were selected for numerical simulation, so as to realize the two-phase flow study of CO₂ displacement of salt water under explicit fracture network. It is found that the fracture network orientations directly change the occurrence of water saturation, thereby affecting the difficulty of CO₂-salt water displacement. When CO₂ reaches the same penetration depth in the caprock, the displacement time decreased by 12.59 times with the fracture angle (0°-60°), but there is no obvious effect after the inclination angle increased to 60°. With the expansion of the permeability range of the fracture network, the permeability of CO₂ in the caprock increases. Among them, orthogonal fractures (90°-180°) have the largest permeability. Therefore, considering the fracture network direction is of great significance for the caprock sealing evaluation.

Graphical abstract

关键词

CO₂地质封存 / 盖层密封性 / 裂隙走向 / 离散裂隙基质模型 / 数值模拟.

Key words

CO₂ geological sequestration / caprock sealing / orientation of fracture / discrete fracture matrix model / numerical simulation

0 引言

根据国际气候科学机构“全球碳预算”（Global Carbon Project）的最新数据，全球CO₂排放量在2022年达到创纪录高位后，如果各国没有大幅减少排放，将在未来9年内耗尽剩余的“碳排放”预算，从而导致全球变暖突破1.5 ℃，引发极端天气频繁出现（王卓妮等，2022）.因此，作为缓解大气中CO₂含量的有效途径之一的CO₂资源化利用和地质封存，即CCUS（Carbon Dioxide Capture， Utilization and Storage）研究迫在眉睫.近些年来，国内外已经有多个CCUS示范工程相继启动，例如挪威的Sleipne、阿尔及利亚的In Salah以及中国鄂尔多斯盆地神华咸水层封存项目.目前全球最大的千万吨级别CCUS产业集群为“休斯顿航道CCUS创新区”，预计建成后可每年储存1亿吨CO₂（邓一荣等，2023；蒋恕等，2023；赵震宇等，2023）.但在研究过程中发现，CO₂长期封存在地下时若发生泄漏不仅会引发显著的气候变化，而且会引发地下水污染、土壤生态系统损害等问题（陈浮等，2016），可能对附近人类的健康和生命造成威胁.因此，盖层密封性的研究是CO₂地质封存安全的关键科学问题.

在CO₂地质封存中，部分CO₂在上浮力的作用下积累在盖层底部，低渗透性的膏岩、盐岩和泥页岩有着良好的密封性，能够在毛细圈闭的作用下阻止或减缓CO₂驱替咸水的两相流过程（Song and Zhang，2013；陈博文等，2023）.相关研究表明，当盖层渗透率小于0.08×10^-3 μm²时，通常不考虑盖层中多尺度孔隙间的连通性甚至裂隙性质，认为CO₂突破压力不会对盖层密封性产生明显的影响（武志德，2011）.周雪峰等（2018）认为影响盖层密封性的因素有盖层厚度、超压系数和脆性指数等；代龙等（2019）对盖层的岩石学特征和孔喉特征进行研究后发现泥岩盖层中多为微孔喉和超微孔喉，气体不易渗流，盖层密封性较好；Shang et al.（2022）研究发现，盖层的密封效率与页岩基质的脱水和CO₂吸附引起的膨胀有着显著的关系；因为CO₂‒咸水‒云母的接触角和CO₂⁃盐水界面张力会影响盖层的毛细密封效率，Arif et al. （2016）探究了不同温度、压力和盐度条件下的接触角和界面张力，发现压力越低、温度越高，密封效果越好；Adu⁃Gyamfi et al.（2022）探究了地层矿物的沉淀和溶解以及应力变化对盖层完整性的影响，发现盖层中虽然发生了岩石物理性质的变化和CO₂⁃水‒岩的相互作用，但盖层密封性良好.但事实上，持续注入CO₂会打破盖层‒储层系统的伪稳态或局部平衡条件，盖层内现存的原生裂隙可能会成为优势渗流通道.此时，如果仍忽略裂隙的影响，将难以真实反映CO₂驱替咸水的两相流渗透前端，从而高估盖层的密封性.

当CO₂在具有单一垂直裂隙的盖层中迁移时，郭兵（2020）将孔隙度和渗透率的改变量作为评价盖层密封性指标，最终发现大部分CO₂会沿着裂隙运移，裂隙处的渗透率和孔隙率均会增大，盖层的泄漏风险增大；Olabode and Radonjic（2017）通过实验模拟了原位应力条件下的矿物/沉淀反应，根据裂隙的导电性判断溶解矿物的反应性运移可以阻塞微裂隙的流动路径，从而在CO₂封存条件下提高盖层密封性，降低泄漏风险；Yin et al.（2017）对单裂隙盖层模型进行了耦合T⁃H⁃M（Thermo⁃Hydro⁃Mechanical）效应的数值计算，发现了裂隙渗透率会随裂隙倾角、覆盖层压力、初始裂隙孔径和流体压差的增大而增大.其中，裂隙的渗透率随裂隙倾角（45°~90°）增大了1.74倍，盖层的密封性能随着裂隙倾角增大而变弱.显而易见，裂隙的走向对盖层的渗透性和两相流的运移特征具有明显的影响.然而，目前主要集中在单一裂隙走向和简化模型的研究，对于复杂的裂隙网络和实际岩石盖层的研究还相对较少.Hyman et al.（2020）为探究水力压裂和结构非均质性对超临界CO₂通过盖层的影响，生成了具有不同裂隙网络密度的多个独立、相同分布的网络，模拟分析表明裂隙网络的存在对超临界CO₂驱水量的控制作用大于局部水力特性的变化；Verdon et al.（2013）对比观察了3个大型CCS（Carbon Dioxide Captureand Storage）站点的地质力学行为，发现在In Salah油田处CO₂的注入刺激了一个延伸到盖层100~200 m的裂隙网络；Currie（2017）通过数值模拟发现，与注入相关的裂隙网络会显著降低盖层封存流体的能力，甚至可能会损害储层.

由此可见，以往对于盖层密封性的研究主要是针对均质盖层，而实际盖层大多会表现出各向异性，即在研究过程中常常忽略了盖层中裂隙走向的作用，尤其是裂隙网络对多孔介质中不混溶流体两相流运移的影响，从而高估页岩盖层的密封性.除此之外，在评价盖层密封性时，大多是通过间接参数表征（如孔隙率、渗透率等），未能从裂隙网络的几何特征角度进行考虑.因此，从几何特征角度研究裂隙页岩中CO₂驱替咸水的渗透规律将能更准确地评价页岩盖层的密封性.

根据裂隙特点，本文将采用离散裂隙基质模型即通过显式部分裂隙来模拟裂隙岩体中的渗流过程.它能够辨识裂隙及基质对整体渗流过程的影响，考虑裂隙网络的复杂性，更贴合实际情况.基于此，本文以深部CO₂封存盖层为研究对象，提出考虑离散裂隙网络中CO₂驱替咸水两相流的数值模拟方法，分析地质封存时裂隙盖层中的CO₂渗透过程，探讨不同走向的离散裂隙网络对CO₂驱替前端和密封效率的影响.

1 数学模型

本文模型耦合了多个物理场间的相互作用，其中包括盖层的基质变形、裂隙中的两相流流动以及基质中的两相流流动.图1具体描述了各物理场间的耦合作用.在CO₂驱替咸水的过程中，两相流的流动会引起基质和裂隙中气液两相压力的变化，变化的压力会导致盖层基质的变形.对于基质而言，基质岩体的变形会影响CO₂的吸附/解吸，进一步影响基质中气液两相流的流动.对于裂隙，基质岩体的变形会改变其孔隙度和渗透率，从而影响裂隙中的两相流流动.与此同时，基质和裂隙中的两相流可以通过渗透作用互换物质，从而改变CO₂的分布和驱替效果.因此，CO₂在整个盖层中的渗透量会因为盖层基质的变形和两相流流动而变化.

1.1　基本假设

a.为易于理解和分析，用直线代表天然裂隙，未考虑裂隙迂曲度带来的影响；b.气体的吸附服从Langmuir吸附；c.流动过程中，没有诱导裂隙的产生，且天然裂隙的物理性质均保持不变；d.基质系统中的流动服从达西定律，裂隙中的流动满足立方定律；e.本模型暂未考虑温度场和化学场的影响.

1.2　基岩渗流方程

储层中的超临界CO₂经过长时间的水平迁移后，会在浮力的作用下累积在上覆盖层的底部.作为封存CO₂的天然屏障，盖层的密封性直接决定了CO₂封存的效率.初始条件下，咸水作为湿润相，占据了岩石的孔隙和裂隙空间；而注入的CO₂作为非湿润相，在毛细压力的作用下，渗透进入盖层，形成了气‒水驱替的两相流前端.根据质量守恒定律，咸水和CO₂的输运方程可以分别表示为：

咸水： \frac{\partial \left({\rho }_{\mathrm{w}}\varphi s_{\mathrm{w}}\right)}{\partial t}+\nabla \cdot \left(-{\rho }_{\mathrm{w}}\frac{k{k}_{\mathrm{r}\mathrm{w}}}{{\mu }_{\mathrm{w}}}\left(\nabla p_{\mathrm{w}}-{\rho }_{\mathrm{w}}\mathrm{g}\nabla H\right)\right)=Q_{\mathrm{w}}，（1）

CO₂： \frac{\partial m}{\partial t}+\nabla \cdot \left(-{\rho }_{\mathrm{g}}\frac{k{k}_{\mathrm{r}\mathrm{g}}}{{\mu }_{\mathrm{g}}}\left(\nabla p_{\mathrm{g}}-{\rho }_{\mathrm{g}}\mathrm{g}\nabla H\right)\right)=Q_{\mathrm{g}}，（2）

式中：ρ_w和ρ_g分别是咸水和CO₂的密度，kg/m³；p_w和p_g分别是咸水和CO₂的压力，MPa；k是盖层的绝对渗透率，m²；φ是盖层的孔隙度，受体积应变、孔隙压力和吸附应变的共同影响，无量纲；s_w代表了液相饱和度，无量纲；k_rw和k_rg分别是液相和气相的相对渗透率，无量纲；μ_w和μ_g分别表示液相和气相的粘度，Pa·s；H表示盖层的埋藏深度，m；g是重力加速度，m/s²；Q_w和Q_g分别代表了咸水和CO₂的源项，kg/（m³·s）.

在两相流的驱替过程中，液相饱和度s_w会随着毛细压力p_c不断变化.其中：

p_{\mathrm{w}}=p_{\mathrm{g}}-p_{\mathrm{c}},（3）

根据Brooks⁃Corey模型，毛细压力p_c可以用液相相对饱和度s_e表示为：

p_{\mathrm{c}}=p_{\mathrm{e}}(s_{\mathrm{e}}{)}^{-\mathrm{1}/\lambda }=p_{\mathrm{e}}{(\frac{s_{\mathrm{w}}-s_{\mathrm{r}\mathrm{w}}}{\mathrm{1}-s_{\mathrm{r}\mathrm{g}}-s_{\mathrm{r}\mathrm{w}}})}^{-\mathrm{1}/\lambda }，（4）

式中：p_e是毛细进入压力，取值随岩石种类不同而变化，MPa；λ是孔隙大小的分布系数，无量纲；s_rw和s_rg分别是液相和气相的残余饱和度，取决于地层岩石的孔隙尺寸分布，无量纲.

咸水在孔隙介质中自由流动，而大部分CO₂以自由态存在于孔隙中，部分CO₂以吸附态的形式存在于基质孔隙表面.本文忽略CO₂在孔隙和基质表面的扩散时间，因此CO₂的质量m_g可以表示为：

m_{\mathrm{g}}={\rho }_{\mathrm{g}}\varphi s_{\mathrm{g}}+(\mathrm{1}-\varphi ){\rho }_{\mathrm{g}\mathrm{a}}{\rho }_{\mathrm{c}}\frac{V_{\mathrm{L}}{p}_{\mathrm{g}}}{p_{\mathrm{g}}+P_{\mathrm{L}}},（5）

式中：s_g是CO₂的饱和度，无量纲；ρ_ga是CO₂在标准状态下的密度，kg/m³；ρ_c是盖层岩石的密度，kg/m³；P_L是Langmuir吸附压力，MPa；V_L是Langmuir吸附体积，m³/kg.

将方程（3）~（5）代入方程（1）和（2）中，可以得到两相流控制方程最终表达式如下：

{\rho }_{\mathrm{w}}\varphi \frac{\partial s_{\mathrm{w}}}{\partial t}+\nabla \cdot (-{\rho }_{\mathrm{w}}\frac{k{k}_{\mathrm{r}\mathrm{w}}}{{\mu }_{\mathrm{w}}}(\nabla P_{\mathrm{g}}-\frac{\partial P_{\mathrm{c}}}{\partial s_{\mathrm{w}}}\nabla s_{\mathrm{w}}-{\rho }_{\mathrm{w}}\mathrm{g}\nabla H))=Q_{\mathrm{w}}-{\rho }_{\mathrm{w}}{s}_{\mathrm{w}}\frac{\partial \varphi }{\partial t},（6）

\begin{array}{l}\left(\varphi s_{\mathrm{g}}\beta +\left(\mathrm{1}-\varphi \right){\rho }_{\mathrm{g}\mathrm{a}}{\rho }_{\mathrm{c}}\frac{V_{\mathrm{L}}{P}_{\mathrm{L}}}{{\left(p_{\mathrm{g}}+P_{\mathrm{L}}\right)}^{\mathrm{2}}}\right)\frac{\partial p_{\mathrm{g}}}{\partial t}-\varphi {\rho }_{\mathrm{g}}\frac{\partial s_{\mathrm{w}}}{\partial t}+\\ \nabla \cdot \left(-{\rho }_{\mathrm{g}}\frac{k{k}_{\mathrm{r}\mathrm{g}}}{{\mu }_{\mathrm{g}}}\left(\nabla p_{\mathrm{g}}-{\rho }_{\mathrm{g}}\mathrm{g}\nabla H\right)\right)=\\ Q_{\mathrm{g}}-\left({\rho }_{\mathrm{g}}{s}_{\mathrm{g}}-{\rho }_{\mathrm{g}\mathrm{a}}{\rho }_{\mathrm{c}}\frac{V_{\mathrm{L}}{P}_{\mathrm{g}}}{p_{\mathrm{g}}+P_{\mathrm{L}}}\right)\frac{\partial \varphi }{\partial t}\end{array},（7）

其中：

\beta =\frac{M}{Z\mathrm{R}T},（8）

式中：M为CO₂的摩尔质量，kg/mol；Z为压缩系数，无量纲；R为摩尔气体常数，J/（mol·K）；T为温度，K.

最终，将式（6）和式（7）在COMSOL实体单元中的PDE模块实现，自变量为CO₂的压力p_g和液相饱和度s_w.在之前的研究中，曾将CO₂压力和水压力看作所求偏微分方程的自变量（Raduha et al.，2016），但这样的求解过程中存在两个不能回避的问题：（1）孔隙介质中初始水饱和的状态无法用初始压力很好地表示，即用压力控制CO₂驱替咸水的过程，需要一个较大的CO₂压力初值；（2）CO₂压力和水压力为自变量的控制方程，导致驱替的前缘会出现压力不稳定的现象，两相流求解时收敛性较差.选用CO₂压力和水饱和度为自变量，在PDE模块中求解，与COMSOL最新版本中两相流模块的控制方程类似，能够增加求解效率，提高收敛性.

1.3　裂隙渗流方程

在CO₂驱替咸水的过程中，将显式裂隙视为渗流区域的内边界，求解裂隙流时可通过对所有裂隙边界增加弱形式的方式表现，其渗流守恒方程为：

咸水： d_{\mathrm{f}}\frac{\partial }{\partial t}({\varphi }_{\mathrm{f}}{\rho }_{\mathrm{w}})+{\nabla }_{\mathrm{T}}\cdot (-\frac{d_{\mathrm{f}}^{\mathrm{3}}}{\mathrm{12}{\mu }_{\mathrm{w}}}{\rho }_{\mathrm{w}}({\nabla }_{\mathrm{T}}{P}_{\mathrm{g}}-\frac{\partial P_{\mathrm{c}}}{\partial s_{\mathrm{w}}}{\nabla }_{\mathrm{T}}{s}_{\mathrm{w}}))=d_{\mathrm{f}}{Q}_{\mathrm{f}\mathrm{w}}，（9）

CO₂： d_{\mathrm{f}}\frac{\partial }{\partial t}({\varphi }_{\mathrm{f}}{\rho }_{\mathrm{g}})+{\nabla }_{\mathrm{T}}\cdot (-\frac{d_{\mathrm{f}}^{\mathrm{3}}}{\mathrm{12}{\mu }_{\mathrm{g}}}{\rho }_{\mathrm{g}}{\nabla }_{\mathrm{T}}{P}_{\mathrm{g}})=d_{\mathrm{f}}{Q}_{\mathrm{f}\mathrm{g}}，（10）

式中：d_f为裂隙宽度，m；φ_f是裂隙的孔隙度，无量纲；K_f是裂隙的渗透率，m²；\nabla_T是指沿裂隙切线方向求导；Q_fw和Q_fg分别指咸水和CO₂的质量源项，kg/（m³·s）.

1.4　裂隙盖层体的变形控制方程

线弹性体的裂隙盖层变形控制方程可以表示为：

G{u}_{i,jj}+\frac{G}{\mathrm{1}-\mathrm{2}\nu }{u}_{k,kj}=K{\epsilon }_{\mathrm{s}\mathrm{1},i}-\alpha p_i-f_i,（11）

式中：u_i 是i方向上的位移，m；G是剪切模量，GPa；ν是泊松比，无量纲.α=1⁃K_c/K_s为比奥系数，无量纲.ε_s1是吸附引起的体积应变，无量纲.上式中的体积力有3个组成部分：重力（f_i ）、孔隙流体流动引起的拖曳力（αp_i ）和吸附引起的体积力（Kε_s1，i ）.此处盖层岩体的变形以及修正的孔隙度模型可参考以往的研究（Raduha et al.，2016），不再详细赘述.

2 研究算例分析

2.1　概念模型及初始条件

以随机离散裂隙中两相流渗流为基础，建立CO₂封存中盖层的密封性评价模型.根据Rutqvist and Tsang（2002）的地质参数，盖层埋深在1 300 m附近，盖层厚度为100 m，封存的CO₂主要以超临界状态存在，因此不需要考虑CO₂相态的变化.假定在超过1 000年的长期封存中，CO₂在盖层中的有效渗透深度不会达到100 m.为集中研究CO₂在随机离散裂隙网络中驱替咸水的渗流机理，本文的模型尺寸简化为20 m×20 m，如图2所示.

模拟区域的温度为315.5 K，初始液相饱和度为0.9.除了顶面以外的所有边界，法向位移都固定为零.下部边界的入口条件设定为CO₂压力注入，注入压力是以阶跃的形式从初始气相压力18 MPa增加到25 MPa.左右两侧边界均为零通量边界.

2.2　生成随机离散裂隙网络

目前广泛使用的裂隙岩体渗流模型通常分为3大类：等效连续介质模型（ECM）、离散裂隙网络模型（DFN）和离散裂隙基质模型（DFM）（Berre et al.， 2019）.其中，等效连续介质模型常用裂隙岩体的非均质各向异性的渗透率代替裂隙.这种理论模型虽然能体现裂隙和基质系统在渗流机理上的差异，被很好地应用于连续介质力学的有限元计算，但是并不能体现出复杂的裂隙网络影响；离散裂隙网络模型只考虑裂隙中的渗流，可以有效表征裂隙岩体的渗流特征，但忽略了基质的透水性，不能表现非连通裂隙对于岩体渗流的影响；离散裂隙基质模型可以显式表达部分支配性的裂隙，并且裂隙和基质之间可以相互交换流体，这样能更好地描述区域渗流问题.因此，本文将引入离散裂隙基质模型，进一步探究CO₂封存时裂隙分布及走向对页岩裂隙盖层中两相渗流的影响规律.

页岩体有众多裂隙网络，结构复杂.统计发现，复杂的随机裂隙网络仍然服从一些空间几何分布规律，如裂隙的中心点服从均匀分布、裂隙迹长服从负指数分布或者对数正态分布.裂隙由中心点、迹长和角度确定，角度的定义为由X轴逆时针旋转至迹线的度数.因此，基于Monte Carlo方法，利用MATLAB软件随机生成一组二维离散裂隙网络，读取每条裂隙的中心点坐标和迹长.

本文欲研究不同宏观方向上裂隙中的流体流动对于盖层密封性的影响，所以在保证中心点坐标和迹长相同的基础上改变裂隙倾角大小，利用CAD生成COMSOL可识别的dxf文件后进行模拟计算.根据裂隙倾角的不同，裂隙可以分为水平裂隙（倾角≤15°）、低角度裂隙（15°<倾角≤45°）、高角度裂隙（45°<倾角≤75°）和垂直裂隙（75°<倾角≤90°）（白雪峰，2009；He et al.， 2022；曾联波等，2023）.参考五峰组‒龙马溪组页岩纵向发育不同的岩相：五峰组⁃S1l1亚段的1~3亚层和4亚层下部是以高角度或垂直裂隙为主；4亚层中上部是以低角度、水平裂隙发育，少量垂直裂隙为主的裂隙（He et al.，2022）.其中，为体现裂隙网络整体的方向性，每一组裂隙的两种走向均对称，共生成6组离散裂隙网络，如图3所示.A组裂隙代表的是以水平裂隙为主的层间离层裂隙（0°~180°）；B组裂隙代表的是低角度裂隙（30°~150°）；C组裂隙代表的是低角度裂隙和高角度裂隙的临界（45°~135°）；D组裂隙代表的是高角度裂隙（60°~120°）；E组裂隙代表的是以垂直裂隙为主的纵向破断裂隙 （90°~90°）；F组裂隙为正交裂隙（90°~180°）.

3 结果与讨论

3.1　裂隙走向对液相饱和度分布的影响

裂隙网络会直接影响CO₂驱替咸水时的驱替前端渗透形态和速度，从而影响盖层密封CO₂的效率.随着储层中的CO₂不断累积到盖层底部，因为页岩基质渗透率远低于裂隙的渗透率，气‒液两相流的驱替首先会发生在裂隙网络中.观察不同随机裂隙网络模型在封存时间为1年、10年和50年时的液相饱和度分布，可以发现液相饱和度的分布能直接反映渗透前端.

在第1年的时候，如图4所示，6组裂隙网络下的CO₂渗透前端已经开始出现较为明显的差别.A组水平裂隙对垂直方向上的CO₂驱替咸水过程有抑制作用，所以该模型驱替前端的移动距离最小.B、C和D组裂隙的驱替前端形态呈“V”状，并且驱替前端随着裂隙角度的增大而增长（即B组30°~150°最小，C组45°~135°次之，D组60°~120°最大）.E组的纵向裂隙虽然有利于CO₂在垂直于盖层方向上的渗透，但是单一的裂隙走向并不能形成裂隙网络，液相饱和度的变化基本都只集中于裂隙处.F组裂隙的驱替前端形态呈“十”字状，驱替前端的移动速度也会受到水平裂隙的阻碍，但其区域影响范围大于单一方向的裂隙网络.另外，可以观察到部分区域中的液相饱和度已经小于初始液相饱和度0.9了，这表明液相已被有效驱替.

第10年时，如图5所示，每一组裂隙基质模型中靠近CO₂注入侧的区域液相饱和度均出现明显下降，其中相交裂隙基质模型相较于单一裂隙基质模型饱和度下降得更快，驱替效果更好.A组裂隙的驱替前端虽然仍然落后于其他5组裂隙，但驱替效果显著——液相饱和度已接近残余饱和度.B、C和D组裂隙基质模型中的渗透前端依次递增，D组裂隙网络的影响范围会更大，驱替效果要优于其他两组.E组裂隙基质模型中裂隙密度大的区域驱替速度大于裂隙密度小的区域，这导致液相饱和度的分布图出现了一定的层次.F组裂隙基质模型由于中间部分水平裂隙的阻碍作用，两侧的驱替前端远超中间部分.总而言之，相交裂隙基质模型的驱替效果更佳，这是因为连通裂隙网络中存在有效的流体流动.

第50年的时候，如图6所示，A组裂隙基质模型中CO₂驱替咸水的渗透前端虽然远远落后其他几组，但其整体的渗透前端仍较为一致，其余5组裂隙基质模型均受到裂隙分布的影响——左右两侧的渗透前端有较大的差距，尤其是在离散裂隙更为密集的地方，饱和度会出现较为明显的改变——裂隙处的饱和度会明显低于其影响区域内的饱和度，没有裂隙甚至离散裂隙连通性较差区域的饱和度变化会稍微滞后于裂隙处的饱和度变化.除此之外，B、C和D组裂隙模型中均有部分区域出现回流现象，这是因为随机生成的裂隙网络具有非均匀性，导致区域间的连通性能不同——存在优势通道.E组和F组裂隙基质模型回流现象不明显，这是因为E组裂隙网络未能形成连通网络，而F组裂隙网络中有水平裂隙的阻碍作用.因此，裂隙网络的走向会直接影响两相流分布情况.

3.2　裂隙走向对CO₂渗透深度的影响

在以往的研究中，通过液相饱和度分布来定义CO₂的渗透深度——当液相饱和度沿Y轴不断增大，直到第一次出现下降时，这个转折点所对应的Y轴坐标被定义为CO₂渗透深度（Wang and Wang， 2018）.但是，考虑到显式裂隙分布的不均匀性会导致纵向渗透速度不同，如果只是利用某一垂线则无法判断非均匀的渗透深度，因此本文将渗透深度定义为：若某一点的液相饱和度开始持续低于域内初始饱和度值，则将该点所对应的Y轴坐标视为该时刻的CO₂的渗透深度.

本文选取B组裂隙网络走向（30°~150°）的裂隙基质模型为研究对象，探究在10 m处的液相饱和度分布情况以判断渗透深度的变化.从图7可以发现，不同位置的液相饱和度存在一定差异，随着时间推移该位置上的液相饱和度会随之减小.第20年时，该高度上所有点的液相饱和度均小于初始值0.9，所以在第20年时CO₂的渗透深度已经超过10 m了.

为探究不同裂隙倾角的裂隙基质模型对于CO₂驱替时间的影响，选取相同的渗透深度10 m，利用数值模拟软件计算并结合10 m处的液相饱和度分布曲线，得出不同倾角裂隙基质模型所需要的驱替时间及驱替速度，具体结果如图8所示.

可以发现，裂隙基质模型中的裂隙倾角越大驱替时间越短，即越有利于纵向上的渗透.随着裂隙角度的增加，驱替时间降低了12.59倍，驱替速度增加1个数量级.其中，0°~0°的裂隙基质模型所需要的驱替时间最长，这是因为单一的水平裂隙不仅不能形成有效的裂隙网络，而且对咸水的驱替过程起到阻碍作用.另外，60°~120°、90°~90°和90°~180°裂隙基质模型所需要的驱替时间相同，这表明当裂隙角度增加到60°后，对于驱替时间的影响逐渐变得不明显.此时，这3种裂隙基质模型的10 m处液相饱和度分布如图9所示.

在达到相同的渗透深度10 m时，不同裂隙基质模型在该处有不同的液相饱和度分布.总而言之，以垂直裂隙为主的纵向破断裂隙（90°~90°）基质模型的饱和度最大，相交裂隙（60°~120°）次之，正交裂隙（90°~180°）最小.结合模型整体的液相饱和度分布云图，整体饱和度越大表明驱替到该深度时液相残余饱和度越大，其驱替效果越差.因此在这3种裂隙基质模型中，正交裂隙（90°~180°）的裂隙基质模型的整体驱替效果最佳，其次是相交裂隙 （60°~120°），最差的是以垂直裂隙为主的纵向破断裂隙（90°~90°）.这是因为90°~180°的裂隙基质模型在形成了有效的连通裂隙网络的基础上还有纵向裂隙促进咸水的渗透，60°~120°的裂隙基质模型虽然也形成了有效的裂隙网络但其裂隙所影响的纵向深度小于90°~180°裂隙模型所影响的纵向深度，而90°~90°的裂隙基质模型未能形成连通裂隙网络.

当裂隙网络的走向平行于盖层厚度的方向时，CO₂的优势流动方向与突破盖层的方向一致，渗透深度快速增加；且CO₂沿优势裂隙侵入后，孔隙流体压力变大，有效应力的减小导致裂隙开度变大.若是垂直于盖层厚度方向的裂隙网络，其优势渗流方向则与突破盖层方向垂直，虽然在孔隙流体压力的影响下也会导致裂隙开度变大，但此时CO₂压力在横向流动中被消耗，反而减小了突破盖层的风险.

3.3　裂隙走向对裂隙盖层密封效率的影响

当大量的超临界CO₂注入到深部储层时，毛细压力封存成为短时间内最有效的封存机制.但是少量的CO₂可以穿透盖层，这种CO₂的逃逸可能带来许多风险，例如地下水酸化、重金属迁移和矿物溶解（Bacon et al.， 2016；付葆春和许媛，2022）.由于很难直接计算出盖层中逸出的CO₂量，如果CO₂在渗透过程中未从数值模型顶部边界突破，评估盖层密封效率的另一种有效方法可以是测算CO₂在盖层中的渗透量.

在不考虑盖层厚度的情况下，本文直接使用以下公式预估渗透到裂隙盖层中的CO₂量：

M_{\mathrm{C}{\mathrm{O}}_{\mathrm{2}}}=C_{\mathrm{c}}\iint \varphi {\rho }_{\mathrm{g}}{S}_{\mathrm{g}}\mathrm{d}x\mathrm{d}y，（12）

式中：C_c是盖层封存效率系数，由页岩的各向异性和CO₂的流动性等多种因素共同决定.根据美国地质勘探局在区域尺度上对192个封存评估单元的长期稳态CO₂资源的评估，本文参考其取值，假定C_c的取值为0.6（Bachu， 2015）.

通过式（12）对单位面积的盖层进行积分，得到离散裂隙基质模型中CO₂的渗透量随时间变化的曲线.如图10所示，每种裂隙基质模型的CO₂渗透量均呈增长的趋势，这表示此时尚未有CO₂从盖层中逃逸.在此情况下，将这6种裂隙基质模型按照CO₂渗透量从大到小依次排序为：F组裂隙网络（90°~180°）>C组裂隙网络（45°~135°）>D组裂隙网络（60°~120°）>E组裂隙网络（90°~90°）>B组裂隙网络（30°~150°）>A组裂隙网络（0°~180°）.这是因为裂隙盖层中由于渗透率分布不均匀，大部分CO₂会先运移到预先存在的裂隙处.当裂隙网络能够有效连通时就能够形成高渗透率区域，这导致更多的CO₂会优先流入并扩散.如果仅从单个裂隙粗略估计，A组和B组裂隙网络在水平方向上的渗透能力要大于竖直方向上的渗透能力，尤其是A组裂隙网络只能促进水平方向上的渗透；D组和E组裂隙网络在竖直方向上的渗透能力要大于水平方向上的渗透能力，而E组裂隙网络只能促进垂直方向上的渗透；另外，C组和F组裂隙网络水平方向上的渗透能力与竖直方向上的渗透能力相同，两个方向的渗透率比值均为1∶1.总而言之，若CO₂在盖层中的渗透量越大，盖层的密封效率越高，但越容易出现CO₂泄漏现象.

因此，若在盖层中形成了较好的随机离散裂隙网络，CO₂的渗透量也会随之增加.在假定CO₂驱替咸水的前端不会突破盖层上部的情况下，更多的CO₂会渗透到裂隙的毛细管内.随着时间的推移，裂隙中的CO₂还会通过缓慢的扩散过程进入到页岩基质中，吸附在其孔隙表面.

4 结论

（1）基于随机离散裂隙网络的盖层模型发现：裂隙网络的走向会直接影响了两相流分布，进而影响液相饱和度的赋存情况.以水平裂隙为主的层间离层裂隙模型（0°~180°）的驱替前端远远落后于其他五组裂隙；相交裂隙模型（30°~150°、 45°~135°和60°~120°）的驱替前端形态呈“V”状，并且驱替前端随着裂隙角度的增大而增长；以垂直裂隙为主的纵向破断裂隙模型（90°~90°）中液相饱和度的变化基本都只集中于裂隙处；正交裂隙模型（90°~180°）的驱替前端形态呈“十”字状，驱替前端的移动速度会受到水平裂隙的阻碍.

（2）考虑到显式裂隙分布的不均匀性，将渗透深度的定义进一步改善后发现：裂隙基质模型中的随着裂隙网络角度（0°~60°）的增大驱替时间降低了12.59倍，但当裂隙倾角增大到60°后不再有明显变化.此时，正交裂隙模型的整体驱替效果最佳，其次是相交裂隙，最差的是以垂直裂隙为主的纵向破断裂隙.

（3）假定CO₂驱替咸水的前端不会突破盖层上部的情况下，有更多的CO₂会渗透到正交裂隙F（90°~180°）的毛细管内.随着裂隙网络影响范围的增大，模型中CO₂渗透量越大，盖层的密封效率越好.在更长的时间内，裂隙中的CO₂还会通过缓慢的扩散过程进入到页岩基质中，吸附在其孔隙表面.

参考文献

[1] Adu⁃Gyamfi, B., Ampomah, W., Tu, J. W., et al., 2022. Assessment of Chemo⁃Mechanical Impacts of CO₂ Sequestration on the Caprock Formation in Farnsworth Oil Field, Texas. Scientific Reports, 12(1): 1-23. https://doi.org/10.1038/s41598⁃022⁃16990⁃x

[2] Arif, M., Al⁃Yaseri, A. Z., Barifcani, A., et al., 2016. Impact of Pressure and Temperature on CO₂–Brine– Mica Contact Angles and CO₂–Brine Interfacial Tension: Implications for Carbon Geo⁃Sequestration. Journal of Colloid and Interface Science, 462: 208–215. https://doi.org/10.1016/j.jcis.2015.09.076

[3] Bachu, S., 2015. Review of CO₂ Storage Efficiency in Deep Saline Aquifers. International Journal of Greenhouse Gas Control, 40: 188–202. https://doi.org/10.1016/j.ijggc.2015.01.007

[4] Bacon, D. H., Qafoku, N. P., Dai, Z. X., et al., 2016. Modeling the Impact of Carbon Dioxide Leakage into an Unconfined, Oxidizing Carbonate Aquifer. International Journal of Greenhouse Gas Control, 44: 290-299. https://doi.org/10.1016/j.ijggc.2015.04.008

[5] Bai, X.F., 2009. Fracture Characteristics of Volcanic Reservoir in Yingcheng Formation, Songliao Basin. Global Geology, 28(3): 318-325 (in Chinese with English abstract).

[6] Berre, I., Doster, F., Keilegavlen, E., 2019. Flow in Fractured Porous Media: A Review of Conceptual Models and Discretization Approaches. Transport in Porous Media, 130(1): 215-236. https://doi.org/10.1007/s11242⁃018⁃1171⁃6

[7] Chen, B.W., Wang, R., Li, Q., et al., 2023. Status and Advances of Research on Caprock Sealing Properties of CO₂ Geological Storage. Geological Journal of China Universities, 29(1): 85-99 (in Chinese with English abstract).

[8] Chen, F., Yang, B.D., Ma, J., et al., 2016.Effects of CO₂ Geological Storage Leakage on Soil Microbial Community. Journal of China University of Mining & Technology, 45(6): 1285-1293, 1299 (in Chinese with English abstract).

[9] Currie, N.D., 2017.Characterizing Injectite Related Fracture Networks and Their Impact on Seal Integrity, Rainbow Rocks Injectite Complex,Grand County,Utah (Dissertation). New Mexico Institute of Mining and Technology, New Mexico.

[10] Dai, L., You, L., Zhao, Z.J., et al., 2019. Petrological Characteristics and Sealing Analysis of Mesozoic Caprock in Dongfang⁃Ledong Area, Yinggehai Basin. Earth Science, 44(8): 2677-2685 (in Chinese with English abstract).

[11] Deng, Y.R., Wang, Y.H., Zhao, Y.J., et al., 2023. Carbon Dioxide Storage in China: Current Status, Main Challenges, and Future Outlooks. Earth Science Frontiers, 30(4): 429-439 (in Chinese with English abstract).

[12] Fu,B.C., Xu,Y., 2022.Preliminary Study on Environmental Monitoring and Assessment System for CO₂ Geological Storage Project. Shanxi Chemical Industry, 42(8):139-140,143 (in Chinese with English abstract).

[13] Guo, B., 2020.Study on the Sealing Ability Evolution of Caprock System in the Process of CO₂ Geological Storage: A Case Study of Ziniquanzi Formation Mudstone Caprock in Xinjiang Junggar Basin (Dissertation). Jilin University, Changchun (in Chinese with English abstract).

[14] He, S., Qin, Q. R., Qin, Z. J., et al., 2022. Natural Fracture Development Characteristics and Their Relationship with Gas Contents─A Case Study of Wufeng– Longmaxi Formation in Luzhou Area, Southern Sichuan Basin, China. ACS Omega, 7(38): 34066-34079. https://doi.org/10.1021/acsomega.2c03318

[15] Hyman, J. D., Jiménez⁃Martínez, J., Gable, C. W., et al., 2020. Characterizing the Impact of Fractured Caprock Heterogeneity on Supercritical CO₂ Injection. Transport in Porous Media, 131(3): 935-955. https://doi.org/10.1007/s11242⁃019⁃01372⁃1

[16] Jiang,S.,Zhang,K.,Du,F.S.,et al.,2023. Progress and Prospects of CO₂ Storage and Enhanced Oil, Gas, and Geothermal Recovery. Earth Science, 48(7): 2733-2749 (in Chinese with English abstract).

[17] Olabode, A., Radonjic, M., 2017. Fracture Conductivity Modelling in Experimental Shale Rock Interactions with Aqueous CO₂. Energy Procedia, 114: 4494-4507. https://doi.org/10.1016/j.egypro.2017.03.1610

[18] Raduha, S., Butler, D., Mozley, P. S., et al., 2016. Potential Seal Bypass and Caprock Storage Produced by Deformation⁃Band⁃to⁃Opening⁃Mode⁃Fracture Transition at the Reservoir/Caprock Interface. Geofluids, 16(4): 752-768. https://doi.org/10.1111/gfl.12177

[19] Rutqvist, J., Tsang, C. F., 2002. A Study of Caprock Hydromechanical Changes Associated with CO₂⁃Injection into a Brine Formation. Environmental Geology, 42(2): 296-305. https://doi.org/10.1007/s00254⁃001⁃0499⁃2

[20] Shang, X., Wang, J., Wang, H., et al., 2022. Combined Effects of CO₂ Adsorption⁃Induced Swelling and Dehydration⁃Induced Shrinkage on Caprock Sealing Efficiency. International Journal of Environmental Research and Public Health, 19(21): 14574. https://doi.org/10.3390/ijerph192114574

[21] Song, J., Zhang, D. X., 2013. Comprehensive Review of Caprock⁃Sealing Mechanisms for Geologic Carbon Sequestration. Environmental Science & Technology, 47(1): 9-22. https://doi.org/10.1021/es301610p

[22] Verdon, J. P., Kendall, J. M., Stork, A. L., et al., 2013. Comparison of Geomechanical Deformation Induced by Megatonne⁃Scale CO₂ Storage at Sleipner, Weyburn, and in Salah. Proceedings of the National Academy of Sciences of the United States of America, 110(30): E2762-E2771. https://doi.org/10.1073/pnas.1302156110

[23] Wang, J. G., Peng, Y., 2014. Numerical Modeling for the Combined Effects of Two⁃Phase Flow, Deformation, Gas Diffusion and CO₂ Sorption on Caprock Sealing Efficiency. Journal of Geochemical Exploration, 144: 154–167. https://doi.org/10.1016/j.gexplo.2013.12.011

[24] Wang, J. G., Wang, H. M., 2018. Sealing Efficiency Analysis for Shallow⁃Layer Caprocks in CO₂ Geological Storage. Environmental Earth Sciences, 77(21): 1-14. https://doi.org/10.1007/s12665⁃018⁃7924⁃2

[25] Wang, Z.N., Yuan, J.S., Pang, B., et al., 2022. The Interpretation and Highlights on Mitigation of Climate Change in IPCC AR6 WGⅢ Report. Climate Change Research, 18(5):531-537 (in Chinese with English abstract).

[26] Wu, Z.D., 2011. Research on the Stability of Cavern in Bedded Rock Salt Based on Permeability and Time Dependent Behavior (Dissertation). China University of Mining and Technology, Beijing (in Chinese with English abstract).

[27] Yin, Q., Jing, H. W., Su, H. J., et al., 2017. CO₂ Permeability Analysis of Caprock Containing a Single Fracture Subject to Coupled Thermal⁃Hydromechanical Effects. Mathematical Problems in Engineering, 2017(1): 1-13. https://doi.org/10.1155/2017/1290748

[28] Zeng,L.B.,Ma,S.J.,Tian,H., et al.,2023. Research Progress of Natural Fractures in Organic Rich Shale. Earth Science, 48(7): 2427-2442 (in Chinese with English abstract).

[29] Zhao, Z. Y., Yao, S., Yang, S. P., et al., 2023. Under Goals of Carbon Peaking and Carbon Neutrality: Status, Problems, and Suggestions of CCUS in China. Environmental Science, 44(2):1128-1138 (in Chinese with English abstract).

[30] Zhou, X.F., Zhang, Y.S., Yan, D.T., et al., 2018. Quantitative Evaluation of Sealing Ability of Tertiary Mudstone Caprock in Lenghu Area of Qaidam Basin. Earth Science, 43(S2): 226-233 (in Chinese with English abstract).

基金资助

国家自然科学基金项目(42202286)

国家自然科学基金项目(62103198)

江苏省博士后科研资助项目(2021K350C)

河海大学自由探索专项‒新引进教师项目(B210201037)

河海大学大型仪器设备共享基金(B220370201)