CO₂大规模地质封存通过将工业过程、能源利用或大气中分离的CO₂永久性注入地下地质构造中，如深部咸水层和枯竭油气田等，实现CO₂长期与大气的隔绝，是实现“碳达峰”“碳中和”目标的关键技术^[1-2]。页岩气储层致密，纳米级孔隙发育，渗透率普遍低于0.001 mD，是CO₂地质封存的理想介质。我国页岩气资源量丰富，有利储层分布范围广，具有巨大的CO₂地质封存潜力^[3]。

页岩气储层具有源储一体化的特征，页岩有机质生成的烃类气体可大量赋存于干酪根基质内部，干酪根内部气体含量在页岩气总储量中的贡献可超过50%^[4-5]。由于干酪根基质纳米孔隙丰富，烃类气体主要以吸附态的形式赋存^[6-7]。页岩气藏废弃后，CO₂地质封存过程中注入的CO₂同样可通过吸附作用大量赋存于页岩干酪根基质内部^[8-9]。页岩干酪根具有柔性的孔隙骨架，CO₂吸附能够诱导干酪根发生变形，使得干酪根基质体积及内部孔隙结构变化，进而影响CO₂吸附能力和渗流特征。明确页岩干酪根中CO₂吸附与孔隙变形的耦合特征和微观机理有助于指导页岩气储层中CO₂封存潜力评估。

目前国内外学者主要通过室内物理模拟实验、分子模拟以及理论模型方法研究页岩中CO₂吸附与孔隙变形的耦合作用特征。室内物理模拟实验通过联合等温吸附装置和应变计量装置同时测试样品的吸附量及体积应变量，进而表征吸附-变形耦合特征^[10-11]。该方法主要应用于形状固定且规则的页岩样品，无法适用于粉末状的干酪根样品。在分子模拟方法中，通过耦合巨正则蒙特卡罗算法和分子动力学算法能够实现非定型材料中气体吸附与孔隙变形的耦合过程^[12-13]。部分学者利用该方法研究了干酪根中CO₂吸附与孔隙变形的耦合行为^[14]，但尚未揭示其微观耦合机理。上述两种算法的耦合模拟效率较低，需要很长的模拟时间达到吸附-变形平衡，并且无法实现确定体积应变率下的吸附-变形平衡。BROCHARD等^[15]建立了一套多孔介质力学理论模型，为多孔介质吸附-变形行为的预测提供了新思路。HUANG等^[16]和SUI等^[17]基于该模型，结合固定孔隙中CO₂吸附量分子模拟结果，预测了页岩干酪根的吸附-变形特征。然而，上述研究均采用煤岩有机质的吸附-变形耦合系数预测页岩干酪根的吸附-变形行为，预测结果的合理性有待验证。

目前，页岩干酪根中CO₂吸附与孔隙变形的耦合系数未见报道。为此，本研究针对页岩不同成熟度的干酪根，通过开展不同应变和压力条件下的CO₂吸附分子模拟，结合多孔介质力学理论，确定干酪根中CO₂吸附与孔隙变形的耦合系数，对比不同成熟度干酪根与煤岩有机质之间的差异。通过分析变形干酪根的孔隙结构特征及其与CO₂的能量作用机制阐明吸附-变形耦合机理。

1.   干酪根模型和模拟方法

1.1   干酪根模型构建

构建合理的干酪根模型是在分子尺度上量化干酪根中CO₂吸附量与干酪根体积应变量的前提。为了研究干酪根成熟度的影响，本研究选取了Ⅱ-A、Ⅱ-C和Ⅱ-D三种不同成熟度的干酪根。这三种类型干酪根的分子结构单元(图 1)由UNGERER等^[18]基于KELEMEN等^[19]的实验表征数据构建，化学结构式分别为C₂₅₂H₂₉₄O₂₄N₆S₃(Ⅱ-A)、C₂₄₂H₂₁₉O₁₃N₅S₂(Ⅱ-C)、C₁₇₅H₁₀₂O₉N₄S₂(Ⅱ-D)。

图  1  不同成熟度干酪根的分子结构单元和凝聚态基质模型

Figure  1.  Structural units and condensed matrix models of kerogen with different maturities

下载: 全尺寸图片 幻灯片

