第一作者简介：张烈辉（1967-），男，四川仁寿人，西南石油大学教授，主要从事油气藏渗流、数值模拟和气藏工程等方面的科研和教学工作。地址：四川省成都市新都区新都大道8号，西南石油大学油气藏地质及开发工程全国重点实验室，邮政编码：610500。E-mail：zhangliehui@vip.163.com

联系作者简介：张涛（1992-），男，四川丹棱人，西南石油大学助理研究员，主要从事微观渗流理论及模拟、气藏工程理论、油气藏数值模拟等方面的科研和教学工作。地址：四川省成都市新都区新都大道8号，西南石油大学油气藏地质及开发工程全国重点实验室，邮政编码：610500。E-mail: tobiascheuing@163.com

项目基金：国家自然科学基金重点基金项目“海相页岩水平井超临界二氧化碳压裂机理与一体化模拟研究”（52234003）；国家自然科学基金优秀青年科学基金项目“页岩气多尺度非线性渗流力学”（52222402）；国家自然科学基金青年项目“CO₂-水-岩系统中溶解-沉淀-沉淀运移机制及微观渗流规律研究”（52304044）

二氧化碳-水-岩作用机理及微观模拟方法研究进展

张烈辉¹，张涛¹，赵玉龙¹，胡浩然^{1, 2}，文绍牧³，

吴建发²，曹成¹，汪永朝¹，范云婷¹

（1. 西南石油大学油气藏地质与开发工程全国重点实验室，成都 610500；2. 中国石油西南油气田公司页岩气研究院，成都 610051；3. 中国石油西南油气田公司，成都 610051）摘要：系统综述CO2-水-岩复杂作用机理、多孔介质反应输运（溶解、沉淀及沉淀运移）微观模拟、CO2-水-岩系统微观模拟3个方面的研究进展，指出目前研究存在的主要问题并提出了关于未来研究方向的建议。CO2注入储集层后，不仅存在常规渗流体系的流动和传质作用，还会产生溶解、沉淀及沉淀运移等特殊物理化学现象，其耦合作用导致多孔介质的孔渗参数变化规律复杂。孔隙尺度的微观渗流模拟，可以得到孔喉三维空间内的详细信息，且能显性观察到多孔介质流-固界面随反应的变化。目前研究主要在复杂作用机理解耦合、多矿物差异性反应表征、沉淀生成机理及表征（晶体成核和矿物脱落）、沉淀-流体相互作用模拟方法、多物理化学过程耦合渗流机制等方面存在局限。未来研究中，需要创新实验方法对“溶解—沉淀—沉淀运移”解耦合，提高矿物地球化学反应相关参数实验测试的准确度，在不同沉淀机理可靠表征的基础上，建立沉淀-流体相互作用模拟方法，并有机耦合各个物理化学过程，最终实现对CO2-水-岩系统中“溶解—沉淀—沉淀运移”的耦合渗流模拟。

关键词：CO2-水-岩系统；溶解；沉淀；沉淀运移；微观模拟；CO2捕集、利用与埋存；碳中和；解耦合

0 引言

近年来碳中和受到世界各国的空前关注^[1]。在众多的碳中和技术中，CO₂捕集、利用与埋存（CCUS）及CO₂捕集与埋存（CCS）是实现碳中和的托底技术^[2]。CO₂-水-岩系统广泛存在于CCUS及CCS相关的各种工业活动中，包括CO₂驱油驱气埋存^[3]、超临界CO₂压裂埋存^[4]、咸水层CO₂埋存^[5]、有水气藏EGR（Enhanced Gas Recovery，提高天然气采收率）等。中国气藏CO₂驱气埋存、油藏CO₂驱油埋存及深部咸水层CO₂埋存潜力均非常大，潜力区域主要分布在松辽、鄂尔多斯、渤海湾、准噶尔等盆地^[1-2]。CO₂注入地下多孔介质的过程中，至少涉及到以下几个方面的基本问题：①多孔介质流动，地下岩石为具有复杂孔喉结构的多孔介质，因而涉及到多孔介质内的流动问题；②传质（多组分传输），地层水离子组分复杂，在与CO₂和岩石的相互作用中，离子组分浓度可能发生变化，即存在多个组分的传输问题；③溶解（非均相化学反应），CO₂可能与岩石部分矿物发生化学溶解反应，导致岩石孔喉结构变化，这属于非均相化学反应问题；④沉淀生成，CO₂与地层水接触后可能导致某些无机盐过饱和结晶析出，此外，溶解作用可能诱发矿物颗粒从主矿物体脱落，生成新的沉淀；⑤沉淀运移，生成的沉淀在水动力的作用下会与流体一起运移，还可能发生沉淀沉积和再启动现象。因此，CO₂-水-岩系统涉及到CO₂、地层水、岩石之间的相互作用，其中，溶解和沉淀作用将导致岩石的孔隙度和渗透率发生改变，而沉淀运移则可能严重影响岩石的渗流能力（如堵塞喉道），这些物理化学过程的耦合作用可能导致储集层物性发生剧烈变化。对这些过程的微观作用机制认识不明确、缺乏有效的模拟表征方法，将导致储集层宏观孔渗参数变化规律预测不准确，进而使得埋存容量、注采能力、埋存封闭性等一系列参数的评价不可靠^[6]，影响CCUS和CCS工业注采方案制定及安全评估的可靠性。溶解—沉淀—沉淀运移是一种孔隙尺度行为，传统岩心渗流实验往往难以直观刻画其作用过程（如离子浓度变化、孔隙形貌演化、沉淀运移追踪等），导致对储集层物性参数变化规律认识不清。基于孔隙尺度微观渗流模拟，可以得到流体、离子、压力等在孔隙空间的三维分布，且能显性观察到多孔介质固体形貌随反应的变化，是连接宏观尺度与微观尺度的桥梁。开展CO₂-水-岩系统的孔隙尺度模拟，有助于确定溶解、沉淀和沉淀运移的发生条件与判据，量化其对渗流规律的影响，为宏观数值模拟提供基础参数。本文以CCUS及CCS技术中出现的CO₂-水-岩系统为背景，系统综述CO₂-水-岩复杂作用机理、多孔介质反应输运（溶解、沉淀及沉淀运移）微观模拟、CO₂-水-岩系统微观模拟3个方面的研究现状，讨论国内外相关研究的优缺点，梳理目前研究存在的主要问题及技术瓶颈，并针对性地提出未来研究方向。

