CO₂驱煤层气地质封存(CO₂-ECBM)通过向富含CH₄的煤储层注入CO₂，能够有效地提高煤层气(CBM)井CH₄产气速度，并将CO₂封存于煤层内，是提高非常规天然气开采、实现“碳达峰”与“碳中和”目标的一条重要路径。当前，世界范围内开展了多次CO₂-ECBM场地试验，有利地证明了CO₂-ECBM的可行性与有效性^[1-2]。自20世纪90年代中期起，美国先后在圣胡安盆地、伊利诺伊盆地、阿巴拉契亚盆地、威利斯顿盆地、黑勇士盆地实施了多次CO₂-ECBM场地试验，加拿大在阿尔伯塔省先后开展了名为FBV 4A和CSEMP的CO₂-ECBM场地试验^[2-3]。此外，德国、波兰和日本也开展了CO₂-ECBM先导性试验^[3-6]。我国从2004年开始在沁水盆地和鄂尔多斯盆地多次实施了不同规模的CO₂-ECBM场地示范试验。其中，最后一次位于沁水盆地柿庄南区块，目标层位为山西组3号煤储层，自2020年6月至2021年6月共向储层内注入CO₂约2 000 t^[1]。

为支持CO₂-ECBM工程实施，国内外开展了许多相关数值模拟工作。凡永鹏^[7]和刘世奇^[8]等基于均质煤储层地质模型和渗流−温度−应力(HTM)多场耦合数学模型分析了CO₂注入煤层后储层孔隙度和渗透率的变化情况，发现受煤基质吸附CO₂膨胀变形和解吸CH₄收缩变形的影响^[9-10]，储层孔隙度和渗透率在CO₂吸附区显著下降，在CO₂未波及区或逐渐升高或先小幅降低再升高。ZHANG等^[11]利用数值模拟发现沁水盆地柿庄北区块SX006井组CH₄采收率在CO₂注入2 a后会增加，整个井组CO₂封存潜力约为12 616 t。FAN等^[12]基于数值模拟研究发现，与初次开采(未向煤层注CO₂)相比，CO₂-ECBM过程中煤层气井的CH₄产气峰值较高但在时间上存在延迟。LIU等^[13]通过数值模拟分析认为以4～8 MPa向柿庄南3号煤层注入CO₂，50 d后单井CO₂存储量可达0.73×10⁵～2.54×10⁵ m³。此外，在CO₂注入优化方面，SAMS等^[14]指出注气优化需要兼顾到CO₂封存量和CH₄产气量；PERERA等^[15]发现CO₂注气压力在CO₂-ECBM工艺优化中起主导作用。

上述模拟研究工作有助于更好地理解CO₂-ECBM过程，为CO₂-ECBM工程实施提供一定帮助，但对CO₂-ECBM实际工程的指导仍存在一些问题，主要表现在：① 地质模型过于简单，缺乏对煤储层构造与属性非均质分布的考虑；② 数学建模不够完善，忽视了CO₂注入对煤储层化学场变化的影响；③ 注气优化不够详实，缺少对CO₂-ECBM实际工程时间和井组关井方式的约束。笔者在前人研究基础上，构建了基于COMSOL Multiphysics的柿庄南3号煤储层三维非均质地质模型，通过CO₂-ECBM数值模拟分析了CO₂注入压力对煤储层CH₄增产、CO₂封存和地球化学环境变化的影响，结合CO₂-ECBM工程时间和井组关井方式，优化了不同工程目标的CO₂注气压力。

1.   基于COMSOL Multiphysics的三维非均质地质模型

1.1   地质建模平台及模型

1.1.1   Petrel地质建模

作为一款专业地质建模软件，Petrel拥有多种先进的地质建模算法，如序贯指示模拟、序贯高斯模拟、和多点地质统计随机模拟等^[16-17]。其建立的地质模型多为三维非均质角点地质模型，这种模型通过网格的节点坐标 (x, y, z) 确定模型的空间位置和构造特征，孔隙度和渗透率等属性参数赋存于每个网格块内。由于网格可以定义成任何立体几何形状，角点地质模型能够较准确地展示地层的三维展布形态和断层、尖灭等构造特征。

基于Petrel平台建立的角点地质模型包含构造模型和属性模型。构造模型数据由每个网格节点的 (x, y, z) 空间坐标组成，网格节点空间排布顺序一般呈垂向纵切、先上后下、先左后右的排列特征。属性模型数据则根据完整网格的顺序依次排列。以图1中2×2×1的模型为例，构造模型的网格节点顺序依次为1(x₁, y₁, z₁) → 2(x₂, y₂, z₂) → 3(x₃, y₃, z₃) → 4(x₄, y₄, z₄) → 5(x₅, y₅, z₅) → 6(x₆, y₆, z₆) → 7(x₇, y₇, z₇) → 8(x₈, y₈, z₈) → 9(x₉, y₉, z₉) → 10(x₁₀, y₁₀, z₁₀) → 11(x₁₁, y₁₁, z₁₁) → 12(x₁₂, y₁₂, z₁₂) → 13(x₁₃, y₁₃, z₁₃) → 14(x₁₄, y₁₄, z₁₄) → 15(x₁₅, y₁₅, z₁₅) → 16(x₁₆, y₁₆, z₁₆) → 17(x₁₇, y₁₇, z₁₇) → 18(x₁₈, y₁₈, z₁₈)，属性模型参数的排列顺序依次为Ⅰ → Ⅱ → Ⅲ → Ⅳ。

图  1  典型角点地质模型

Figure  1.  Typical corner grid geological model

下载: 全尺寸图片 幻灯片

1.1.2   COMSOL Multiphysics地质建模

与专业地质建模软件Petrel不同，COMSOL Multiphysics的功能以数值模拟计算为主，地质建模能力则相对较弱。在三维地质建模方面，该软件缺乏相应的专业建模算法，不能独立完成非均质地质建模工作。尽管其地质建模能力相对较弱，但其内置的一些基础建模方法仍具有较大的意义。在构造建模方面上，COMSOL Multiphysics软件可以直接定义圆柱体、圆台、长方体、球体等多种立方体(图2)。在模型修饰方面，该软件具有回转、拉伸、扫掠、放样、布尔操作与分割、倒角、圆角和删除实体等多种处理方法。此外，它还能将导入的外部数据(如文本数据)通过线性内插、最邻近插值、参数化曲面等方法构建各种数据的三维空间分布特征。