基于不同成熟度干酪根的分子结构单元，参照课题组提出的干酪根结构弛豫方法和流程^[20]，构建页岩凝聚态的干酪根基质模型。首先，对干酪根结构单元进行几何优化和退火模拟，得到具有最低能量结构的干酪根结构单元。然后，将优化后的一定数量干酪根结构单元整合到不定形晶胞中，构建密度为0.05 g/cm³的干酪根初始三维模型。考虑模拟经济性和效率性，参考先前研究工作^[21]，本研究在干酪根初始三维模型构建中分别选取6个Ⅱ-A型、7个Ⅱ-C型和12个Ⅱ-D型干酪根结构单元。然后采用NVT系综(恒定原子数、体积、温度)动力学方法，在800 K条件下对干酪根初始三维模型开展400 ps高温弛豫。随后，在20 MPa条件下，采用NPT系综(恒定原子数、压力和温度)控制温度逐步下降，该过程分为800、700、500、300 K四个阶段，每个阶段模拟200 ps。最终，在温度328.15 K和压力20 MPa条件下开展1 ns的干酪根结构弛豫模拟，得到密度收敛、热力学性质稳定的干酪根基质模型，如图 1所示。

1.2   模拟方法和计算细节

为了确定干酪根中CO₂吸附与孔隙变形的耦合特征，需获取干酪根基质在不同应变条件下的吸附量。本研究基于未变形的干酪根基质模型，提出了一套新的模拟方法(图 2)。首先结合干酪根目标体积应变率(-10%~10%)，将模型盒子调整至变形后的尺寸；然后通过在NVT系统下进行1 ns的结构弛豫，得到变形平衡后的干酪根模型；最后基于变形平衡后的干酪根模型，开展CO₂吸附分子模拟，得到干酪根在不同应变下的CO₂等温吸附曲线。

图  2  干酪根中CO₂吸附与孔隙变形耦合行为的模拟流程示意图

Figure  2.  Schematic diagram of the simulation method for coupling CO₂ adsorption and pore deformation in kerogen

下载: 全尺寸图片 幻灯片

干酪根中CO₂吸附与孔隙变形的耦合特征可由耦合系数量化表征。基于多孔介质力学理论，当干酪根中CO₂吸附量与体积应变符合近似线性关系时，CO₂吸附量可以表示为^[15]

其中，n(ε, p)为变形干酪根的CO₂吸附量(mmol/g)；n₀(p)为未变形干酪根的CO₂吸附量(mmol/g)；C(p)为干酪根中CO₂吸附与孔隙变形的耦合系数(无量纲)。

ε为干酪根的体积应变率(%)：

其中，V为变形后的干酪根体积(nm³)，V₀为未变形的干酪根体积(nm³)。干酪根在不同应变条件下的CO₂吸附量采用巨正则蒙特卡罗(GCMC)方法获取，每个压力点下的总模拟步数为4×10⁷步，其中前2×10⁷步用于吸附平衡，后2×10⁷步用于吸附量的系综统计，吸附模拟的温度恒定为328.15 K。

采用PCFF力场描述干酪根模型、CO₂分子及其相互作用。静电作用力和范德华力分别通过Ewald和Atom based方法计算，截断半径为1.55 nm，控温和控压方法分别采用Nose和Berendsen方法，模拟步长为1 fs。

采用Connolly方法表征干酪根的孔隙结构。该方法在干酪根基质中随机插入特定半径的探针分子，探针分子在干酪根基质骨架的van der Waals表面滚动，通过探针的滚动轨迹可确定干酪根骨架表面，干酪根骨架表面之间的区域为自由孔隙体积。结合干酪根中的自由孔隙体积和干酪根骨架表面积，可以计算得到干酪根基质的孔隙度和比表面积：

其中，φ为干酪根孔隙度(%), V_f为干酪根自由孔隙体积(nm³), S为干酪根基质的比表面积(m²/g), S_a为干酪根骨架表面积(nm²), ρ为干酪根模型密度(g/cm³)。基于Connolly方法，通过采用不同半径的探针分子探测干酪根自由孔隙体积，可以确定干酪根的孔径分布。在页岩吸附解吸实验中，通常采用氦气表征页岩纳米孔隙结构^[22]。因此采用氦气的分子动力学直径作为探针直径。

干酪根的膨胀/压缩变形会对干酪根吸附能力造成影响，干酪根中气体吸附能力可以根据相互作用能来表征。吸附质(CO₂)与吸附剂(干酪根)之间的相互作用能的计算公式为

其中，E_int为气体-干酪根的相互作用能，E_total为CO₂-干酪根平衡后的体系能量，E_gas为平衡体系中CO₂分子内部的能量，E_kerogen为平衡体系中干酪根内部的能量。能量单位均为kJ/mol。