1 CO₂–水-岩复杂作用机理CO₂-水-岩系统是一个耦合了流动、传质、溶解、沉淀、沉淀运移等过程的复杂渗流体系，如图1所示。CO₂注入地下与地层水接触后，首先在气相压力的作用下溶解在水中。CO₂在水中的溶解度与压力呈正比，与温度、地层水矿化度呈反比^[7]。溶解度主要采用Henry常数计算，较常用的还有一些经验公式^[8]。CO₂溶解于水后，主要以自由分子的形式存在，仅部分CO₂与水结合形成H₂CO₃（见（1）式）。CO₂在水中的溶解反应是一个平衡的过程，在没有催化剂时该反应需要很长时间才能达到平衡^[9]。该溶解反应会释放出H⁺，导致溶液的pH值降低，且降低程度随着CO₂溶解度的增加而增加^[10]。在常温常压下，CO₂饱和溶液的pH值为5.6。在酸性环境下，大多数矿物的化学性质不稳定，H⁺与这些矿物进行一系列的化学反应，产生溶解现象，并可能伴随释放出新的成垢离子（如Ca²⁺、Mg²⁺等）。同时，当CO₂溶于水形成H₂CO₃后，由于解离作用，会进一步生成HCO₃^–和CO₃^2-（见（2）式、（3）式），这些离子与成垢离子互相作用，产生沉淀现象。因此，CO₂-水-岩相互作用会导致溶解和沉淀过程同时发生。

图1CO₂-水-岩系统多物理化学过程耦合的输运示意图

CO₂-水-岩石系统作用机理研究主要通过静态溶解实验和动态岩心驱替实验开展^[11]。通过矿物组分分析、扫描成像、流体组分分析等手段，监测反应前后或反应过程中的矿物组分变化、孔喉变化及最终导致的岩心孔隙度、渗透率变化等信息。

1.1 矿物溶解

由于储集层中的许多岩石矿物（方解石、白云石、长石等）和胶结物化学性质不稳定，在与CO₂、地层水形成的弱酸性流体接触时，会在化学反应的作用下发生溶解，如表1所示。部分学者的实验结果表明，CO₂-水-岩反应导致岩石孔隙度和渗透率增加。那金等^[12]、Cui等^[13]、丁茜等^[14]采用碳酸盐岩岩心开展了一系列储集层条件下的溶解实验和岩心驱替实验，结果表明CO₂在地层水的作用下对岩石矿物产生了强烈的溶蚀现象，增加了岩心的孔隙度和渗透率。孔隙度、渗透率增加的现象不仅出现在碳酸盐岩实验中，在砂岩实验中也有发现。李四海等^[15]将砂岩岩心用CO₂水溶液浸泡后，岩心产生了强烈的溶蚀作用，表现为方解石和白云石含量显著降低、溶蚀孔隙数量增多、原孔径变大，最终导致孔隙度和渗透率也相应增大。需要说明的是，不同矿物的反应速率会有几个数量级的差异^[16-17]，因此，目前开展的相关室内实验受到时间尺度限制，一些速率较慢（几百至几千年）的反应在室内实验中可能无法观测到，不利于评估地质时间尺度下CO₂对储集层孔隙度、渗透率的影响。

1.2 矿物沉淀及运移

