二氧化碳地质封存中盖层力学完整性数值模拟研究综述

刘苗苗，孟令东，王海学，吴 桐

(1.东北石油大学，河北 秦皇岛 066004；2.黑龙江省油气藏及地下储库完整性评价重点实验室，黑龙江 大庆 163318；3.黑龙江省石油大数据与智能分析重点实验室，黑龙江 大庆 163318)

0 引 言

CO2地质封存是控制大气中CO2浓度最科学有效的手段之一[1]。将大量CO2持续注入地质封存体内，会引发一定范围内的压力传播、流体运移等现象[2-3]，对盖层力学稳定性产生重要影响[4]。要使CO2持续注入并稳定封存几百或几千年不逸散，就要保证地层压力不超过盖层突破压力，因此，需定量表征封存过程中盖层岩石的变形损伤程度，正确评价其力学完整性[5]。针对该问题，目前的研究方法主要有理论分析法[6-7]、室内实验[8]或场地浅层含水层的注气试验法[9]、数值模拟法[10-14]。相比之下，数值模拟法能有效弥补其他研究手段的缺陷，已成为当前地质工程学中定量化研究的有效方法之一。鉴于此，对CO2地质封存中盖层力学完整性分析相关问题的数值模拟研究现状进行了综合分析，以期为CO2地质封存的相关研究提供技术支持与引导。

1 CO2地质封存研究中的数值模拟方法

大量CO2的持续注入会导致岩石所处的应力状态、流体的渗流路径等发生变化[15]，对这一过程的刻画是分析岩层应力-应变规律的关键。室内实验方法只能分析该过程产生的机理，而数值模拟方法能充分考虑各个过程间复杂的耦合关系，可直观有效地对地质模型进行数学和力学分析。目前，已有许多学者从不同角度，使用不同的软件、模型及方法对CO2地质封存的相关问题进行了数值模拟研究。

1.1 多场耦合原理

研究显示，将CO2以超临界态注入咸水层封存是最为理想的方式，而CO2咸水层封存则是一个复杂的热力学-流体力学-岩石力学-化学多场相互作用的过程[16]，如图1所示。由于CO2注入过程一般只有几十年，这期间化学场对盖层力学完整性的影响作用几乎可以忽略。此外，研究显示，CO2注入引起的温度变化一般仅局限于注入井周围的几十米[17]。因此，许多盖层力学完整性数值模拟研究均忽略温度的变化，仅研究渗流-应力耦合过程。根据前人研究[18-19]，流体模块基础模型中渗透率纵横比通常取0.1，相对渗透率多采用VGM模型，而毛细压力模型采用VG函数。对于力学模型，一般的数值模拟研究都假定CO2注入引起的垂向总应力不变，而地层侧向变形受限，只在垂直方向上变形。

图1 CO2咸水层封存多场耦合原理

1.2 力学本构模型

用于CO2地质封存力学模拟的岩体本构模型主要有弹性、塑性和非连续介质模型3类[20]。

(1) 弹性本构模型。不能描述复杂应力下的塑性变形状态，但其简单快捷，计算结果收敛，易于进行高耦合程度的多场耦合计算，因此，在CO2地质封存力学模拟中仍占较大比例[21]。

(2) 塑性本构模型。主要包括弹塑性和黏塑性本构模型。复杂应力条件下，储层和盖层岩体的力学特征并不遵循弹性介质的应力-应变曲线，此时选择弹塑性本构模型能更好地刻画其变形屈服过程，尤其当模拟储层和盖层的长期蠕变特征时，可使用黏塑性本构模型[22]。

(3) 非连续介质模型。针对CO2的注入可扰动初始应力场，导致新裂隙产生的问题，主要利用非连续介质模型进行储层及盖层的岩体力学模拟[23]。非连续介质模型能更好地刻画CO2封存系统中储层及盖层的裂隙特征，因此，使用其进行裂隙启闭以及扩展过程的模拟是未来数值模拟领域的发展趋势。