2.   结果讨论

2.1   干酪根模型验证

结合干酪根基质的密度、孔隙度和比表面分析干酪根模型的合理性。Ⅱ-A型干酪根模型的密度范围为(1.12±0.02) g/cm³，与新阿尔巴尼页岩未成熟Ⅱ型干酪根的密度范围相近^[23]。Ⅱ-C型和Ⅱ-D型干酪根的密度范围分别为(1.15±0.01) g/cm³和(1.17±0.01) g/cm³。干酪根基质的密度随着成熟度的增加而增加，该趋势与前人报道的趋势^[21]相同。另外，Ⅱ-A、Ⅱ-C和Ⅱ-D型干酪根基质模型的孔隙度分别为(12.96±0.24)%、(17.47±0.32)%、(24.84±0.86)%，与Barnett泥岩有机质的孔隙度范围(4.45%~22.50%)基本接近^[24]。不同成熟度干酪根模型的比表面积介于1 621.34~1 993.06 m²/g，与前人报道的模拟结果(1 512~1 845 m²/g)^[16-17]较为吻合。因此，本研究构建的干酪根基质模型具有合理的密度、孔隙度和比表面积，能够用于CO₂吸附与孔隙变形耦合行为的研究。

2.2   吸附-变形耦合特征

气体吸附会诱导干酪根发生膨胀或压缩变形，而干酪根变形会影响气体吸附能力，干酪根变形与气体吸附通过相互耦合作用直到平衡。干酪根在不同应变和压力条件下的CO₂吸附量如图 3所示。

图  3  干酪根在不同应变和压力条件下的CO₂吸附量

Figure  3.  CO₂ adsorption capacity in kerogen with varying stresses and pressures

下载: 全尺寸图片 幻灯片

Ⅱ-C和Ⅱ-D型干酪根的CO₂吸附量在整个变形范围内与体积应变率呈近似线性关系，CO₂吸附量随体积应变率增加而增大，其吸附-变形耦合规律与不定型链的结果^[15]一致。Ⅱ-A型干酪根在低压缩和膨胀变形范围(ε处于-6%~10%)的CO₂吸附量与体积应变率近似线性相关，与高压缩变形范围(ε处于-10%~-6%)的吸附-变形耦合规律存在差异。高压缩变形下，Ⅱ-A型干酪根基质内部微小孔隙大量闭合，导致干酪根吸附-变形特征显著变化。

不同应变条件下干酪根中CO₂吸附量受压力和成熟度影响。低压条件下(< 4 MPa)，CO₂分子优先占据高能吸附位点，吸附量随压力增加而迅速增大，在更高压力下，干酪根中吸附相空间逐步被CO₂分子填满，吸附量趋于饱和。随着干酪根成熟度增加，CO₂吸附量逐渐增大。在体积应变率10%、气体压力20 MPa条件下，过成熟Ⅱ-D型干酪根内CO₂吸附量达到7.06 mmol/g，而未成熟Ⅱ-A型干酪根内CO₂吸附量为1.90 mmol/g。这是因为随着干酪根热演化程度加深，干酪根芳构化程度加大，微小孔隙尺寸增大，孔隙度增大，使得吸附相空间增加，进而导致干酪根吸附能力增强。

图 4为不同成熟度和应变条件下干酪根基质中的CO₂密度分布剖面。CO₂分子在干酪根基质中非均匀分布，随成熟度增加，高密度区域的CO₂密度值增大，进一步证实干酪根对CO₂的吸附能力增强。当干酪根由压缩变形状态转变为膨胀变形状态，基质内部的CO₂高密度区域位置发生变化，表明干酪根体积应变使得基质内部孔隙结构改变。另外，相比压缩变形状态，膨胀变形状态下干酪根内部CO₂高密度区域的范围呈增大趋势，而CO₂低密度区域呈减小趋势。因为膨胀变形使得干酪根内部微小孔隙尺寸增大，不仅增大了CO₂吸附相空间，也改善了孔隙连通性，增加了CO₂与吸附位点的接触。

图  4  10 MPa下变形干酪根孔隙中CO₂密度分布剖面图

Figure  4.  CO₂ density distribution contour in the pores of deformed kerogen at 10 MPa

下载: 全尺寸图片 幻灯片