矿物溶解可能伴随产生次生矿物和溶解诱发的矿物颗粒脱落，这些沉淀可能在水动力作用下运移堵塞孔隙，导致岩石渗透率降低，且影响程度远大于对孔隙度的影响。于志超等^[23]开展了地层条件下（100 ℃，24 MPa）饱和CO₂地层水的驱替实验，发现实验过程中砂岩岩心的孔隙度基本不变但渗透率持续降低，认为是新生成的次生矿物和碳酸盐胶结物溶解释放出的黏土颗粒运移堵塞孔隙造成。Mohamed等^[24]对砂岩的CO₂驱替实验结果表明岩心的渗透率降低了55%，认为有如下两种因素导致岩心渗透率降低：①CO₂-水-岩反应产生沉淀物运移堵塞孔隙；②部分胶结物发生反应溶解导致附着的黏土颗粒运移堵塞孔隙。Chen等^[25]在60 ℃、20 MPa条件下对渗透率为100×10^-3 μm²的砂岩开展饱和CO₂水的驱替实验，发现岩心的孔隙度增加了1.8个百分点、渗透率降低了5.1%，这是碳酸盐胶结物溶解后诱发的颗粒运移引起的。Nunez等^[26]通过串联2个白云岩岩心研究注CO₂后近井筒和远井筒区域的溶解和沉淀特征，发现第1个岩心的孔隙度、渗透率都有所增加，而第2个岩心孔隙度基本不变，第1个岩心所形成的次生矿物运移到第2个岩心产生孔隙堵塞现象，最终使得远井筒区域渗透率降低。Wang等^[27]研究发现CO₂饱和水注入砂岩后，碳酸盐矿物溶解既会导致渗透率增加，也会诱发颗粒运移降低渗透率，地层水盐度较低时，溶蚀作用占主导地位，渗透率增加，而地层水盐度较高时，颗粒脱落运移显著，渗透率降低。

1.3蒸发结晶

当干燥的超临界CO₂注入矿化度较高的储集层时，会导致地层水蒸发进入CO₂，当地层水中的盐浓度超过其溶解度极限时，易结晶析出，生成无机盐沉淀^[28-29]。具体而言，存在如图2所示的4个机理^[30]：①CO₂注入储集层形成两相流，即超临界CO₂相和含盐地层水相。超临界CO₂驱替地层水后，地层水以水膜、水桥等形式形成残余水。在干燥的超临界CO₂持续流动作用下，这些含盐残余水开始蒸发。随着水的蒸发，CO₂的相对渗透率增加，将进一步促进蒸发。②由于驱替前缘的干燥蒸发作用，会与远端高含水区产生毛管压力梯度，当该梯度超过注入压力梯度时，地层水回流至蒸发前缘，导致蒸发持续发生^[31]；③随着水被蒸发进入超临界CO₂中，残余水中的盐浓度增加，也会导致盐从干燥前缘向地层水侧扩散^[32]。一旦盐浓度因蒸发达到其溶解度极限，盐将从溶液中析出，形成沉淀物。④沉淀的盐对盐水有很强的亲和力，诱导远距离的盐水向蒸发前缘移动、蒸发，进一步导致沉淀增加^[33]，类似于毛细管抽吸作用。因此，当干燥的CO₂注入到高矿化度储集层中时，一般从注入端到储集层端会形成3个典型区域：产生大量沉淀的干燥区域、CO₂驱替地层水的两相流区域、地层水单相流动区域。

图2CO₂驱替—蒸发—沉淀机理（据文献[30]修改）大量的室内实验与矿场实践都发现了超临界CO₂注入导致无机盐沉淀的现象^[33-34]，如图3、图4所示。在室内实验方面，Berntsen等^[35]在饱和地层水的圆柱体砂岩样品中注入超临界CO₂（40 ℃, 10 MPa），通过扫描成像和注入压力监测发现，注入井附近区域的岩石基质中有盐析出并产生了堵塞。在矿场实践方面，2020年，SaskPower公司在Aquistore CCS示范项目中发现，由于储集层中地层水矿化度高（330 g/L），在CO₂注入井周边发现了无机盐晶体，导致注入性降低，经分析认为这些盐沉淀是因为CO₂的持续干化形成^[34]。

图3可视化实验中的无机盐沉淀^[33]

图4井下摄像观察到的无机盐沉淀^[34]

2 多孔介质反应输运的微观模拟

相对于室内实验，基于等比例微观模型的数值模拟具有可重复性好、灵活性高、成本低等优势。CO₂-水-岩系统是一种典型的物理化学过程耦合的输运系统，溶解、沉淀及沉淀运移是该系统区别于其他流体运移系统最显著的特征，下面分别对这些机理的微观模拟研究进展予以阐述。

2.1溶解的模拟

对于包含溶解作用的流体输运系统，需要求解的方程包括：流体的流动方程、流体组分的对流扩散方程（传质方程）、固相发生溶解的化学反应方程以及描述固体溶解的几何形状演变方程。学者们提出了一系列方法来求解这些方程（见表2）：第1类为传统计算流体力学（CFD）的离散方法^[36-39]，例如有限体积（FV）和有限差分（FD）方法；第2类为采用拉格朗日方法描述流动与传质的粒子法^[40-42]，例如光滑粒子流体动力学方法（SPH）和移动粒子半隐式法（MPS）；第3类为混合计算方法^[43-44]，例如采用FD求解流动方程、采用随机移动法（RW）刻画溶质的传输^[43]。传统CFD方法的计算量较大，但是近年来随着计算机性能的提高及非结构化网格的成熟，该类方法也逐渐用于模拟反应流体的输运过程。粒子法的无网格特性使其非常适合模拟紊流、自由表面流、多相流等，但对反应流输运模拟的相关报道还较少。对于混合类计算方法，学者们往往利用一些成熟的软件包或者某种方法在计算效率上的优势，采用不同计算方法求解不同的场方程，但是随着计算机性能的提高及算法的成熟，越来越多地采用单一计算方法求解多物理化学场问题。近些年还发展了一种非贴体网格模拟方法—格子玻尔兹曼法（LBM）^[45-48]，这种方法不需要随着各相区域的运移变化而重构或调整计算网格。相对于常规的CFD方法，LBM具有不受限于流体连续性假设、易于考虑微观力学特征、对复杂多孔介质适应性好、并行简单等优势，因此，越来越多学者开始采用LBM模拟多孔介质中的流动和反应过程。