当前的数值模拟方法主要按照以下规则分析盖层的力学完整性：对于弹性本构模型，盖层单元应力状态满足剪切或张拉破坏准则[24]；对于弹塑性本构模型，盖层底部单元出现塑性屈服[25]；对于黏塑性本构模型，盖层单元的塑性应变达到设定的某个临界值[22]；对于非连续介质模型，盖层单元节点断开或已有裂隙发生张开或滑移[26]。

1.3 耦合分析方法

针对CO2地质封存中的多场耦合问题，存在多种数值分析方法，根据求解过程耦合程度的不同，可大致归结为全耦合、弱耦合与单向耦合3类。

(1) 全耦合。需要同时求解热力学、渗流力学和岩石力学这3类特征差异较大的方程，所有参数在每个迭代步计算完成后同时更新代入刚度矩阵，用于下一时间步的迭代计算。该方法结果准确，但需花费大量时间，因此，适用于尺度较小、相对简单的几何模型，例如仅用于热-水动力场之间的耦合，且一般采用弹性和塑性本构模型[27]。

(2) 弱耦合。将偏微分方程拆解为流体模块和力学模块，每个模块在一个时间步内先按照一定的顺序序列独自完成计算，再将结果代入对方下一时间步的计算中修正输入参数，如此反复迭代执行直至结束。该方法计算结果精度有所降低，但算法实现相对容易，速度快，收敛性较好，已成为当前CO2地质封存研究中应用最广泛的耦合分析方法。

(3) 单向耦合。各个模块独立运行，计算效率最高，收敛性最好，但力学与流体2个独立模块的计算结果不反馈给对方。因此，该方法多用于只关注力学而不注重渗流过程的模拟，且往往用于分析非连续介质模型[28]。

1.4 数值模拟分析软件

能够用于CO2地质封存流固耦合过程分析的数值模拟软件有很多(表1)。渗流场分析可采用TOUGH2、ECLIPSE和CMG-GEM等[29]，应力场分析可采用FLAC3D、Abaqus和Petrel等。针对CO2地质封存流固耦合数值模拟研究，使用最多且效果较好的是FLAC3D和TOUGH2。FLAC3D可对CO2注入后盖层岩石应力-应变的整个动态过程进行追踪与反演，实现不同时空尺度的模拟及定量化研究，具有极强的可塑性与应用灵活性，同时可节约成本，兼顾各个因素的共同影响。而TOUGH2则拥有通用的多相流动方程及热对流传导数学模型，能够对多相流、多组分、非等温流体流动中的热和水动力耦合过程进行分析，且模拟结果准确可靠。由于TOUGH2不能刻画由水动力变化而导致的岩体变形过程，而力学过程的模拟在CO2地质封存问题研究中至关重要。因此，一些学者设计了TOUGH2-FLAC3D接口程序，对CO2地质封存中的相关问题进行了模拟研究。此外，为了对比分析不同数值模拟器的模拟结果，学者们分别采用了TOUGH2-FLAC3D、TOUGH2/ECO2N、TOUGHREACT等多种模拟器进行了研究[30]。目前，基于TOUGH2-FLAC3D的耦合分析已被成功用于模拟CO2注入咸水层后压力的抬升、盖层的力学性质变化等问题。

表1 数值模拟软件对比分析

2 CO2咸水层封存中盖层力学稳定性数值模拟研究现状

2.1 用于岩体力学模拟分析的数值计算方法

CO2地质封存中储层的力学响应依赖于温度、流体压力及应力变化，而盖层及远离注入点的区域力学响应则主要依赖于应力变化。注入初期，流体压力相对较小，岩体力学响应是弹性可逆的。随着注入压力增大，岩体可能出现不可逆的力学损伤。当压力足够大时，还有可能导致剪切或拉张破坏。为模拟这一物理过程，尤其是力学方面的变化情况，需选择合适的数值计算方法。针对地质力学问题，常用的数值计算方法有3类：适用于连续介质的(如有限单元法FEM、有限差分法FDM)、适用于非连续介质的(如离散元法DEM、非连续变形分析DDA)，以及用于连续-非连续介质的数值方法(如混合FEM/BEM)。对于这些数值方法，不能武断地说某种方法优越于另一种。在地质力学建模中，数值方法的选择需考虑较多因素，如研究目的、模拟的裂缝几何体尺度等。研究显示，基于连续介质的数值方法在CO2封存的力学建模分析中占主导地位[31]。相关研究中常将数值计算方法组合起来或链接到相应的模拟软件，并使用弹性或弹塑性本构模型进行耦合分析，刻画岩体地质力学行为。以连续介质的数值计算方法为例，THM耦合分析时序方案如图2所示。