吸附-变形耦合系数反映干酪根的应变敏感性和变形能力。根据图 3中CO₂吸附量与干酪根体积应变率的线性段耦合关系，结合等式(2)计算干酪根中CO₂吸附与孔隙变形的耦合系数(图 5)。为了分析不同应变阶段的吸附-变形耦合规律，分别计算干酪根在压缩变形、膨胀变形和压缩-膨胀变形范围的耦合系数。低压条件下(< 2 MPa)，在压缩或者膨胀变形阶段，Ⅱ-A型和Ⅱ-C型干酪根的吸附-变形耦合系数随压力的增加呈增大趋势，更高压力条件下，吸附-变形耦合系数基本不随压力改变。Ⅱ-D型干酪根的吸附-变形耦合系数在压缩或者膨胀变形阶段均不受压力影响，可近似为常数。

图  5  不同变形范围内干酪根的吸附-变性耦合系数

Figure  5.  Sorption-deformation coupling coefficient of kerogen within different strain stages

下载: 全尺寸图片 幻灯片

表 1为3种干酪根在不同变形阶段的吸附-变形耦合系数。随着成熟度增加，干酪根的吸附-变形耦合系数逐渐降低，在相同CO₂吸附量条件下，干酪根的变形程度减小，变形能力下降。这与不同成熟度干酪根的化学结构差异相关，未成熟的Ⅱ-A型干酪根结构中包含大量长脂肪链，增强了其柔性结构特性；Ⅱ-C型干酪根结构中长脂肪链变短，芳香环数量增多，使得其变形能力下降；Ⅱ-D型干酪根芳构化程度加大，结构中以多环芳香簇为主，导致其变形能力较低。

表  1  不同变形区间内干酪根中CO₂吸附与孔隙变形的耦合系数

Table  1.  Coupling coefficient between CO₂ adsorption and pore deformation in kerogen within varying strain conditions

参数	Ⅱ-A	Ⅱ-C	Ⅱ-D
压缩变形	8.22±0.25	7.02±0.49	4.54±0.08
膨胀变形	16.22±1.05	8.18±0.79	4.84±0.19
压缩-膨胀变形	12.26±0.44	7.55±0.24	4.86±0.12

下载: 导出CSV 

| 显示表格

不同成熟度干酪根的吸附-变形耦合系数在膨胀变形范围与压缩变形范围内存在不同程度差异。其中，Ⅱ-A型干酪根在膨胀变形范围内的C(p)约为压缩变形范围内的1.97倍。因此，对于未成熟阶段的干酪根，建议根据其应变状态选择相应的耦合系数。Ⅱ-C型和Ⅱ-D型干酪根在膨胀变形范围内的C(p)与压缩变形范围内的结果接近，压缩和膨胀变形阶段可采用相同的耦合系数。

目前未见页岩干酪根中CO₂吸附与孔隙变形的耦合系数报道，主要采用煤岩有机质中CO₂吸附-孔隙变形耦合系数(图 5中C(p)=7.60)预测页岩干酪根的吸附-变形耦合行为^[15]。根据表 1结果，相比报道的煤岩有机质中CO₂吸附-孔隙变形耦合系数，Ⅱ-C型干酪根中结果较为相近，Ⅱ-A型干酪根中结果偏大，Ⅱ-D型干酪根中结果偏小。因此，页岩干酪根中CO₂吸附与孔隙变形的耦合系数不能直接采用煤岩有机质中结果，基于煤岩有机质结果预测页岩干酪根吸附-变形耦合行为的合理性和准确性有待检验。

2.3   变形干酪根的孔隙结构特征

干酪根变形会导致其纳米孔隙结构发生变化，进而对气体吸附能力造成影响。图 6为干酪根在不同应变下的孔隙结构参数。图 7为干酪根在不同应变下的孔径分布曲线。由图 6可知，干酪根基质的孔隙度和比表面积均与体积应变率呈正相关关系。随体积应变率增加，干酪根孔隙度线性增大，且不同成熟度干酪根孔隙度的增大幅度基本接近。干酪根比表面积在压缩变形和膨胀变形范围内均分别与体积应变率呈线性关系。Ⅱ-D型干酪根比表面积与体积应变率在整个变形范围内具有较好线性关系。相比压缩变形范围，Ⅱ-A型和Ⅱ-C型干酪根的比表面积在膨胀变形范围内随体积应变率增加的幅度降低。随着干酪根膨胀变形，干酪根中微小孔隙数量增加可使得比表面积增大，而微小孔隙尺寸增大可使得比表面积增大幅度减小。