图  2  COMSOL Multiphysics软件内置的立方体模型

Figure  2.  Cube model built into COMSOL Multiphysics software

下载: 全尺寸图片 幻灯片

1.2   拆分−重构建模方法

基于对Petrel平台所构建的角点地质模型特征和COMSOL Multiphysics建模能力的分析，笔者提出拆分−重构法将基于Petrel平台建立的沁水盆地柿庄南区块CO₂-ECBM工程研究区煤储层三维非均质地质模型在COMSOL Multiphysics内进行重构。拆分−重构法的核心思想是将角点地质模型的构造模型和属性模型进行分拆，利用编程处理成COMSOL Multiphysics软件识别的数据体(如文本文件)。然后将数据体导入COMSOL Multiphysics软件，结合软件内置的基础建模方法，重构三维非均质地质模型，具体流程如图3所示。

图  3  拆分−重构三维非均质地质模型流程

Figure  3.  Split-reconstruction process of 3D heterogeneous geological model

下载: 全尺寸图片 幻灯片

1.2.1   模型拆分与数据处理

首先，将柿庄南区块煤层气井的井头、井斜、地层分层、测井曲线解释成果等数据导入Petrel软件，利用序贯高斯模拟方法建立煤储层三维非均质角点地质模型；通过关键点坐标构建Polygen边界将CO₂-ECBM工程研究区分割出来。然后，从Petrel内将分割出来的三维非均质地质模型以.GRDECL的格式导出，并将其格式后缀修改为文本.txt形式，利用编程软件分别将描述构造模型的网格节点坐标和描述属性模型的参数值分别提取并存储。最后，通过编程方法，一方面将网格节点坐标转换成描述储层顶面和底面的坐标点阵，为建立层面模型做准备；另一方面根据网格节点坐标计算每个网格块的中心坐标，并与属性参数值结合构建伴有空间坐标的属性参数点阵。

1.2.2   数据导入与模型重构

首先，将包含有煤层顶、底面坐标点阵的文本文件通过插值函数功能导入COMSOL Multiphysics，使用参数化曲面功能构建煤储层的顶面模型和底面模型(图4)，它们共同构成煤储层的构造模型。然后，再将伴有空间位置信息的属性参数点阵导入软件，通过线性内插方法构建在三维空间上呈非均质分布的属性模型。最后，结合煤储层顶、底层面所占空间范围及各煤层气井的空间位置与井径尺寸，利用COMSOL Multiphysics软件内置的立方体建模、布尔操作、分割和删除实体等一系列几何处理功能，建立煤储层和煤层气井^[18]。经过上述处理后，一个完整的煤储三维非均质地质模型即可构建完成。

图  4  煤储层层面模型

Figure  4.  Surfaces of coal reservoir

下载: 全尺寸图片 幻灯片

1.3   模型重构结果、检验与网格划分

1.3.1   地质模型重构结果

利用拆分−重构法在COMSOL Multiphysics内重构的CO₂-ECBM研究区煤储层非均质地质模型如图5和图6所示^[1]。模型内共有9口煤层气井，其中W3井是CO₂-ECBM工程的CO₂注入井，其余井为生产井。煤储层厚度主要为5.00～7.00 m，平均6.26 m，厚煤层主要分布于东北部地区，局部厚度超过7.00 m；薄煤层主要位于模型东南部和西北部边缘地区，最薄处低于5.00 m(图5)。煤储层CH₄含气量主要为5.00～25.00 m³/t，平均14.12 m³/t；W2、W4、W7和W9井所在位置附近区域煤储层CH₄含气量较高，局部地区超过18 m³/t(图6(a))。CH₄吸附时间主要为20～50 d，平均31.35 d；高吸附时间区域主要集中于模型南部及边界，低吸附时间主要分布于模型的中偏北部(图6(b))。

图  5  煤储层厚度变化

Figure  5.  Change of coal reservoir thickness

下载: 全尺寸图片 幻灯片

图  6  煤储层属性模型

Figure  6.  Property model of coal reservoir

下载: 全尺寸图片 幻灯片

孔隙度主要为4.00%～6.00%，平均5.03%。相较于其他位置，模型中部区域煤储层的孔隙度较大，W5井所处位置及其周围区域的孔隙度超过5.55%。模型北部的W6井所处位置及周边区域煤储层孔隙度较低，甚至低于4.50%(图6(c))。由于缺乏煤岩真实渗透率的测试结果，没有在Petrel内直接对渗透率进行建模。但可以根据孔隙度与渗透率之间的数学关系(式(1))^[19-20]，初步确定渗透率的分布情况。

其中，$ k $为煤储层渗透率，10⁻¹⁵ m²；$ {k_{\text{b}}} $为煤储层基准渗透率，暂取10⁻¹⁵ m²；$ \phi $为煤储层孔隙度；$ {\phi _{\text{a}}} $为基准孔隙度，取0.05。由于渗透率与孔隙度呈立方关系，其与孔隙度的分布具有相同的特征(图6(d))。

煤岩弹性模量主要为2.50～3.50 GPa，平均2.80 GPa；高弹性模量主要分布于东部和南部边缘处，最大值超过3.60 GPa(图6(e))。泊松比主要为0.33～0.38，平均0.36；高值泊松比占据了模型的大部分区域，局部地区超过0.38，低值区呈零星斑状分布于模型的东北、东南、西南等区域(图6(f))。

1.3.2   重构模型检验

利用拆分−重构法在COMSOL Multiphysics内重构的煤储层三维非均质地质模型既能描述煤储层在地下空间的三维起伏情况，又能描述储层内属性参数的非均质变化情况。为了检验重构模型的准确性，将角点地质模型与重构地质模型的模型参数进行对比分析(表1)。结果显示2种模型中的煤储层厚度、CH₄含气量、CH₄吸附时间、孔隙度、渗透率、弹性模量和泊松比的最小值、最大值和平均值差异非常小或无差异，煤储层厚度、孔隙度、渗透率、弹性模量和泊松比的平均值相对误差为0，煤储层CH₄含气量和CH₄吸附时间的平均值相对误差分别为0.007%和0.070%，小于0.1%。这表明利用拆分−重构法构建的煤储层三维非均质地质模型能够高度还原基于Petrel平台建立的角点地质模型。