图2 THM耦合分析时序方案示意图

2.2 储、盖层应力变化数值模拟

针对CO2咸水层封存中储、盖层的力学稳定性问题，一些学者从大尺度角度，使用FLAC3D、TOUGH等模拟了实际封存工程中储层流体压力及上覆盖层的应力变化情况。

(1) 随着大量CO2的注入，储、盖层孔隙流体压力均逐渐增加，压力响应对空间离散更敏感。CO2开始注入后，压力在注入点周围迅速积聚，形成以注入点为中心、向四周逐渐减弱传播的环型高压集聚带。压力影响的波及速度很快，且早期压力上升范围快于后期[32-33]。

(2) 储层流体压力随注入时间的增加而增大，随径向距离的增大而减小；水平和垂向的总应力都增加，且前者增量大于后者；水平和垂向有效应力都减小。CO2注入阶段，储层流体压力初始阶段变化较快，随后逐渐稳定，且注入井附近压降较小，当径向距离增大到一定程度时，压降突然增大直到降至地层初始压力[34]。随着CO2持续注入，多个环形压力集聚带逐渐相连、干扰重叠，形成一定面积的压力抬升区[35]。

(3) 储层压力抬升很快传递至盖层，盖层渗透率对压力传播有重要影响。随着注入量的增加，注入点附近压力明显大于远处区域压力；整个注入阶段，储、盖层交界处有效应力变化及岩石变形最大，最易发生岩石断裂[36]。垂向压力的变化很大程度上依赖于盖层的参数，如渗透系数及压缩性[37]。从时间角度而言，低渗透盖层的压力积聚速度快于高渗透盖层，注入井附近压力的最大值随着盖层渗透率的增加而减小[38]。盖层渗透率的提高虽然能够避免储层内部出现过高压力，但却增加了盖层泄漏风险[39]。

2.3 储、盖层压力增量敏感性数值模拟

对CO2注入咸水层后储、盖层压力增量敏感性进行了数值模拟分析，结果如下。

(1) 以注入井处的压力为响应变量时，孔隙度的敏感性最高[40]。

(2) 储层渗透率纵横比越小，注入井附近垂向压力的增量越大；高渗透率储层可以使压力波动传播得更快，且在远离注入井的区域观测到的压力峰值更大[11]。

(3) 盖层渗透率及压缩系数的变化对储层压力增量有较大影响。当盖层渗透系数变大时，垂向压力变化较大；当盖层厚度较小且渗透率相对较高时，局部封盖效果不佳，渗透率和渗透率纵横比的下降会导致压力上升幅度增加[18]。

(4) 孔渗非均质的岩层中孔隙压力的增量要高于均质岩层中孔隙压力的增量，并且随时间增加，二者之间的差异也逐渐增大；孔隙度和渗透率的非均质性会导致孔隙流体压力大幅增加，提高CO2泄漏的风险[41]。

(5) 温度低不利于压力消散，当注入的CO2温度与目标地层温度差超过一定值时，冷却效应可能导致注入区域出现拉张裂缝[42]。

除以上因素外，储盖层的几何尺寸、岩体的本构模型及参数选择、边界范围设定等也会对应力增量产生影响[43-44]。例如：边界范围设定引起的压力增量随范围的增大而逐渐减小；定压力边界范围设置过小会低估储层压力增量，设定过大又会增加计算量。

2.4 盖层水力破裂风险评价数值模拟

当前，CO2地质封存中盖层力学完整性评价主要包含2个方面：一是注入压力扰动引起的储集层局部高压引发的储、盖层张拉破坏风险；二是复杂的地质环境及岩石物性变化导致局部应力集中，由此引发剪切破坏风险。