在孔隙尺度的多孔介质反应输运模拟中，首先要解决显性表征流体相与固体相空间分布的问题。若采用传统CFD方法，需要采用合适的技术追踪或捕获相界面，例如水平集法（LS）^[36]、相场法（PF）^[37]、流体体积法（VOF）^[39]和任意拉格朗日-欧拉方法（Arbitrary Lagrangian-Eulerian，ALE）^[38]。当采用这些方法模拟溶解时，需要根据反应不断地更新流固界面或重构网格，从而捕捉固体边界随时间空间的演化。但是，由于多孔介质内固体边界形状复杂，传统数值方法处理变化的固体边界难度较大。LBM是以分子运动论为基础的介观数值方法，能高效处理任意形状的固体边界，无需追踪、捕获或重构复杂的相界面，模拟多孔介质中的复杂流动优势明显，笔者团队将传统LBM拓展至页岩纳米多孔介质流体传输模拟^[49-51]。当采用LBM模拟溶解时，发展了一些简单高效的界面更新法，例如固相平衡临界值法、像元体积法（VOP）等^[45-48]。在化学反应输运模拟中，需要解决的另一个问题即为化学反应对流体组分浓度的影响。当反应发生在溶液体相不同组分之间时（均相反应），通过在LBM方程中添加源汇项处理；当反应发生在溶液组分与固体表面之间时（非均相反应），在LBM中有3种代表性的处理方法：①Kang等^[46-47]在流固界面采用He等^[52]的热力学边界表征非均相反应边界；②Patel等^[53]提出了一种伪均质方法，通过在流固边界相邻的流体节点添加源汇项表征非均相反应对流体组分的影响；③郭照立团队^[54]提出了一种处理非均相反应的局部反应边界格式，该格式在反弹格式的基础上考虑了化学反应消耗与反应可逆的影响。

2.2 沉淀的模拟

在多孔介质的模拟中，沉淀是溶解的一个反过程，在模拟中的处理方法是相似的，包括边界处理方法、化学反应对流体组分浓度影响的处理方法等。因此，沉淀模拟的关键在于对不同物理化学作用导致的沉淀过程予以刻画。对于CO₂-水-岩系统，沉淀的生成主要包含无机盐的析出和非反应矿物的脱落，沉淀在多孔介质中将严重改变流体的流动路径^[55]，必须对沉淀的生成位置精确模拟。因此，需要根据无机盐析出和非反应矿物脱落的机理，深入分析沉淀生成机理，从而实现对沉淀的合理模拟。无机盐的析出是溶液过饱和后结晶核形成和结晶体生长的过程。成核可以在液体内发生，称为均相成核，也可以在液体和固体界面发生，称为异相成核。成核是控制矿物沉淀的初步生长过程，是一个概率性过程，晶体可以在存在类似条件（如表面性质、过饱和度和温度）的任何地方发生成核^[56]。概率性的成核导致沉淀物在多孔介质不同空间生成，影响流体渗流能力，需要采用概率性方法或基于合理物理基础的尺度化表示方法来表征。目前，针对无机盐沉淀的数值模拟中，一般假设无机盐的成核发生在浓度达到一定阈值时^{[48, 57]}，未考虑概率性的成核过程。近年来，基于经典成核理论（CNT）中的指数成核速率方程（见（4）式），建立了一些新的动力学模型和公式来表征成核过程。有学者尝试利用CNT在沉淀模拟中预测概率性的成核行为。Li等^[58]将CNT耦合进连续尺度的模拟器中刻画方解石的沉淀，模拟了CO₂储存过程中井筒水泥的化学变化。Prasianakis等^[59]将孔隙尺度的反应传输求解器与CNT耦合，展示了均相和异相沉淀动力学对多孔介质演化的影响。Hellevang等^[60]通过实验对玄武岩柱上次生矿物生长的空间分布进行了研究，认为需要开发新的概率成核方法表征无机盐的随机成核过程。Fazeli等^[61]建立了一个新的基于CNT的概率成核模型，并将其耦合进孔隙尺度的反应输运求解器中，模拟了多孔介质中的矿物成核和生长，利用单组分矿物反应对不同的过饱和度、生长速率和流速进行了模拟。Masoudi等^[62]进一步拓展该模型来模拟富含CO₂相中微米级无机盐晶体的形成，以及它们对孔隙形态和渗透率变化的影响。因此，针对无机盐的析出，需要模拟无机盐的成核过程，表征成核相与原固相的差异性生长。确定晶体的成核与生长表征方法后，与溶解类似，不仅需要更新无机盐析出后的流-固几何界面，还要通过添加源汇项至浓度扩散方程，考虑沉淀（晶体成核与生长）对流体离子浓度的影响。