表  1  角点地质模型与重构地质模型参数对比

Table  1.  Parameter comparison between corner grid geological model and reconstructed geological model

项目	模型来源	最小值	最大值	平均值
煤储层厚度/m	角点地质模型	4.912	7.304	6.259
	重构地质模型	4.911	7.305	6.259
煤储层CH4含气量/(m3·t−1)	角点地质模型	5.927	25.639	14.119
	重构地质模型	5.926	25.640	14.118
CH4吸附时间/d	角点地质模型	6.932	62.391	31.353
	重构地质模型	6.930	62.391	31.331
孔隙度/%	角点地质模型	3.938	6.252	5.033
	重构地质模型	3.937	6.253	5.033
渗透率/10−15 m2	角点地质模型	0.489	1.955	1.020
	重构地质模型	0.488	1.956	1.020
弹性模量/GPa	角点地质模型	2.205	3.798	2.804
	重构地质模型	2.205	3.799	2.804
泊松比	角点地质模型	0.326	0.384	0.359
	重构地质模型	0.326	0.384	0.359

下载: 导出CSV 

| 显示表格

1.3.3   重构模型网格划分

地质模型网格划分质量直接影响数值模拟结果的准确性、模型的运算速度以及计算机的内存需求，过于稀疏或过于密集的网格均不利于CO₂-ECBM数值模拟工作的开展。以往模拟工作中，由于地质模型在几何上多为规则的二维矩形或三维长方体，属性上多为均质模型，常将地质模型进行对称简化，取整个模型的1/4进行网格划分^[21-23]。然而，重构的煤储层地质模型垂向上呈三维起伏变化，储层内属性参数呈非均质分布，模型内有9口间距不等的煤层气井(图7(a))，使得在网格划分时不能将其进行对称简化。

图  7  CO₂-ECBM研究区煤储层煤层气井分布及地质模型网格划分

Figure  7.  Distribution of CBM wells and grid division of coal reservoir geological model for CO₂-ECBM study area

下载: 全尺寸图片 幻灯片

根据重构地质模型的构造和属性特征及煤层气井的分布情况，在COMSOL Multiphysics软件内通过自由三角形和扫掠定义功能对其进行网格划分(图7(b))。定义自由三角形最小单元为尺寸0.02 m、最大单元尺寸为40 m，最大单元增长率为1.50，曲率因子0.60。经过网格划分后的重构地质模型拥有8 370个域单元、864个边单元、80个顶点单元；最小网格质量为0.40，平均网格质量为0.75。网格特征总表现为：在煤层气井附近较为密集，在模型边界处较为稀疏。重构模型的网格具有较高的质量和较少的网格数量，既能保证数值模拟结果的准确性，又能提高模型的运算速度和降低计算机的内存需求，有利于后续相关数值模拟研究的开展。

2.   CO₂-ECBM数值模拟

2.1   数学模型

当前，用于CO₂-ECBM数值模拟研究的数学模型多为渗流−温度−应力(Hydraulic-Thermal-Mechanical，HTM)耦合模型^[21-23]，它可以有效描述CO₂-ECBM过程中流体在煤储层内的赋存与运移以及煤储层温度与应力的状态变化。然而，HTM模型欠缺对CO₂注入煤层后引发储层地球化学环境变化的考虑。在前人研究基础上，笔者使用一个自主推导且获得验证的能更全面表征CO₂-ECBM过程的吸附−渗流−温度−应力−化学(Adsorption-Hydraulic-Thermal-Mechanical-Chemical，AHTMC)多场耦合模型^[24]。该模型包含吸附场(A)、渗流场(H)、温度场(T)、应力场(M)、化学场(C)及储层孔−渗方程，不仅考虑了CO₂-ECBM工程中煤储层内的气体吸附、流体渗流、温度和应力变化，还考虑了CO₂注入煤储层后引起的一系列地球化学反应，如CO₂在煤储层裂隙水中溶解形成的碳酸(H₂CO₃)电离生成碳酸氢根离子(${\rm{HCO}}_3^-$)与氢离子(H⁺)、碳酸氢根离子电离生成氢离子(H⁺)与碳酸根离子(${\rm{CO}}_3^{2-} $)、水电离生成氢离子(H⁺)与氢氧根离子(OH⁻)以及方解石在水中溶解生成碳酸根离子(${\rm{CO}}_3^{2-} $)与钙离子(Ca²⁺)(式(2))。同时，该模型还能够描述这些化学离子在储层裂隙水内的运移(式(3))。

式中，${C_{{j}}}$为化学离子组分$j$在裂隙水中的浓度，mol/m³；$ {S_{\text{w}}} $为裂隙含水饱和度；$ t $为时间，d；$ {v_{\text{w}}} $为水的渗流速度，m/s；$ \tau $为孔隙迂曲度；$ {D_{{j}}} $为化学离子组分$j$在水中的扩散系数，m²/s；$ {Q_{{j}}} $为化学反应源汇相，mol/(m³·s)。

此外，CH₄、CO₂和水在储层内赋存与运移时，储层孔隙度与渗透率呈动态变化。随着气体在煤基质内吸附/解吸、储层压力动态变化、储层温度波动以及方解石溶解/沉淀，煤储层孔隙度和渗透率的变化规律可以表示为

式中，$ \alpha $为Biot有效应力系数；$ {\varepsilon _{\text{v}}} $为煤储层应变；下角0代表初始状态；$ {P_{\text{f}}} $为煤储层孔隙压力，MPa；$ K_{\text{s}} $为煤基质块体积模量，GPa；$ {\varepsilon _{\text{s}}} $为煤基质吸附气体膨胀应变；$ {\alpha _{\text{T}}} $为煤岩热膨胀系数，K⁻¹；$ T $为储层温度，K；$ {\varepsilon _{\text{m}}} $为储层固体部分方解石浓度，mol/m³；$ {M_{{\text{vm}}}} $为方解石摩尔体积，m³/mol。