(1) 盖层张拉破坏风险评价。开展CO2咸水层封存工程之前，需要测试封存场所的地应力情况，以准确控制注入压力的上限，保障盖层的力学完整性。通常情况下，通过开展水力压裂等场地地应力测试和室内岩心实验，来准确测试储、盖层的最小水平主应力，以此评估现场注入压力是否会引发储、盖层的张拉破坏。

(2) 盖层剪切破坏风险评价。CO2咸水层封存中，盖层剪切破坏风险评价是在岩石力学实验的基础上进行的，一般采用数值模拟方法，通过研究封存系统应力场及其动态变化特征，准确模拟初始地应力场及CO2注入产生的压力扰动，依据剪切破坏准则计算破坏指标，从而评估盖层的剪切破坏风险。

(3) 基于地层初始应力状态分析的盖层水力破裂形式数值模拟。CO2封存中，初始的地应力状态是影响岩石破坏形式的主要因素。当采用弹性本构模型时，无论地层初始应力处于拉伸还是压缩应力状态，发生剪切破坏的可能性都高于张拉破坏的可能性；在压缩应力状态下，剪切破坏面不会穿透盖层，不影响其封存性；在拉伸应力状态下，剪切破坏面通常会贯穿整个盖层区域，破坏面的角度(大约为60 °)比压缩应力状态下(大约为30 °)更大[45]；且拉应力状态下，盖层底部发生拉张破坏的可能性更大。雷宏武等[11]的研究显示，CO2的注入引起流体压力急剧增加，地层有效应力减小，可能发生剪切破坏的方向在注入过程中变化不大，剪切破坏最有可能出现在最大速率注入点上部的盖层，尤其是在紧靠第2注入层上部的盖层，其次为靠近地表的位置，这与Rutqvist等[45]的分析结果相似。

3 当前研究中存在的主要问题

从总体现状来看，流固耦合下盖层力学稳定性研究取得了一定成果，但有关地质封存中盖层水力破裂机理的定量表征及力学完整性评价方面的数值模拟研究还存在以下问题。

(1) 相关研究所建立的地质模型大多为概化模型，受软、硬件条件的限制，模型中对盖层地质特征的描述精度有待进一步提高。由于地震和测井精度一般大于50 m，而盖层厚度通常小于50 m，因此，很难通过地震和测井资料对盖层进行精确描述。此外，通过室内实验和露头得到的地质资料具有分散性、扰动性，不能完整精确表征盖层的地质特点。由于缺少野外实际资料，目前所建数值模型对地层描述通常进行了不同程度的简化。例如，将实际地层概化为水平地层，非均质介质概化为均质介质，且大部分模型忽略了断层的封闭-开放性的影响、毛细压力函数和相对渗透率模型取值变化的影响等。这些影响因素都给建模过程中对盖层特性的精确描述带来了很多不确定性，直接影响数值模拟结果的准确性。

(2) 由于缺乏充足的实验数据，目前还没有较权威的能够准确表征CO2地质封存中盖层岩石应力应变的本构关系。当前的数值模拟研究多采用传统的多相流体在多孔介质中的渗流理论与控制方程，由于相对渗透率模型中参数的取值难度较大，研究者只能根据现场监测数据及为数不多的室内实验数据推断现有模型中各参数的取值范围，来近似表征岩石的应力应变及多相渗流特性。虽然已有大量学者对超临界CO2的性质进行了一定研究，但对超临界CO2流动所具有的特殊性质方面的研究较少。且模型中的方程或定律都只适用于特定的范围及条件，当某场地的实际地质条件与模型的适用范围相差较大时，已有模型将不再适用。

(3) 就总体现状而言，目前仍以搭接不同的模拟器来研究CO2地质封存中的热流固耦合过程为主，缺乏单一的THM耦合模拟器。当前，已有的力学问题耦合效应数值模拟研究中，有的单从流场角度分析储层封存的可行性，研究了裂缝对于渗流场的影响，而很少同时关注其对应力场的影响；有的单从力学角度考虑盖层是否会失稳，关注了CO2注入可能诱发的裂缝产生问题，而较少关注裂缝产生后对渗流场的影响。相关研究缺乏多种机理相互影响和作用下的系统性研究，难以精确描述CO2地质封存中盖层发生水力破裂的影响因素、裂缝扩展及CO2泄漏的真实过程，导致模拟结果与实际封存效果常常存在较大差异。