非反应矿物的脱落可以分为两种类型：类型Ⅰ，非反应矿物周边矿物溶解导致直接脱落；类型Ⅱ，矿物溶解导致胶结薄弱面在水动力作用下发生剪切脱落，如图5所示。对于类型Ⅰ，可以利用标定函数判断，若非反应矿物在某时间步不与主岩体连接，则认为非反应矿物脱落，形成了沉淀物，例如Liu等^[63]使用Hoshen-Kopelman算法^[64]来判断脱落矿物与主岩体的关系，从而识别反应过程中溶解导致的矿物脱落。对于类型Ⅱ，首先对全局空间进行位置搜索，判断潜在的胶结薄弱点，分析矿物胶结薄弱点的受力情况。当流速达到一定值（临界流速）时，流体对矿物颗粒的拖曳力（F_d）的扭矩足以克服接触力（F_c1）、胶结力（F_c2）和吸附力（F_a）合力的扭矩，则认为矿物发生脱落、形成沉淀物，如图5所示。因此，首先根据不同矿物的反应速率，采用矿物溶解模拟方法开展模拟计算，在每个计算时间步判断是否有类型Ⅰ或类型Ⅱ的沉淀物生成，若有，则对沉淀物进行运移、沉积与再启动的后续模拟。对于类型Ⅱ的模拟目前尚未见到相关报道。

图5反应矿物溶解后导致非反应矿物脱落形成沉淀示意图

2.3沉淀运移的模拟

CO₂-水-岩系统可能产生沉淀，沉淀在水动力的作用下运移，可能堵塞孔喉，严重降低多孔介质渗流能力。沉淀运移与溶液中颗粒运移现象在本质上一致，而颗粒运移的模拟方法在文献中的报道较多，一般分为Euler-Euler法和Euler-Lagrange法^[65]。Euler-Euler法将流体和颗粒都视为连续相，Euler-Lagrange法则是将流体处理为连续相而将颗粒视为离散相，通过追踪每个颗粒的运动轨迹予以精细表征。在Euler-Lagrange法中，离散元方法（DEM）在研究颗粒流时考虑颗粒与流体之间的作用力，具备明显优势，被广泛应用^[66-68]。DEM的基本思想是将离散介质当作单元体，不同单元体之间通过角、面或边等发生接触，单元体的运动（包括平动和转动）受牛顿第二定律控制，单元体之间的接触力由合适的接触模型计算^[69]。线性接触模型是最基本的接触模型，在该模型中计算颗粒间接触力和力矩时均将其简化为正比于变形量，大幅简化了整个接触力及力矩变化的模拟。对于岩石颗粒系统，常用的接触模型还有线性平行键模型、平面节理模型和软键模型^[65]，这些模型简化了颗粒间的复杂非线性接触过程、降低了计算量，适合特定的工程问题。目前大部分的DEM模拟基于球形颗粒展开，而CO₂-水-岩系统产生的沉淀一般为非球形颗粒。相比于球形颗粒，非球形颗粒DEM模拟的难点主要体现在以下两方面^[70]：①非球形颗粒之间的接触判断更为复杂，并且当颗粒形状发生变化时，其碰撞接触算法也随之改变；②非球形颗粒之间的受力解析难以达到较高精度。目前对非球形颗粒形状的描述方法可分为多体组合方法和精确几何形状描述方法^[71]。多体组合方法采用球体或者其他形状的基本单元构造复杂颗粒形状，该方法的优势在于构造单元间的接触判断相对简单，需要协调几何形状的描述精度与模拟计算效率。精确几何形状描述方法是严格基于颗粒的几何形状进行颗粒碰撞判断，接触判断更加复杂，其优点在于颗粒接触判断精度高，其计算量与颗粒数目及接触判断算法的复杂程度相关，能够较好地拓展到存在尺寸差异的非球形颗粒体系模拟中。与溶解或沉淀类似，DEM也可以与传统CFD、粒子法（如SPH）、LBM等方法耦合模拟颗粒流^[72]，如图6所示。考虑到LBM模拟多孔介质流动的优势，这里仅对LBM-DEM流固耦合模拟的研究现状予以综述。在早期LBM研究中^[73]，为建立流体-颗粒的相互作用体系，当颗粒边界与流体网格相交时，固体颗粒的数值边界被认为处于LBM晶格的中点，但该方法无法准确耦合和计算流固相互作用力，可能在颗粒快速移动时导致流体曳力数值的振荡。Noble等^[74]提出了浸入式移动边界方法（IMB），认为固体颗粒上的拉格朗日点位于颗粒边界和流体网格的交叉点处，应采用流体网格单元被颗粒覆盖的面积分数（即固体体积分数）来计算流固相互作用力。在处理复杂移动边界问题时，近年来，越来越引人注意的方法是由Peskin^[75]提出的用于血液细胞流动模拟的浸入边界法（IBM），在该方法中一个颗粒由一组独立的点表示，拉格朗日点的位置可以简单地通过颗粒的位置和取向得出，颗粒的总作用力和扭矩通过加和颗粒上所有拉格朗日点的动量变化得到。该方法最初应用于球形颗粒的LBM-DEM流固耦合模拟，后来逐渐拓展到二维非球形颗粒^[76]和三维非球形颗粒^[77]的流动模拟中。