2.2   关键参数与模拟方案

2.2.1   关键参数

笔者使用的关键参数主要源于工程测试、室内实验和前人研究文献。根据室内实验测试，沁水盆地3号煤储层CO₂兰氏体积约为CH₄兰氏体积的1.8倍，CO₂初始吸附时间约为CH₄初始吸附时间的1.5倍。同时，假设储层CH₄含气量的非均质分布受煤基质对气体吸附能力的非均质性控制。用于数值模拟研究的其他关键参数见表2。

表  2  地质模型主要关键参数^{[1, 24]}

Table  2.  Main key parameters of geological model^{[1, 24]}

关键参数	数值		关键参数	数值
初始温度/K	298.15		孔径分布指数	1.00
初始H+浓度/(mmol·L−1)	4.57×10−6		煤岩密度/(kg·m−3)	1450
初始${\rm{HCO}}_3^- $浓度/(mmol·L−1)	5.76		水的密度/(kg·m−3)	1000
初始${\rm{CO}}_3^{2-} $浓度/(mmol·L−1)	1.40×10−2		有效应力系数	0.50
初始OH−浓度/(mmol·L−1)	2.19×10−3		CH4导热系数/(W·(m·K)−1)	3.01×10−2
初始Ca2+浓度/(mmol·L−1)	0.32		CO2导热系数/(W·(m·K)−1)	1.37×10−2
CH4兰氏压力/MPa	2.15		水导热系数/(W·(m·K)−1)	5.99×10−2
CO2兰氏压力/MPa	1.40		煤岩导热系数/(W·(m·K)−1)	0.19
压力系数/MPa−1	7.1×10−2		CH4比热容/(kJ·(kg·K)−1)	2.22
温度系数/K−1	2.10×10−2		CO2比热容/(kJ·(kg·K)−1)	0.84
CH4动力黏度/(Pa·s)	1.34×10−5		水比热容/(kJ·(kg·K)−1)	4.19
CO2动力黏度/(Pa·s)	1.84×10−5		煤比热容/(kJ·(kg·K)−1)	1.35
水动力黏度/(Pa·s)	1.01×10−5		CH4吸附热/(kJ·mol−1)	16.40
CH4指数常数	−15.07		CO2吸附热/(kJ·mol−1)	19.20
CO2指数常数	−12.46		CH4溶解热/(kJ·mol−1)	14.12
CH4温度常数/K	1 545.00		CO2溶解热/(kJ·mol−1)	19.99
CO2温度常数/K	2 258.50		残余水饱和度	0.05
标准温度/K	273.15		残余气饱和度	0
标准压力/MPa	1.01×10−2		方解石中性活化能/(kJ·mol−1)	23.50
CH4扩散活化能/(kJ·mol−1)	34.00		方解石酸性活化能/(kJ·mol−1)	14.40
CO2扩散活化能/(kJ·mol−1)	18.00		中性反应速率常数/(mol·(m2·s)−1)	1.55×10−6
摩尔气体常数/(J·(mol·L)−1)	8.31		酸性反应速率常数/(mol·(m2·s)−1)	0.50
气体滑脱系数/MPa	0.75		方解石比表面积/(m2·kg−1)	9.8×10−3
排驱压力/MPa	1.00×10−2		方解石密度/(kg·m−3)	2710

下载: 导出CSV 

| 显示表格

2.2.2   数值模拟方案

CO₂-ECBM工程开展前各煤层气井已经开采多年，且前述构建的三维非均质地质模型中缺少准确的渗透率。在开展CO₂-ECBM数值模拟研究前，有必要通过生产历史拟合确定真实的储层初始渗透率，并进一步检验三维非均质地质模型和AHTMC数学模型的准确性。基于此分析，根据工程先后顺序将整个模拟分为生产历史拟合和CO₂-ECBM模拟2个阶段(表3)。由于煤储层含气量、孔隙度和渗透率等经过多年生产后已经发生很大变化，在CO₂-ECBM模拟时需要继承生产历史拟合的地质模型。

表  3  数值模拟方案

Table  3.  Scheme for numerical simulation

模拟阶段	排采压力/MPa	注入井	CO2注入压力/MPa	CO2注入温度/℃	模拟时长/a
生产历史拟合	实测	—	—	—	实测
CO2-ECBM模拟	0.1	W3	5.5、6.5、7.5、8.5、9.5	20	15

下载: 导出CSV 

| 显示表格

2.3   生产历史拟合结果与分析

图8展示了各煤层气井CH₄实测产气速度和模拟产气速度的对比情况。各井CH₄的生产拟合曲线均能够较好地匹配实测生产曲线，拟合误差为7.56%～18.67%，平均为12.14%。拟合效果最好的井为W8井，该井在约1 100 d之前产气速度低于500 m³/d，1 100 d后产气速度快速增加至约900.00 m³/d后逐渐降低。煤层气井的生产历史拟合结果表明，构建的三维非均质地质模型和AHTMC多场耦合数学模型对柿庄南煤储层的三维空间展布、属性参数分布和储层内气体吸附、渗流等规律的描述较为准确。

图  8  煤层气井CH₄生产历史拟合

Figure  8.  CH₄ production history matching of CBM wells

下载: 全尺寸图片 幻灯片

CH₄生产历史拟合是根据收集到的生产排采数据，通过调整式(1)中的基准值$ {k_{\text{b}}} $进行的。当$ {k_{\text{b}}} $取1.146×10⁻¹⁵ m²时，基于数值模拟获得的煤层气井CH₄产气速度与实测产气速度的变化趋势基本一致，误差较小(图8)。此时，煤储层初始渗透率为0.488 5×10⁻¹⁵～2.252 2×10⁻¹⁵ m²，平均1.165×10⁻¹⁵ m²(图9)。

图  9  煤储层初始渗透率分布

Figure  9.  Initial permeability distribution of coal reservoir

下载: 全尺寸图片 幻灯片