(4) 缺乏考虑地层参数不确定性的盖层力学完整性评价模型。CO2地质封存系统中，地质参数、储层连续性、实际工程情况等均存在较大不确定性[46]。通常情况下，储、盖层的物理性质、初始地应力场以及岩石力学参数是依据测井数据通过数学模型计算获得，然而由于地质环境的复杂性以及人为因素导致的测井数据资料的不完备性，使得地层参数存在大量的不确定性。因此，模型中的参数多为特定条件下的测定结果，并不能真实反映自然条件下的复杂特征。而目前的研究中缺少地层参数不确定性对储、盖层力学响应影响的研究，关键参数取值的不确定性会直接导致数值模拟的预测结果与实际情况出现偏差。

4 未来研究趋势

(1) 精确的地质建模对于模拟结果的正确性至关重要，为了得到更接近于实际情况的模拟结果，精确表征盖层地质特性成为了盖层地质建模中最基础和最重要的问题，在后续研究中有待进一步探索，以提高对实际模型中地质特征的描述精度，为CO2地质封存提供更为精确的参考。

(2) 亟需开展能够精确表征CO2咸水层封存中盖层岩石应力应变的本构关系的研究，以便更加准确地描述超临界CO2注入之后岩石物性参数的变化特征，以及注入压力对岩石应力-应变损伤特征的影响，为盖层水力破裂机理定量表征及密封性失效风险评价奠定良好的基础。

(3) 有待设计单一的能够精确表征CO2地质封存特征的THM耦合模拟器，开展多尺度条件下盖层力学完整性数值模拟研究，定量表征水力破裂机理，基于原地应力场特征和破裂机理，获取极限注入压力等参数，有效预测断裂和裂缝启闭性，为水力破裂风险的有效评估以及盖层力学完整性的精细描述与定量评价提供准确的数据支撑。

(4) 针对地质参数的随机性，后续工作中开展敏感参数的不确定性分析与离散化研究是很有必要的。相关研究能够帮助构建有效的模型，对CO2地质封存中盖层发生水力破裂的风险概率进行全面准确预测，为盖层力学完整性评价及CO2安全、高效封存提供指导。

5 结 论

(1) 用于CO2地质封存数值模拟研究的力学本构模型有弹性、塑型以及非连续介质模型3类。当采用弹塑性和黏塑性模型进行分析时，地层初始呈拉伸应力状态不利于盖层的力学完整性。

(2) 用于模拟CO2地质封存的耦合分析方法有全耦合、弱耦合以及单向耦合3类，其中，弱耦合应用最为广泛。相应的流固耦合过程分析数值模拟软件有TOUGH2、ECLIPSE、FLAC3D等，基于TOUGH2-FLAC3D的渗流-应力耦合分析在盖层力学完整性研究中应用较为广泛。

(3) 针对CO2地质封存力学问题模拟，常用的数值计算方法有适用于连续介质的、适用于非连续介质的，以及用于连续-非连续介质的数值计算方法。其中，基于连续介质的数值方法在模拟岩石破坏的过程，特别是破坏前的损伤演化分析中非常有效。

(4) CO2注入过程中，储、盖层岩石的孔隙流体压力响应对空间离散较敏感；储层流体压力随注入时间增加而增大，随径向距离增大而减小；储层压力抬升很快会传递至盖层，储、盖层交界处最易发生破裂；盖层渗透率对压力传播有重要影响，而孔隙度、渗透率、非均质性、模型边界范围设置等均对压力增量有一定影响。

(5) 深部咸水层封存CO2是地质封存技术中最有前景的途径，从总体现状来看，流固耦合作用下盖层力学完整性数值模拟研究取得了一定成果，但有关多种耦合机理和影响因素共同作用下盖层应力应变本构关系的模拟研究、模型中边界条件的设置与参数的准确获取，以及模型精度的提高等都是后续研究中有待进一步探讨的问题。