图6颗粒流动的DEM耦合方法（据文献[72]修改）

3 CO₂–水-岩系统的微观模拟

由于LBM在模拟多物理化学反应输运系统方面具有明显优势，近些年，已经有学者开始尝试采用LBM模拟CO₂-水-岩系统，并得到了一些重要的进展与认识。

对于溶解的模拟，通常采用溶解反应动力学模型预测溶解速率，孔隙尺度矿物溶解/沉淀速率可以用过渡态理论（TST）来描述。所有矿物的反应速率都可用如下一般形式表示^[78]：

根据（5）式，室温25 ℃条件下方解石矿物在CO₂酸性环境中的溶解/沉淀速率可以表示为：

（6）式中，k₁，k₂，k₃的值分别为8.9×10^-1，5.0×10^-4，6.5×10^-7 mol/（m²·s），K_eq,_CaCO3的值为3.47×10^-9[79-80]。除了上述计算方法以外，也有学者直接采用地球化学模拟器PhreeqcRM获取反应速率。不论是直接采用LBM模拟还是采用LBM与其他CFD离散方法结合（如FV、FE、FD等），再结合数字岩心技术与并行处理技术，基本可以在岩心尺度上实现对CO₂-水-岩系统溶解反应的三维孔喉空间模拟，定量评价溶解导致的孔喉变化、孔隙度变化及渗透率变化。在静态溶解模拟方面，Gao等^[81]将砂岩中的方解石胶结物考虑为可溶解、其他矿物考虑为不可溶解，利用LBM模拟了CO₂-水-岩静态溶解过程，还与前人^[82]开展的溶解实验进行了对比，如图7所示。模型实现了碳酸溶液与方解石长达1 440 d的反应过程（直到达到反应平衡状态）模拟，由于方解石被溶解，岩心孔隙度从1.1%增加到10.7%。他们发现方解石在岩心中的空间展布位置非常重要，由于部分方解石未与酸性流体直接接触，导致仍有超过三分之一的方解石未溶解。在动态溶解模拟方面，采用Patel等^[53]建立的LBM伪均质反应边界，Fazeli等^[83]模拟了碳酸盐岩岩样注入CO₂饱和盐水后的裂缝溶解演化过程，结果表明，流动路径上存在的非反应性矿物可以限制渗透率的增加，单矿物多孔介质中常用的孔隙度和渗透率立方关系在多矿物多孔介质中失效。An等^[84]为了模拟岩心尺度的CO₂-水-岩反应，将LBM与地球化学模拟器PhreeqcRM耦合，评价了不同注入速率、压力和温度下孔隙度和渗透率随时间的演变关系，随着流速的增加，岩石中的溶解作用更加均匀，认为高注入速率下的渗透率增幅要比低注入速率下的更大。Xie等^[85]利用LBM处理流体流动和组分输运，利用VOP法处理反应造成的边界演化，模拟CO₂埋存过程中碳酸钙的溶解过程，阐明压差对碳酸钙溶解的影响随温度的增加而减弱，根据温度和压力之间的关系绘制了碳酸钙溶解模式的图版，认为高温和高注入速率下注入CO₂有助于其实现溶解埋存。Ju等^[86]采用考虑了化学反应消耗与反应可逆影响的局部反应边界格式，模拟了CO₂地下咸水层封存中的对流混合过程，分析了流固溶解反应及H₂S杂质组分对系统稳定性的影响。

图7CO₂-水-岩静态溶解实验及对应的LBM模拟^[81]在矿物溶解的化学反应模拟中，重要矿物地球化学反应相关参数（平衡常数、反应速率常数等）的数值不确定，是限制数值模拟准确性的主要因素。例如，Gray等^[87]指出在利用LBM模拟CO₂饱和溶液时，需要考虑溶解反应的非线性动力学模型（（5）式中的n不为1）和多组分扩散系数模型，即便如此，在与Menke等^[88]碳酸盐岩实验结果对比时，发现模拟结果快几倍（见图8），认为是模型一些关键参数的不确定性引起的。此外，实验中溶解速率与时间呈平方根关系，而模拟中溶解速率与时间近似线性相关，认为可能是CT扫描精度无法捕捉一些更小孔隙中的离子扩散作用，导致模拟中的溶解速率更高。在孔隙尺度的反应模拟中，往往对一些反应速率较慢的过程考虑为不反应，可能低估部分矿物的反应速率。而在长时间、宏观尺度的模拟中，往往对一些快速反应的矿物（如CaCO₃）处理为平衡相，这又可能高估实际多孔介质中的反应速率。