2.4   CO₂-ECBM模拟结果与分析

2.4.1   CH₄生产与增产

图10展示了自然生产和向煤储层注CO₂时W6和W7井的CH₄产气速度变化情况。在自然生产情况下，W6和W7井CH₄产气速度均表现出随生产时间延长而逐渐降低的变化趋势。与自然生产相比，2口井的CH₄产气速度在CO₂注入初期没有明显提升，随着时间的延长，其数值逐渐增加，且增幅不断变大。当较自然生产有明显提升后，W6和W7井的CH₄产气速度随CO₂注入压力的增加而显著升高。此外，与W7井相比，W6井的CH₄产气速度明显偏低，这主要是由于煤储层CH₄含气量、孔隙度和渗透率的非均质分布造成的。W6井所处位置及周围煤储层CH₄含气量、孔隙度和渗透率远小于W7井，缺乏充足的CH₄气源和较好的渗流通道。

图  10  生产井CH₄产气速度

Figure  10.  CH₄ production rate of CBM wells

下载: 全尺寸图片 幻灯片

表4展示了CO₂-ECBM工程实施15 a后不同CO₂注入压力下生产井组CH₄增产情况。5.50～9.50 MPa注气压力下，CH₄增产率为4.93%～14.39%。随着CO₂注入压力的增加，CH₄增产量及增产率均显著增加。这表明，在CO₂-ECBM过程中，较高的CO₂注入压力更有利于CH₄增产。

表  4  不同注入压力下井组CH₄增产情况

Table  4.  CH₄ stimulation of well group under different injection pressures

CO2注入压力/MPa	5.5	6.5	7.5	8.5	9.5
CH4增产量/106 m3	1.11	1.52	2.00	2.57	3.23
CH4增产率/%	4.93	6.76	8.93	11.46	14.39

下载: 导出CSV 

| 显示表格

2.4.2   CO₂注入与封存

图11展示了在不同注气压力下注入井W3的CO₂注入速度变化情况。总体上，CO₂注入速度表现为随注气压力的增加而升高，随注入时间的延长先逐渐降低再趋于稳定。根据达西定律，储层内流体运移速度主要受控于流体压力梯度和渗透率。CO₂-ECBM过程中，CO₂注入压力越大，则储层内的CO₂压力梯度越大，使得CO₂的渗流速度越大；吸附于煤基质的CO₂诱导煤岩膨胀，随着时间的延长，煤基质膨胀逐渐增加并趋于稳定，导致储层孔隙度和渗透率逐渐降低并趋于稳定。

图  11  W3井CO₂注入速度

Figure  11.  CO₂ injection rate of W3 well

下载: 全尺寸图片 幻灯片

经W3井注入煤储层的CO₂大部分吸附于煤基质内。图12展示了不同CO₂注入压力(P_in)下的煤储层基质内CO₂含气量的动态变化情况。由于煤岩对气体的吸附能力是非均质的，储层基质内的CO₂含气量在CO₂-ECBM过程中是呈非均质分布的。随着时间的延长和CO₂注入压力的增加，煤基质内CO₂含气量逐渐增加并由注入井向周围不断扩展，高含气量区域面积(35～50 m³/t)也不断扩大。

图  12  不同CO₂注入压力下煤储层基质内CO₂含气量变化

Figure  12.  Change of CO₂ content in coal matrix with different CO₂ injection pressures

下载: 全尺寸图片 幻灯片

部分注入煤储层的CO₂在压力驱动下经裂隙运移至生产井井底并产出，图13展示了W6和W7井CO₂产气速度变化情况。在CO₂注入的初始6 a内，2口井几乎没有CO₂产出，这是由于注入储层的CO₂大部分被煤基质吸附。此后，随着注入时间的延长，W6和W7井CO₂产气速度逐渐增加，并随着注气压力的增加而升高。此外，2口井的CO₂产气速度存在明显差异，同一时刻W6井的CO₂产气速度明显高于W7井，这主要是由煤储层的非均质性造成的。三维空间上，2口井与注入井W3所处位置的煤层倾角非常小，基本处于同一水平；平面上，W6与W3的距离较W7与W3的距离更远一些，如果没有储层非均质性的影响，相同时刻W6的CO₂生产速度应该偏低^[24]。然而，受W7所处位置周围煤储层对气体吸附能力较W6井高出很多的影响，CO₂在向W7井运移的过程中会被煤基质大量吸附，促使W7井的CO₂产气速度低于W6井。

图  13  生产井CO₂产气速度

Figure  13.  CO₂ production rate of CBM wells

下载: 全尺寸图片 幻灯片

注入压力对注入井CO₂注入速度和生产井CO₂产出速度的影响直接控制CO₂在煤储层内的封存。CO₂-ECBM工程实施15 a后，随着CO₂注入压力的增加，CO₂累计注入量、累计产出量以及封存量均增加，最大可分别达约8.71×10⁷、3.95×10⁶、8.31×10⁷ m³。受CO₂从生产井产出的影响，CO₂封存率随注气压力增加而减小，但各注气压力下其均保持在95%以上(表5)。这表明，柿庄南3号煤储层拥有较好的CO₂封存能力，低压注入有利于提高CO₂在煤储层内的封存率，高压注入有利于提高CO₂在煤储层内的封存量。

表  5  煤储层CO₂注入与封存情况

Table  5.  CO₂ injection and storage of coal reservoir

注入压力/MPa	累计注入量/m3	累计产出量/m3	封存量 /m3	封存率/%
5.5	32 853 745	745 129	32 108 616	97.73
6.5	44 035 152	1 208 946	42 826 206	97.25
7.5	56 801 441	1 863 003	54 938 437	96.72
8.5	71 155 914	2 757 821	68 398 093	96.12
9.5	87 100 860	3 953 589	83 147 271	95.46

下载: 导出CSV 

| 显示表格

2.4.3   煤储层地球化学环境变化