图  6  干酪根在不同应变下的孔隙结构特征

Figure  6.  Porous structural characteristics of deformed kerogen

下载: 全尺寸图片 幻灯片

图  7  干酪根在不同应变下的孔径分布曲线

Figure  7.  Pore size distribution curve of deformed kerogen

下载: 全尺寸图片 幻灯片

不同成熟度干酪根基质中(图 7)，Ⅱ-D型干酪根的孔径范围最大(0.26~1.22 nm)，其中CO₂可吸附的有效孔隙(>0.33 nm)为主要孔隙类型；而Ⅱ-A型干酪根的孔径范围最小(0.26~0.58 nm)，其中CO₂不能进入的无效孔隙(< 0.33 nm)为主要的孔隙类型。随着体积应变率增加，干酪根在压缩变形状态下的微小孔隙主峰降低，更大尺寸的孔隙峰增高，孔径分布曲线向更大孔隙方向偏移。干酪根基质内部微小孔隙逐渐转变为更大尺寸孔隙，通过增大吸附相空间和改善孔隙连通性，增大CO₂吸附能力。

2.4   变形干酪根的能量作用机制

干酪根基质中CO₂分子与孔壁原子的非键相互作用是导致孔隙变形的主要原因。由不同应变下干酪根与CO₂的非键作用能(图 8)可知，范德华作用能、静电作用能以及总非键作用能均为负值，表明CO₂分子与干酪根之间的范德华力作用和静电力作用均能促进CO₂分子在干酪根基质中的吸附。不同类型非键作用能的绝对数值越大，CO₂与干酪根之间的作用强度越大。在3种不同成熟度干酪根中，范德华作用强度、静电作用强度以及总非键作用强度从强到弱的顺序均为：Ⅱ-D型、Ⅱ-C型、Ⅱ-A型。CO₂与干酪根之间作用强度与干酪根体积应变率呈正相关关系。随着CO₂与干酪根之间作用强度增加，干酪根基质孔壁受到的膨胀应力增大，干酪根孔隙度与比表面积相应增加，使得干酪根体积应变率增大。

图  8  干酪根与CO₂的非键作用能与体积应变率的关系(10 MPa)

Figure  8.  Correlation between non-bonded interaction energy and kerogen volumetric strain ratio at 10 MPa

下载: 全尺寸图片 幻灯片

图 8D量化了范德华力作用和静电力作用在CO₂与干酪根之间总非键作用中的贡献。范德华力作用贡献率介于63%~80%，静电力作用贡献率介于20%~37%。CO₂与干酪根之间的相互作用以范德华力为主，同时静电力也发挥着重要作用。一方面，干酪根基质中包含大量微小孔隙，孔壁原子对CO₂分子的范德华势重叠，产生微孔充填效应，使得范德华力作用强度较大。另一方面，CO₂分子具有较大的原子电荷，能够与干酪根中的杂原子官能团发生较强的静电相互作用。对比3种成熟度干酪根，Ⅱ-A型具有更大的范德华力作用贡献率，Ⅱ-C型具有更大的静电力作用贡献率。

3.   结论

本研究采用分子模拟方法获取了页岩干酪根中CO₂在不同应变和压力下的吸附量，结合多孔介质力学理论分析了干酪根中CO₂吸附与孔隙变形的耦合特征，联合孔隙结构和非键作用能揭示了CO₂吸附和孔隙变形的耦合机制。主要结论如下：

(1) 页岩干酪根中CO₂吸附与孔隙变形的耦合系数不能直接采用煤岩有机质中结果。随着干酪根成熟度增加，受干酪根化学结构影响，吸附-变形耦合系数逐渐降低，干酪根变形能力降低。

(2) 干酪根中CO₂吸附与孔隙变形的耦合系数在一定压力范围内(>2 MPa)可近似为常数。未成熟干酪根在压缩变形和膨胀变形阶段需分别确定对应的耦合系数，成熟和过成熟干酪根在压缩和膨胀变形阶段可采用相同的耦合系数。

(3) 随干酪根成熟度增加，CO₂与干酪根之间的吸附作用强度增大，吸附作用以范德华力为主(贡献率介于63%~80%)，同时静电力也发挥着重要作用(贡献率介于20%~37%)。

(4) 干酪根中CO₂吸附使得基质孔隙度和比表面积增大，微小孔隙数量减小，更大尺寸孔隙数量增加，孔径分布曲线向更大尺寸孔隙方向偏移，导致干酪根体积应变率增大。