图8CO₂-水-岩动态驱替-溶解实验及对应的LBM模拟^[87]在沉淀生成及运移模拟方面，除了前文提到的无机盐结晶沉淀模拟研究以外，还有如下学者开展了CO₂-水-岩系统的沉淀模拟。Jiang等^[57]采用颜色LBM求解两相流动方程、FV求解组分输运方程，考虑无机盐的沉淀出现在浓度达到一定阈值时，模拟发现沉淀引起的孔隙结构变化显著降低了多孔介质的绝对渗透率，且对非润湿相相对渗透率的影响大于对润湿相相对渗透率的影响。Liu等^[63]利用LBM求解流动方程、FV求解组分输运方程，通过Stokes定律的简化方法，将受溶解释放的颗粒迁移耦合至LBM反应输运模型，探索性地模拟了固体沉淀的运移过程，如图9所示。通过与在线X射线断层扫描仪（Micro-CT）监测的动态岩心驱替实验结果对比，发现模拟结果可以很好地匹配不同时间步下沿驱替方向的孔隙度变化。通过对碳酸盐岩多孔介质的大量模拟，发现了渗透率和孔隙度之间的幂律关系，同时也探讨了由于溶解释放而引起的固体颗粒迁移对渗透率的影响。Parvan等开展了相似的模拟研究^[89]，但是在模拟中假设所有固体与流体之间的表面都可以进行反应，忽略了岩石表面的化学非均质特征以及有效的反应接触面。除了针对单一矿物沉淀的研究存在不确定性以外，还有学者指出，溶解反应后生成的部分次生矿物的种类难以确定、生成量难以量化、反应动力学参数缺失，导致部分次生矿物的模拟难度大，也可能对化学反应模拟结果产生较大影响^[90]。

图9CO₂-水-岩动态驱替沉淀-运移实验及对应的LBM模拟^[63]

4 存在的问题及建议

①CO₂-水-岩石系统中复杂作用机理（溶解、沉淀与沉淀运移等）的解耦合。CO₂注入储集层打破地层水与岩石间的物理化学平衡，产生溶解、沉淀和沉淀运移等现象，它们的综合作用导致储集层渗透率可能出现增大、减少、波动变化等复杂结果。溶解—沉淀—沉淀运移是一种孔隙尺度行为，传统岩心渗流实验往往难以直观刻画其作用过程（如离子浓度变化、孔隙形貌演化、沉淀运移追踪等），导致对储集层物性参数变化规律认识不清、预测难度大。因此，还需要进一步深入剖析CO₂-水-岩石系统中溶解、沉淀与沉淀运移的耦合渗流作用机理。具体而言，可以通过不同实验设备及可视化手段对各种复杂机理（溶解、沉淀和沉淀运移）解耦合，系统研究各个机理的作用特征及规律。例如，借助静态溶解实验，分析多矿物差异溶解规律；借助刻蚀玻璃实验，捕捉地层水过饱和沉淀（结晶）及沉淀运移规律；借助岩心在线驱替实验，分析驱替条件下的溶解、沉淀及沉淀运移规律。②多矿物差异性反应的反应系数测量及微观表征方法。总体而言，相对于传统CFD或拉格朗日描述的粒子法等方法，LBM在模拟多孔介质内的复杂流动方面具有明显优势，包括易处理复杂边界、易捕捉相界面及易并行等。针对反应输运流，可以通过在LBM方程中添加源汇项并结合恰当的边界更新法，实现溶液中均相反应和流固接触面非均相反应的表征，已在各种物理化学反应流中得到了广泛应用。在CO₂-水-岩系统中，CO₂水溶液与不同地层水组分、不同矿物之间可能发生反应，且反应速率存在差异。需要针对性地开展不同温度、压力环境下，CO₂水溶液与离子组分、矿物组分之间的反应速率标定实验，提高实验测试的准确性，降低主要矿物地球化学反应相关参数（平衡常数、反应速率常数等）的数值不确定性，最终增加表征组分交换流动方程的数值预测可靠性。③晶体成核和矿物脱落作用下的沉淀生成机理及表征。沉淀的生成位置、大小及运移能力直接影响多孔介质的渗流能力，因此，沉淀生成表征是CO₂-水-岩系统微观模拟成功与否的关键。沉淀的生成包含无机盐的析出和非反应矿物的脱落，其中，无机盐的结晶成核受概率成核理论控制，且晶体生长也遵循一定的规律，但目前的相关研究报道对结晶沉淀的表征过于简化，将结晶成核和晶体生长简单处理为一个随机过程。受周边矿物溶解和胶结薄弱处破坏的诱导，非反应矿物从主固相体脱落是生成沉淀的另一个重要机理，但目前尚未发现相关机理及表征方法的研究。因此，需要基于晶体成核理论预测无机盐的生成过程、提出非反应矿物脱落的判别方法，从而建立可靠的沉淀生成表征方法。④不规则、多类型、变数量、变形状沉淀条件下的沉淀-流体相互作用模拟方法。沉淀（颗粒）在多孔介质中的运移是影响CO₂-水-岩系统渗透率的关键因素。DEM是多孔介质颗粒流模拟最常用的处理方法，结合IBM法，可以实现LBM体系下的颗粒-流体相互作用模拟。CO₂-水-岩系统的颗粒运移模拟主要存在3个方面的特殊性：沉淀的产生方式多样（次生矿物、矿物颗粒脱落）；沉淀的形状一般不规则；生成的沉淀在运移过程中可能消失（溶解），也即沉淀的数量和大小在运移过程中可能变化。因此，需要针对不同的沉淀生成机理和沉淀的几何特征，建立不规则、多类型、变数量、变形状沉淀条件下的沉淀-流体相互作用模拟方法。⑤CO₂-水-岩系统中多物理化学过程的耦合渗流机制。近年来，已有学者采用LBM模拟CO₂-水-岩系统的反应流，但建立的模拟方法主要针对某一个或某几个物理化学过程，然而CO₂注入储集层中是一个多相（气/水/油）、多组分（CO₂水溶液，H⁺，HCO₃^–，CO₃^2-和OH^–等）、非均相/均相化学反应的过程，对这些复杂过程实现全耦合模拟是一个具有挑战性的问题。CO₂-水-岩系统中的各个物理化学过程之间相互耦合、相互影响，例如，矿物的溶解可能诱导周边矿物脱落产生沉淀，沉淀运移后暴露的新鲜矿物面接触CO₂后可能产生新的溶解。根据文献调研结果，目前针对该问题的研究主要针对上述某一个或某几个物理化学过程，若未有效协调各机理之间的耦合关系，将导致模拟结果与室内实验结果存在较大差异。因此，需要有机衔接各个物理化学过程之间的模拟边界条件、协调各机理的耦合关系，从而定量评价溶解、沉淀、沉淀运移对孔隙度、渗透率影响。