部分注入煤储层的CO₂会溶解于裂隙水中，并与水中已有的离子发生化学反应，改变储层的地球化学环境。H⁺浓度常被用于判断煤储层裂隙水的酸碱性，298.15 K时，当H⁺浓度超过1.0×10⁻⁷ mol/L时，裂隙水呈碱性；当H⁺浓度小于1.0×10⁻⁷ mol/L时，裂隙水呈酸性。图14展示了不同CO₂注入压力下储层内H⁺浓度的变化情况。相较于初始裂隙水的碱性环境(H⁺初始浓度为4.57×10⁻⁹ mol/L)，注入井附近储层裂隙水在CO₂注入后快速转变为酸性(H⁺浓度大于5.0×10⁻⁵ mol/L)，且随CO₂注入压力的升高，含酸性裂隙水的煤储层范围逐渐增大。随着时间的延长，含酸性裂隙水的煤储层区域不断向外扩展并逐渐波及到生产井，可以预见，生产井产出水中的H⁺浓度也会逐渐增加，需要对其进行处理才能排放。

图  14  不同CO₂注入压力下煤储层裂隙水内H⁺浓度变化

Figure  14.  Change of H⁺ concentration in water of coal fracture with different CO₂ injection pressures

下载: 全尺寸图片 幻灯片

3.   CO₂注入优化

CO₂-ECBM的主要目的是提高煤层气井CH₄产量并将CO₂封存于煤储层中。在实际工程中，随着生产井产出气中CO₂体积分数不断增加，CO₂封存效果会逐渐变差^[25-26]，不得不关闭生产井。此外，煤层气井的寿命存在上限，并非无限延长。这些使得CO₂-ECBM的工程效果会受到一定制约。因此，需要结合时间因素，综合CH₄增产和CO₂封存效果，优化CO₂注入以促使CO₂-ECBM的工程效果达到最优。

3.1   注入优化方法

3.1.1   目标函数构建

工程耗时、CH₄增产和CO₂封存是评价CO₂-ECBM工程效果的基本指标。一般情况下，工程耗时越短、CH₄产量越大、CO₂封存量和封存率越大，CO₂-ECBM工程效果越好。基于此，构建包含工程时间(CO₂开始注入煤储层至所有生产井关闭时CO₂-ECBM工程所消耗的时间)、CH₄产气量、CO₂封存量和封存率等指标在内的目标函数用以评价CO₂-ECBM工程效果(式(5))：

其中，$ {V_{\text{a}}} $为CO₂-ECBM工程效果值；$ {t_{\text{p}}} $为CO₂-ECBM工程时间，d；$ {q_{\text{p}}} $为CH₄产气量，m³；$ {q_{\text{s}}} $为CO₂封存量，m³；$ {r_{\text{s}}} $为CO₂的封存率，%；${\left[ {a - {\text{min}}\left( a \right)} \right]} \mathord{\left/ {\vphantom {{\left[ {a - {\text{min}}\left( a \right)} \right]} {\left[ {{\text{max}}\left( a \right) - {\text{min}}\left( a \right)} \right]}}} \right. } \left[ {\text{max}}\left( a \right) - {\text{min}}\left( a \right) \right]$为对基本指标$ a $进行线性归一化，以便统一所有指标的单位^[27]；$ {\beta _1} $、$ {\beta _2} $、$ {\beta _3} $分别为工程时间、CH₄生产、CO₂封存的权重系数，$ {\beta _1} + {\beta _2} + {\beta _3} = 1 $。在相同的权重系数下，$ {V_{\text{a}}} $越大表示CO₂-ECBM工程效果越好。结合CO₂-ECBM工程对时效性、CH₄增产和CO₂封存3个基本指标的侧重，定义4种CO₂-ECBM工程目标：① 省时提效、CH₄增产和CO₂封存综合兼顾；② 省时提效；③ CH₄增产；④ CO₂封存。其目标函数对应的各参数权重取值见表6。

表  6  目标函数权重系数

Table  6.  Weight coefficient of objective function

工程目标	β1	β2	β3
综合兼顾	1/3	1/3	1/3
省时提效	3/5	1/5	1/5
CH4增产	1/5	3/5	1/5
CO2封存	1/5	1/5	3/5

下载: 导出CSV 

| 显示表格

3.1.2   样本采集

研究区煤储层局部地区的最小主应力不超过9.5 MPa，持续高压注入CO₂很可能会压裂煤储层并在注入井和生产井间建立高渗通道，导致注入煤储层的CO₂沿裂缝直接运移至生产井并产出^[28]。而过低的注气压力不能有效将CO₂注入煤储层，无法取得较好的CO₂-ECBM工程效果。笔者将CO₂注入压力限制在3.50～9.50 MPa。除选择3.50、4.50、5.50、6.50、7.50、8.50、9.50 MPa等7个注气压力外，还利用均匀随机分布函数模拟产生了30个3.50～9.50 MPa的其他注气压力。

3.1.3   关井方式

由于研究区井组有8口生产井，有必要确立生产井的关井方案。笔者以生产井产出气中CO₂体积分数达到10%为阈值设立以下2种关井方式：① 同时关井，指以整个研究区8口生产井为研究对象，当8口产气井瞬时累计产出气的CO₂体积分数达到10%时，同时关闭所有井；② 序贯关井，指在CO₂-ECBM过程中单井产出气中CO₂体积分数达到10%时，关闭该生产井，当所有生产井都关闭后，CO₂-ECBM工程结束。

3.2   同时关井CO₂注入优化

3.2.1   井组产出气CO₂体积分数与CH₄产气速度

同时关井方式下，CO₂注气压力对生产井组产出气中CO₂体积分数变化影响较大。当CO₂从生产井产出后，井组产出气中CO₂体积分数随注气时间延长不断增加，并随CO₂注入压力增加而升高(图15(a))。至井组产出气中CO₂体积分数达到10%时，CO₂注入时间由小到大对应的注气压力依次为9.5、6.5、3.5 MPa。

图  15  生产井组产出气CO₂体积分数及CH₄产气速度变化

Figure  15.  Change in CO₂ concentration and CH₄ production rate of gas produced by well groups with different CO₂ injection pressures

下载: 全尺寸图片 幻灯片

CO₂注气压力的差异也影响生产井组CH₄的生产。图15(b)展示了3.5、6.5和9.5 MPa注气压力下研究区生产井组CH₄产气速度随CO₂注气时间的变化情况。随CO₂注气时间延长，井组CH₄产气速度逐渐降低；随着CO₂注气压力增加，井组CH₄产气速度升高。至CO₂-ECBM工程结束时，3.5、6.5和9.5 MPa注气压力下的井组CH₄产气速度分别约为4 750、4 000和2 500 m³/d，即单井平均CH₄产气速度分别约为593.75、500.00和312.50 m³/d。

图16展示了3.5、6.5和9.5 MPa注气压力下CO₂-ECBM工程结束时煤储层CH₄含气量分布情况。受煤岩吸附气体非均质性的影响，储层CH₄含气量在煤储层内呈非均质分布，数值范围随CO₂注气压力的增加而扩大。工程结束时，3.5、6.5和9.5 MPa注气压力下的储层CH₄含气量分别为0.63～19.47、0.63～22.80、0.63～24.01 m³/t，平均值分别为9.42、11.45、12.31 m³/t。根据初始和工程结束时煤储层CH₄含气量的数值，可以计算出3.5、6.5和9.5 MPa注气压力下储层煤层气采收率分别为33.28%、18.90%和12.81%，即CH₄采收率随注气压力增加而降低。

图  16  CO₂-ECBM工程结束时不同CO₂注气压力下煤储层CH₄含气量

Figure  16.  CH₄ content in coal reservoir with different CO₂ injection pressures at the end of CO₂-ECBM project

下载: 全尺寸图片 幻灯片

3.2.2   基本指标拟合

图17展示了研究区同时关井方式下CO₂-ECBM工程时间、CH₄产气量、CO₂封存量和封存率随注气压力的变化情况。随着CO₂注入压力的增加，工程时间、CH₄产气量和CO₂封存率逐渐升高，而CO₂封存量则先增加后减小。它们之间的关系可以用指数或多项式函数进行表征：

图  17  同时关井情况下基本工程效果指标随CO₂注入压力变化情况

Figure  17.  Change of basic engineering effect index with CO₂ injection pressure under simultaneous shut-in of wells

下载: 全尺寸图片 幻灯片

3.2.3   CO₂注气压力优化

图18展示了基于式(5)和式(6)计算的同时关井情况下4种目标的工程效果值随CO₂注入压力变化情况。以综合兼顾、省时提效和CO₂封存为目标的工程效果值随CO₂注入压力增加先升高后降低；以CH₄增产为目标的效果值则随CO₂注入压力增加持续减小。此外，以省时提效为目标的效果值在注气压力5.37 MPa以下时小于0，说明采用低压向煤储层注入CO₂不会取得好的省时提效工程效果。在实际CO₂-ECBM工程中，对CO₂注入压力的控制不可能非常精细，以0.05 MPa为基础变量求取最优CO₂注入压力。结果显示，综合兼顾、省时提效、CH₄增产和CO₂封存对应的最优CO₂注入压力分别为6.25、8.05、3.50和6.35 MPa。

图  18  同时关井情况下工程效果值随CO₂注入压力变化情况

Figure  18.  Change of engineering effect value with CO₂ injection pressure under simultaneous shut-in of wells

下载: 全尺寸图片 幻灯片

3.3   序贯关井CO₂注入优化

3.3.1   井组产出气CO₂体积分数与CH₄产气速度

序贯关井方式下，生产井内CO₂体积分数达到10%即会被关闭，导致生产井组产出气的CO₂体积分数在CO₂-ECBM过程中剧烈波动(图19(a))。在整个工程周期内，CO₂注入压力为3.5和6.5 MPa时，井组共发生8次关井，3.5 MPa下生产井关井顺序为W5→W9→W2→W6→W7→W4→W8→W1，6.5 MPa下生产井关井顺序为W5→W9→W2→W6→W4→W7→W8→W1；CO₂注入压力为9.5 MPa时，井组共发生7次关井，生产井关井顺序为W5→W9→W2/6(W2和W6同时关井)→W4→W7→W8→W1。序贯关井方式的关井特征表明CO₂注气压力对CO₂-ECBM过程中生产井组的关井次数和顺序存在影响。

图  19  生产井组产出气CO₂体积分数及CH₄产气速度变化

Figure  19.  Change in CO₂ concentration and CH₄ production rate of gas produced by well groups with different CO₂ injection pressure

下载: 全尺寸图片 幻灯片

序贯关井方式下CO₂-ECBM过程中生产井的关闭直接影响井组CH₄产气速度变化(图19(b))。与同时关井相似，序贯关井时井组CH₄产气速度总体上也呈现出随CO₂注入时间增加而降低的变化趋势。然而，序贯关井方式下井组CH₄产气速度的变化是更复杂的，生产井每关闭一次，井组CH₄生产速度就会下降一次，且低压注入下的CH₄产气速度总体上高于高压注入的情况。值得注意的是，当井组内有效生产井数量相同时，CO₂注入压力越高，井组CH₄产气速度越大。例如，在未关井期间，井组CH₄产气速度由大到小对应的注气压力依次为9.5、6.5和3.5 MPa。此外，在序贯关井方式下CO₂-ECBM工程的最后一次关井过程中，3.5、6.5和9.5 MPa注气压力下的W1井CH₄产气速度分别高于220、390和520 m³/d。

与同时关井方式下的煤储层CH₄含气量分布特征相似，序贯关井方式下采用3.5、6.5和9.5 MPa注气压力时的煤储层CH₄含气量在CO₂-ECBM工程结束时也呈非均质分布，数值范围也随CO₂注气压力的增加而扩大(图20)。3.5、6.5和9.5 MPa注气压力下的储层CH₄含气量分别为1.13～16.60、1.13～21.90、1.13～25.60 m³/t，平均值分别为8.24、10.02、10.87 m³/t。根据初始和工程结束时煤储层CH₄含气量的数值，可以计算出3.5、6.5和9.5 MPa注气压力下储层煤层气采收率分别为41.63%、29.03%和23.01%，即CH₄采收率随注气压力增加而降低。

图  20  CO₂-ECBM工程结束时不同CO₂注气压力下煤储层CH₄含气量