5 结论CO₂注入地下后，打破地层水与岩石之间的物理化学平衡，产生溶解、沉淀及沉淀运移等特殊现象，影响CO₂注入性、盖层封堵性及埋存可靠性。目前基于宏观高温高压静态实验、岩心渗流实验等方法，难以从本质上深入认识CO₂-水-岩系统的溶解—沉淀—沉淀运移耦合作用机制，无法有效预测宏观孔渗参数的复杂变化。开展这些特殊物理化学过程的耦合渗流机制及微观模拟研究，是量化物理化学反应导致的多孔介质孔喉结构变化、预测复杂渗透率变化规律的有效手段。CO₂-水-岩系统溶解—沉淀—沉淀运移的耦合渗流机制及微观模拟的研究，尚需要解决一些关键的科学问题。例如，考虑CO₂溶于地层水与不同矿物之间发生的差异溶解反应，建立非均相矿物溶解的模拟方法；针对不同类型的沉淀生成机理（如化学反应产生的次生矿物、溶解诱发的矿物颗粒脱落），建立合理的沉淀生成表征方法；建立不规则、多类型、变数量、变形状沉淀条件下的沉淀-流体相互作用模拟方法；有机衔接各个物理化学过程之间的模拟边界条件，协调各机理的耦合关系。合理解决这些问题，可以实现碳埋存多孔介质内溶解—沉淀—沉淀运移的可视化模拟。在未来研究中，还需要深入结合室内可视化实验与微观模拟，通过静态溶解实验、刻蚀玻璃实验及岩心在线驱替实验对溶解—沉淀—沉淀运移解耦合，为各过程的微观模拟方法提供验证数据。例如，以LBM和DEM为基本模拟手段处理流动、传质和沉淀运移过程，在其中耦合基于溶解动力学的非均相矿物溶解模拟方法、基于晶体成核理论和矿物脱落判别法的沉淀生成模拟方法，形成溶解—沉淀—沉淀运移耦合渗流的微观模拟技术，实现CO₂-水-岩系统渗流相关参数的定量预测，为CCUS和CCS工业活动的安全、高效实施提供理论依据与技术支撑。符号注释：

参考文献：

———END———

中国石油大学（北京）联合国家自然科学基金委中国21世纪议程管理中心、国际二氧化碳捕集利用封存技术创新合作组织（筹）、联合国教科文组织“碳中和与气候变化驱动绿色转型”教席、中国岩石力学与工程学会以及中关村绿色矿山产业联盟，定于2025年5月28-29日在北京市召开“2025第三届二氧化碳捕集、利用与封存（CCUS）国际会议”。

为给论文作者提供充足的准备时间，鼓励广大作者撰写高质量论文。会议将开启第二轮征稿，论文摘要截止时间延长至2025年 5月18日。

同时， 为满足参会者在学术交流方面的多样化需求，大会组委会在征集英文稿件（Scopus收录的CCUS 会议论文集、Applied Energy、Engineering、Journal of Energy Resources Technology 的优先评审）的基础上新增中文稿件征集。中文稿件可优先推送至《钻采工艺》的”油气藏开发提高采收率-CCUS 技术”专刊及《石油勘探与开发》、《中国石油勘探》。论文摘要截止时间：2025年5月18日；论文全文截止时间：2025年7月20日;

论文投递方式：

将论文摘要或全文

（注：论文摘要与参会回执表需同时投递）

发到邮箱paper@iccusc.com

联系电话：马老师15822801766

（扫码添加会议联系人企业微信咨询）扫码下载会议通知

会议咨询丨马老师 联系电话15822801766（微信同）行研及合作咨询丨李老师 联系电话13581630352（微信同）