Figure  20.  CH₄ content in coal reservoir with different CO₂ injection pressures at the end of CO₂-ECBM project

下载: 全尺寸图片 幻灯片

3.3.2   基本指标拟合

图21展示了序贯关井时CO₂-ECBM工程时间、CH₄产气量、CO₂封存量与封存率随注气压力变化情况。随着CO₂注气压力的增加，工程时间、CH₄产气量逐渐降低，而CO₂封存量和封存率逐渐升高。它们之间的关系可以用指数函数进行表征：

图  21  序贯关井情况下基本工程效果指标随CO₂注入压力变化情况

Figure  21.  Change of basic engineering effect index with CO₂ injection pressure under sequential shut-in of wells

下载: 全尺寸图片 幻灯片

3.3.3   CO₂注气压力优化

基于式(5)和(7)计算了序贯关井情况下4种目标的效果值随CO₂注入压力变化情况(图22)。以综合兼顾、省时提效和CO₂封存为目标的效果值随CO₂注入压力的增加而增加；以CH₄增产为目标的效果值随注入压力的增加而持续减小。此外，与同时关井情况相似，序贯关井时以省时提效为目标的工程效果值在6 MPa以下时小于0，说明在序贯关井时采用低压向煤储层注入CO₂也不会获得较好的省时提效工程效果。综合兼顾、省时提效、CH₄增产和CO₂封存对应的最优CO₂注入压力均处于注气压力区间的端点，分别为9.5、9.5、3.50和9.5 MPa。综合兼顾、省时提效和CO₂封存对应的最优CO₂注入压力相同，说明相同的注气压力能够实现不同的CO₂-ECBM工程目标。

图  22  序贯关井情况下工程效果值随CO₂注入压力变化情况

Figure  22.  Change of engineering effect value with CO₂ injection pressure under sequential shut-in of wells

下载: 全尺寸图片 幻灯片

3.4   CO₂注入优化分析

表7展示了采用最优CO₂注入压力时2种关井方式下不同CO₂-ECBM目标的工程耗时、CH₄产气量、CO₂封存量和CO₂封存率。采用同时关井方式时，工程时间的顺序为CH₄增产 > 综合兼顾 > CO₂封存 > 省时提效；CH₄产气量的顺序为CH₄增产 > 综合兼顾 > CO₂封存 > 省时提效；CO₂封存量的顺序为综合兼顾 > CO₂封存 > 省时提效 > CH₄增产；CO₂封存率顺序为省时提效 > CO₂封存 > 综合兼顾 > CH₄增产。采用序贯关井方式时，由于综合兼顾、省时提效、CO₂封存3种目标的最优CO₂注入压力相同，它们的工程时间、CH₄产气量、CO₂封存量和CO₂封存率也相同，分别为12 262 d、4.819 5×10⁷ m³、1.821 9×10⁸ m³和99.66%；而以CH₄增产为目标的工程时间、CH₄产气量、CO₂封存量和封存率分别为40 218 d、8.613 1×10⁷ m³、0.999 6×10⁸ m³和98.28%。2种关井方式下的各项基本指标特征表明：欲在CO₂-ECBM工程中节省时间，就无法获得较好的CH₄增产效果；欲获得较好的CH₄增产效果，需以较长的工程时间和较差的CO₂封存效果为代价。此外，在最优CO₂注入压力下，相较于各井同时关井，各井序贯关井以较长的工程时间为代价可以获得较高的CH₄产气量、CO₂封存量和封存率，这表明如果CO₂-ECBM工程时间足够长，采用序贯关井方式可以同时获得较好的CH₄增产和CO₂封存效果。

表  7  最优CO₂注入压力下不同工程目标的基本指标

Table  7.  Basic indexes of different engineering objectives under optimal CO₂ injection pressures

工程目标	同时关井					序贯关井
	综合兼顾	省时提效	CH4增产	CO2封存		综合兼顾	省时提效	CH4增产	CO2封存
工程时间/d	5 370	3 085	14 488	5 194		12 262	12 262	40 218	12 262
CH4产气量/107 m3	4.309 5	3.396 8	6.995 8	4.245 8		4.819 5	4.819 5	8.613 1	4.819 5
CO2封存量/108 m3	0.397 2	0.371 0	0.367 6	0.396 4		1.821 9	1.821 9	0.999 6	1.821 9
CO2封存率/%	97.42	98.33	94.95	97.48		99.66	99.66	98.28	99.66

下载: 导出CSV 

| 显示表格

4.   结　　论

(1)利用拆分−重构法在COMSOL Multiphysics内重构了基于Petrel建立的沁水盆地柿庄南CO₂-ECBM研究区煤储层三维非均质地质模型；经检验，重构的地质模型与原始地质模型差异非常小，能够有效表征煤储层的三维空间展布和储层CH₄含气量、CH₄吸附时间、孔隙度、渗透率、弹性模量和泊松比等属性参数的非均质分布特征；经网格划分，煤储层地质模型被划分为8 370个域单元，网格平均质量为0.75。

(2)注气压力对沁水盆地柿庄南3号煤储层CO₂-ECBM存在显著影响，CO₂-ECBM工程实施15 a后，CO₂注入压力越大，CH₄增产效果越明显、CO₂封存量越大、高CO₂含气量和含酸性裂隙水的煤储层范围越大；受生产井CO₂产出的影响，CO₂封存率随注气压力升高而降低；受煤储层CH₄含气量、CO₂吸附能力和孔渗非均质性的影响，W7井比W6井的CH₄增产效果更好且CO₂产出速度更低。

(3) CO₂注气压力对生产井组产出气CO₂体积分数、CH₄产气速度及煤储层CH₄采收率影响显著；以省时提效、CH₄增产、CO₂封存及3者综合兼顾为工程目标的CO₂最优注入压力在采用同时关井方式时分别为8.05、3.50、6.35和6.25 MPa，在采用序贯关井方式时分别9.50、9.50、3.50和9.50 MPa；好的CH₄增产效果需以较长的工程时间和较差的CO₂封存效果为代价；在最优CO₂注入压力下，与同时关井相比，序贯关井通过延长工程时间能够获得较好的CH₄增产和CO₂封存效果。