CO2 地质封存风险分析的多场耦合数值模拟技术综述1)

于恩毅 邸 元,2) 吴 辉 曹小朋 张庆福 张传宝

* (北京大学工学院,北京 100871)

† (北京大学地球与空间科学学院,北京 100871)

** (中国石油化工股份有限公司胜利油田分公司勘探开发研究院,山东东营 257015)

引言

二氧化碳(CO2)被认为是主要的温室气体之一[1].由于大量排放,CO2在大气中的含量不断攀升,引起全球气候变化,严重威胁人类的生存环境,减少CO2的排放已成为当今世界关注的焦点问题[2-3].在可预见的未来,化石能源仍将在能源消费体系中扮演重要角色,CO2捕集、利用和封存(CO2capture,utilization and storage,CCUS)是实现化石能源低碳化利用、减少CO2排放最直接和最关键的技术途径[4],能够为实现碳中和目标提供重要支撑[5].

CO2地质封存是指通过工程技术手段将捕集的CO2注入深部地质储层,实现CO2与大气长期隔绝的过程[5].如图1 所示,封存点一般选择在深度800～1000 m 以下,渗透性好且上方存在低渗盖层的储层[6].储层类型主要包括咸水层、油气藏和煤层等.在所有封存类型中,深部咸水层分布广泛,封存容量占比约为98%[5].油气藏也是较为理想的CO2封存场所,因为存在详细的地质勘探资料和较完备的地面设施等基础条件[5],并且将CO2注入油气藏能够利用CO2驱替达到提高原油采收率的目的[7-9].

图1 CO2 地质封存的途径[10]Fig.1 Approaches of CO2 geological storage[10]

CO2在储层岩体中的封存过程较为复杂,受到地层特性和地球化学反应等影响[11],涉及温度场-渗流场-力学场-化学场(T-H-M-C)的耦合作用[12],如图2 所示.CO2注入地层后会导致较大的流体压力积聚,导致有效应力场发生变化,应力场的改变会影响岩体孔隙度、渗透率和毛管压力等,从而影响CO2的注入和流动[13-14].注入的CO2为超临界态,温度要明显低于周围地层温度,二者的温度差导致局部区域内的温度场发生变化,从而改变CO2流体的密度、黏度和溶解度等,影响其渗流特性,温度变化引起的热应力也会直接改变岩体的应力状态[13].CO2易溶于水,进而与周围岩石矿物发生化学反应,溶解岩石或钙化沉淀,并可能与盖层有机质发生反应,导致盖层渗透率和孔隙度等性质发生变化[11,15].在多孔介质中的化学反应速率又受到温度、压力、渗流速度和CO2扩散速率等因素的影响[16].

图2 CO2 地质封存多场耦合原理[13]Fig.2 Multi-field coupled principle of CO2 geological storage[13]

CO2地质封存的主要机制包括: (1)构造封存,即低渗的密封盖层限制CO2的迁移[17-18];(2)毛细封存,即CO2被孔隙结构中的毛管力束缚[19];(3)吸附封存,即CO2吸附在黏土矿物表面[20];(4)溶解封存,即CO2与地层水或油混溶[21-22];(5) 矿化封存,即CO2与岩石发生化学反应生成碳酸盐矿化物[23-24].当CO2注入到地层后,首先圈闭在油气藏或深层盐水层中,CO2的运移由多相对流过程主导,主要由注入压力和流体密度差异引起的浮力驱动[25].随着时间的推移,保持自由相的CO2通过岩石中的渗流通道向上迁移,剩余的CO2移动非常缓慢,直到最终被残留在孔隙中、溶解在地层水中,或作为碳酸盐岩矿物沉淀[10,26].

CO2地质封存有其自身的风险.若CO2垂向运移至近地表区域,会造成浅层饮用水污染[27-30]、逃逸至大气环境[31]等风险.所以,在CO2注入作业前,必须对储层-盖层系统进行力学分析和稳定性评估[32-33].盖层一般是未受干扰、低渗透、厚且横向分布广泛的地层,常见为页岩、泥岩和碳酸盐岩等[34],具有高毛细管力和突破压力,以防止注入流体泄漏[35].在CO2封存过程中,盖层必须承受短期的过量注入压力和长期的浮力驱动压力[36].断层和裂缝带是盖层岩体的软弱带,在CO2注入过程中,盖层中的裂缝可能发生起裂、扩展,断层可能被活化,发生滑移[37-38].CO2注入还可能会导致储层中的岩体发生压裂,裂缝可向上延伸至盖层,甚至穿透盖层[39].天然裂缝/断层和压裂裂缝形成了CO2逸出的潜在通道,对盖层完整性构成威胁[35,40-41].

储层-盖层系统多场耦合数值模拟是分析盖层完整性和断层稳定性并进行预警控制的关键核心技术[25,42].CO2地质封存过程的数值模拟需要考虑大时间尺度和大空间尺度下的多物理过程、多场耦合作用、储层非均质性和流体相态变化等因素[43-45],需要研发计算效率高、计算能力强的CO2地质封存数值模拟程序.本文对CO2地质封存储层-盖层系统多场耦合数值模拟方面的研究进行了全面综述,介绍了CO2地质封存过程中T-H-M-C 多场耦合问题的数学模型、耦合问题数值解法和封存风险分析等方面国内外的研究进展,对当前我国CO2地质封存数值模拟技术面临的主要问题进行了总结.

1 数学模型

CO2地质封存多场耦合问题的数值模拟涉及的研究对象是岩石固体骨架与孔隙流体(CO2、水、油等)组成的多相多组分系统,其数学模型主要包括基本控制方程和相平衡计算方法.基本控制方程包括应力平衡方程、质量守恒方程、能量守恒方程和化学反应等.

1.1 应力平衡方程

根据可变形多孔介质多相流孔隙热弹性理论,应力平衡方程为[12,46]

式中,ω 和 λ 为拉梅常数,u为位移,P为流体压力,Ks为岩石骨架体积模量,Km为基质颗粒体积模量,β为骨架的体积热膨胀系数,T为温度.本文公式中的变量均为国际单位制.

用平均应力表述,应力平衡方程也可写为[47-48]

式中,v为泊松比,σm为平均应力,α 为Biot 系数.

1.2 质量守恒方程

采用组分模型对CO2地质封存进行数值模拟,α组分的质量守恒方程为

对于 α 组分,其质量累积项和传输项的表达式分别为

式中,下标 β 表示相(包含油相、气(CO2) 相和水相),ϕ 为孔隙度,ρ 为密度,S为饱和度,表示 α 组分在 β 相中的质量分数,Fadv为对流项,Fdif为扩散项,V为质量流速,j为扩散速度.

岩石基质中流体的流动符合达西定律[49],其表达式为

式中,k为渗透率,kr为相对渗透率,µ 为黏度,g为重力加速度,D为深度.

考虑CO2扩散,扩散项可采用 Fick 定律描述

式中,D为扩散系数.

对于低速非达西渗流和高速非达西渗流,达西定律描述的渗流速度与水力梯度之间的线性关系不再成立[50].对于低渗透-致密储层中的低速非达西渗流,只有驱动压力梯度 ∇P-ρg∇D大于拟启动压力梯度L时,流体才能流动,其表达式为[51]

对于裂缝中的高速非达西渗流,可采用Forchheimer定律描述,其表达式为

式中,χ为高速非达西渗流系数.

1.3 能量守恒方程

能量守恒方程与质量守恒方程具有相同的形式

式中,ρR为岩石密度,CR为岩石比热,u为比内能.

忽略辐射传热,热量的流量由热传导项与热对流项组成,其表达式为

式中,κ 为热传导系数,h为比焓.

1.4 化学反应

CO2注入地层后发生的化学反应主要包括CO2与地层流体之间的离子交换和溶解的CO2与地层岩石之间的矿物反应[21].对于T-H-M-C 多场耦合问题,化学反应过程通过质量和能量守恒方程中的化学反应源项R来表述,其表达式为[12,26]

式中,ξ 为化学反应计量系数,r为化学反应速率.矿物溶解和矿化反应的化学反应速率计算公式为[12]

式中,A为反应接触面积,kR为速率常数,Q为化学亲和性,Keq为平衡常数.化学反应速率r若为正,表示矿化沉淀;若为负,表示矿物溶解.

根据不同的反应机理,矿物溶解和沉淀不仅受到纯H2O 的催化,而且受到H+和OH-的影响.因此,速率常数kR的计算公式包含3 项[47]

上述平衡常数和速率常数对化学反应动力学模型的准确性至关重要,然而由于地球化学过程的复杂性,特别是相关机制的持续时间极长,准确的测定较为困难[26].

1.5 岩石应力-应变本构关系

采用双孔隙弹性理论描述岩石的应力-应变关系[52]

式中,G为剪切模量,为裂缝(i=1)和基质(i=2)中的流体和热效应耦合系数,δij为Kronecker符号.

基于等效连续介质模型,裂缝性岩石的本构关系为

式中,ε 为正应变,γ 为剪应变,σ′为有效应力,τ 为剪应力,GM和GF分别为基质和裂缝的剪切模量,EM和EF分别为基质和裂缝的杨氏模量,ECF为裂缝的压缩模量,αM和 αF分别为岩体中基质和裂缝的体积分数.

由于裂缝性岩石中基质占据的体积远大于裂缝占据的体积,因此可认为 αM≈1,式(17)中的裂缝性岩石柔度矩阵C可以写为基质与裂缝柔度矩阵的加权平均

式中,CM为基质柔度矩阵,CF为裂缝柔度矩阵.

裂缝的力学性质与岩石基质存在很大差别,其力学参数和渗流能力会随施加在裂缝壁面上的法向载荷产生动态变化,需要通过裂缝柔度矩阵获得裂缝壁面位移与岩石应力间的关系,以对裂缝孔隙度、渗透率等参数进行更新[53].对于天然裂缝在法向加载下的变形过程,Goodman 等[54]指出裂缝在法向受压下的变形本构关系近似为双曲线型模型.Bandis等[55]对大量存在天然裂缝的岩块进行实验,总结出裂缝受压变形的经验性双曲线型本构公式.

1.6 相平衡计算及辅助方程

地层条件下油-水-气(CO2)系统的相平衡计算是油气藏地质封存CO2数值模拟研究的核心问题之一.油-水-气(CO2)的相平衡计算主要有平衡常数法、基于Rachford-Rice 函数的闪蒸计算法和Gibbs自由能最小化法3 种计算方法[56].

平衡常数法采用基于实验的经验公式来计算平衡常数,计算速度快、稳定性良好,虽然避免了复杂、耗时的相平衡计算,但由于缺乏严密的热力学理论基础,计算结果往往偏差较大.基于Rachford-Rice 函数的闪蒸计算法主要用于油-气(CO2)两相体系的相平衡计算,同平衡常数法相比较,有较高的精度,但计算量较平衡常数法大.适用于油-水-气(CO2) 三相体系的闪蒸计算法还有待研究和发展.Gibbs 自由能最小化方法理论严密、适应性广、计算稳定和收敛性好,比较适合于油-水-气(CO2)三相的相平衡计算,但是因为需要反复迭代,计算量比平衡常数法和闪蒸计算法都要大很多.

由于渗流场、温度场、应力场和化学场之间存在耦合作用,在求解T-H-M-C 多场耦合问题的过程中,需要实时更新计算的物性参数,主要包括毛细管力、渗透率、相对渗透率、孔隙度、流体黏度和密度等.相对渗透率曲线可以采用Brooks-Corey 模型.考虑应力场、温度场和化学场对孔隙度的影响[12,53]

式中,ϕ0为0 应力下的初始孔隙度大小,σm为平均总应力(压为正),Kb为岩体的体积模量,βT为热膨胀系数,T为当前时刻温度,cs为固相s组分浓度,下标 0 表示初始状态.

基质渗透率和毛细管力可根据孔隙度变化进行更新

式中,k0和Pc0分别为0 应力下的初始渗透率和初始毛细管力,d为实验测定的参数.

温度变化会造成流体黏度的变化,一般来说随着温度的升高流体黏度降低,可采用经验公式描述[53]

式中,A,B,C和D均为实验测定参数.

2 耦合问题数值解法

对于CO2地质封存的T-H-M-C 多场耦合问题,可采用有限元和有限体积混合的方法进行离散,即对于式(1)的应力平衡方程采用有限元法离散,对于式(3)质量守恒方程和式(10)能量守恒方程采用有限体积法进行离散.离散后得到的方程可采用全耦合、迭代耦合、弱耦合、显式耦合和拟耦合等数值方法进行求解,如图3 所示.

图3 常用的耦合问题数值解法Fig.3 Commonly used coupling solutions

全耦合方法即为全隐式的数值求解方法,该方法对离散后的应力守恒方程、质量守恒方程、能量守恒方程和化学反应方程同时求解,所有参数变量需要在每个迭代步内进行更新[13].全耦合方法的收敛性最好,计算精度最高,但由于耦合问题的方程数量庞大,全耦合方法的计算效率往往较低.

为了提高多场耦合问题数值模拟的计算速度,通常采用半隐式的方法进行求解,即优先对一部分控制方程(一般是质量守恒方程和能量守恒方程)进行隐式求解,之后将求解得到的未知量作为已知量代入另一部分尚未求解的控制方程中,求解其余的未知量.显式耦合方法、弱耦合方法和迭代耦合方法均属于半隐式求解方法.

迭代耦合方法中,在每个时间步的牛顿迭代内首先对质量和能量守恒方程进行隐式求解,将计算得到的平均孔隙压力传递至应力平衡方程中以求解岩石的变形情况,根据岩石变形程度对孔隙度、渗透率等关键物性参数进行更新,开始下一迭代循环,分别对质量守恒方程、能量守恒方程和应力平衡方程进行计算,直至达到收敛条件进入下一时间步.与全耦合方法相比,迭代耦合方法显著提高计算效率[57].

在显式耦合中,渗流场、温度场和应力场的计算独立进行,每个时间步中信息仅从渗流场、温度场向应力场单向传递一次,即在每一时间步中仅进行一次应力平衡方程计算和参数更新,而应力场的信息则不会反馈给渗流场、温度场.显示耦合方法是耦合程度较低的方法,计算精度有所降低,但计算效率高、收敛性好,适用于大尺度、复杂地质条件和复杂岩体本构模型的力学模拟[13].

弱耦合方法是在显式耦合方法的基础上进一步减少应力平衡方程的计算频率,在每个时间步内对渗温场和应力场进行独立计算,每隔数个时间步才进行一次应力平衡方程的计算和参数更新,如此反复循环迭代直至结束.该方法计算结果精度有所降低,但计算速度快、收敛性好、适用范围广,能够耦合不同用途的程序[13].

拟耦合方法并不进行应力场求解,仅通过经验公式或数据表建立孔隙度、渗透率等与孔隙压力的关系.该方法无法对岩体的变形过程进行模拟,因此精度有限.

3 CO2 地质封存风险分析

分析CO2长久封存过程的安全性是多场耦合数值模拟的关键应用.图4 展示的是在CO2地质封存中储层-盖层可能发生的典型地质风险,主要包括CO2泄漏、地表变形、盖层破坏和地震诱发等[39,58-61].盖层破坏的常见形式包括断层活化和盖层压裂(裂缝扩展),断层活化主要与断层剪切破坏有关,盖层压裂(裂缝扩展)主要与盖层拉伸破坏有关[58].

图4 CO2 封存储层-盖层中的地质风险示意图[39]Fig.4 Geomechanical risks during CO2 geological storage[39]

3.1 CO2 泄漏

CO2通过盖层泄漏的主要机理包括压力驱动的流体对流[62]和浓度驱动的分子扩散[63],如图5 所示.断层/裂缝和岩石基质存在巨大的渗透率差异[64],断层核心可视为阻碍流动的屏障,而断层两侧破碎的损伤区可视为沿断层的流动通道[35,65].在CO2的注入过程中,流动路径主要受裂缝网络几何形状控制,CO2在遇到裂缝时优先沿裂缝方向运移,在裂缝周围扩散并集中分布在裂缝周围,如图6 所示.随着CO2的积累,碳酸盐岩矿物在裂缝交叉点沉淀,降低了裂缝网络的整体连通性,低渗透带内的物理圈闭效应会显著增强[66].裂缝组合导致储层渗透率各向异性,相对渗透率和毛管压力效应也会影响储层各向异性,从而显著影响CO2的分布[67-69].

图5 通过断层泄漏到地表的潜在途径[70]Fig.5 Potential CO2 leakage path through a fault to ground surface[70]

图6 裂缝对CO2 运移的影响(CO2 饱和度分布)[67]Fig.6 Effect of fractures on CO2 migration (CO2 saturation distribution)[67]

CO2在多孔介质中的扩散速度相对于对流速度较慢.CO2的扩散速度对流体-岩石反应较为敏感,流体-岩石反应可能提高(矿物溶解孔隙度增加)或降低(矿化反应消耗CO2)扩散速度[70].但是现有研究表明,扩散并非深层储层CO2的主要泄漏机制.根据Wang 等[71]的数值模拟结果,将CO2以27 MPa的注入压力注入储层(储层压力为19 MPa),如图7所示,稳定注入1000 年后CO2扩散距离仅为0.93 m.此外,根据Busch 等[70]的估算,纯水通过扩散透过10 m 厚度的薄盖层(有效扩散系数为2.0×10-9m2/s),大约需要105年.

图7 CO2 扩散运移数值模拟结果[71]Fig.7 Numerical simulation of CO2 diffusion[71]

3.2 地表变形

阿尔及利亚In Salah 二氧化碳地质封存项目已经观测到由CO2注入导致的地表抬升现象[72-74],如图8 所示.地表变形的分析方法主要包括解析法、半解析法和数值模拟方法[33].Fjaer 等[75]针对简化的一维横向薄储层提出计算储层垂直膨胀的解析式,Xu 等[76]针对半无限空间的问题推导了压力源地面变形的解析解,Selvadurai[77]提出地表岩层和深部岩体之间相互作用的基本模型,并给出地表岩层隆起的半解析式.对于实际CO2地质封存问题中复杂岩体的变形,需要采用数值模拟方法进行计算.

图8 In Salah KB-502 井地面垂直位移的演化[72]Fig.8 Transient evolution of vertical ground displacement at the In Salah KB-502 well[72]

Rinaldi 等[72]通过数值模拟分析了KB-502 地表位移演化情况,并与监测数据进行了比对,如图9 所示.郝术仁等[78]模拟研究了不同流量和不同压力注入CO2引起地表位移的变化规律.Siriwardane 等[45]模拟分析了盖层断层/裂缝带对地表位移的影响,无断层/裂缝时地表变形较小且对称,有断层/裂缝时地表位移较大且非对称.裂缝带渗透率较低时,地面垂直位移较小;盖层裂缝带离注入源越近,地表变形幅度越大.

图9 In Salah KB-502 井2006 年12 月23 日地表变形模拟结果[72]Fig.9 The simulation results of ground displacement at the In Salah KB-502 well on December 23,2006[72]

3.3 剪切破坏及断层活化

流体压力的演化直接影响盖层的地质力学稳定性,盖层的低渗透性阻止了超压向盖层快速传播,但在靠近储层界面的盖层部分,由于流体压力的积聚导致有效应力降低,应力状态向失效状态演化[79-80],如图10 所示.根据Mohr-Coulomb 强度失效准则(如下式),当作用在断层面上的剪应力超过断层的剪切强度时,断层活化[81-82]

图10 Mohr-Coulomb 强度失效准则[59]Fig.10 Mohr-Coulomb failure criterion[59]

表1 国外部分CO2 封存储层的力学参数[59]Table 1 Stress state at several CO2 injection sites abroad[59]

CO2注入地层后流体压力积聚、温度场的变化,会改变有效应力场,化学反应会降低岩石的强度[84].孔隙压力的增加会导致莫尔圆半径减小并左移(参见图11 中的红色莫尔圆).CO2到达注入井底部时温度低于储层温度,注入井周围的岩石发生冷却,温度变化诱发温度应力[59],应力状态向剪切破坏状态转移(参见图11 中的蓝色莫尔圆).温度应力主要取决于温差、岩石基质的体积模量和体积热膨胀系数[85],由温度变化 ΔT引起的孔隙压力变化 ΔP可以根据下式进行估算[86]

图11 应力状态耦合效应示意图[59]Fig.11 Coupled effects on stress state[59]

式中,Δ ζ 为流体应变,βf为流体体积热膨胀系数,为岩石基质体积热膨胀系数,Kf为流体体积模量,Kb为岩石体积模量.

断层除非位于注入井附近,一般不会直接受到温度差的影响.但在注入后期,储层岩石冷却收缩,水平应力减小,从而引起远场应力的变化,导致远离注入井的断层可能处于不稳定状态(参见图11 中的绿色莫尔圆)[87].

3.4 拉伸破坏及裂缝扩展

盖层由于CO2注入发生拉伸破坏时,裂缝起裂和扩展方向依赖于最小主应力方向和储层非均质性[88].当孔隙压力P超过盖层的最小主应力 σ3时,发生拉伸破坏,裂缝开始萌生并扩展[89-90]

当储层或盖层岩石剪应力超过剪切强度时,储层或盖层已有的裂缝会发生剪切滑移[91],如图10 所示.如果储层内部形成压裂裂缝、裂缝不向盖层扩展,或者裂缝仅在盖层有限的范围内延伸时,则不会发生泄漏,且能够提高CO2的注入能力[90-93].

4 数值模拟方法研究需解决的主要问题

在地下能源和资源领域,国际上已经研发出了诸多用于多孔介质多相流数值模拟的程序和软件,能够针对CO2地质封存进行T-H-M-C 多场耦合或部分多物理场耦合的模拟计算,表2 给出了目前国外主要的可用于CO2地质封存的数值模拟器[25-26,32].针对CO2地质封存,近期国际石油工程师协会发布了第11 个标准算例(the 11th SPE CSP)[94].该算例包括一个实验室尺度的二维模型、一个油藏尺度的二维和一个油藏尺度的三维地质模型,可用于对不同的数值模拟方法和模拟器进行对比研究.由于CO2地质封存需要对目标封存场地储层-盖层体系的三维地质结构进行精细化建模,需要考虑大时间尺度和大空间尺度下的多场耦合作用,所以在计算 效率和计算能力方面,对数值模拟技术提出了更高的要求.除CO2地质封存过程中复杂的相态计算和地球化学反应涉及的本构关系等问题以外,在数值模拟方法研究方面,还需要解决如下问题.

表2 国外部分可用于CO2 地质封存的数值模拟器Table 2 Overview of numerical simulators for CO2 geological storage

4.1 基于角点网格的有限单元法

角点网格可以准确描述地层界面、断层面和尖灭等复杂地质情况以及储层的空间展布和属性分布特征,因此通常都采用角点网格进行三维地质体精细化建模[95-96].角点网格系统中,通常存在重复节点,其节点无法传递力、位移和温度等连续性变量,相邻网格之间允许错动、断开和夹层,因此基于角点网格建立的地质模型无法直接用于三维地应力的有限元模拟.基于角点网格数据进行有限元网格的剖分,几何处理上较为困难,且有限元网格数量十分巨大.因此,需要结合角点网格和有限元方法的优点,建立基于角点网格的地应力有限元模拟方法,这样既可以完整准确地保留地质模型复杂的几何特征,也可以利用有限单元准确地计算三维应力场[97].

4.2 多场耦合的复合介质模型

CO2注入涉及的地层区域往往包含多尺度的裂缝系统,从长度米级的小裂缝到公里级的大断层[33].对于多尺度裂缝系统,仅采用多重介质模型,会丧失断层或大裂缝准确的几何信息;若仅采用离散裂缝模型,计算量极大,不适用于实际封存问题的模拟;采用嵌入式离散裂缝模型,连续介质虽然能够采用结构化网格,但是需要细分裂缝,并对所有不同尺度的裂缝进行独立表征,将导致网格数量急剧增加,造成极大的计算负担,且不适用于各向异性的连续介质.

为了更准确表征裂缝型储层和提高数值模拟计算效率,应分别对不同尺度的裂缝进行建模: 小型裂缝缝长远小于网格尺度,对于渗流特征的影响局限在基质网格内部,可采用等效单重介质模型建模;中等裂缝缝长与网格尺度相当,集中分布常以裂缝簇的形式赋存,其方向性与连通性对集中分布区域的渗流特征造成较大影响,可采用双重介质模型建模;大型裂缝和断层数量有限,但对于整体渗流场的分布具有十分重要的影响,为了准确描述每一条大型裂缝或断层的水力特征,应采用整体嵌入式离散裂缝模型建模[98].由此形成复合介质模型,将离散裂缝模型、多重介质模型和整体嵌入式离散裂缝模型结合在一起,进而在其基础上建立相应的T-H-M-C 多场耦合复合介质数学模型.

4.3 分区域的有限元-有限体积混合数值离散方法

CO2地质封存数值模拟需要考虑储层的上覆盖层、周边地层和下部地层,建立基层-储层-盖层-地表全地层模型[87,99].CO2地质封存场地尺度的地质建模,往往达到百万级以上的网格数量[100-101],模拟计算的规模十分庞大.图12 是胜利油田G89-1 区块的CO2地质封存全地层模型,其网格数量约240 万,待解的自由度数已达约1500 万.

图12 胜利油田G89-1 区块全地层模型Fig.12 The full synthetic field model of G89-1 block in Shengli oilfield

CO2地质封存的应力计算区域通常远大于渗流和温度计算区域.CO2地质封存数值模拟中对于渗流和温度计算具有较高的要求,储层内网格要保证一定的细密程度.而全地层模型应力场的计算则并不需要采用相同细密程度的网格划分,由于应力平衡方程采用有限元法计算,其主变量个数也远大于渗流场和温度场的计算,采用粗网格可以避免运算资源的浪费[54].因此,对于全地层模型的多场耦合问题,在渗流和传热计算区域内,采用细密网格进行剖分(如图13 所示),基于有限体积法对质量守恒方程和能量守恒方程进行数值离散;在应力计算区域内,采用粗网格进行剖分,基于有限元法对应力平衡方程进行数值离散.根据粗细网格间的映射关系,将离散方程联立求解,从而形成分区域的有限元-有限体积混合数值离散方法.这一方法能够在保证储层渗流场和温度场计算精度的同时,降低应力场求解带来的计算负担,有效提高计算效率.

图13 分区域网格划分示意图[54]Fig.13 Sub-regional grid division diagram[54]

4.4 高效的并行求解技术

CO2地质封存问题涉及的时间尺度和空间尺度较大.在多场耦合作用的数值模拟方面,需要发展高效的并行求解技术.采用并行计算技术,可将一个复杂的大型计算问题分成若干相互独立可以同时计算的子任务[102-103].常见的分解方法大体分为时间和空间上的并行,时间并行方法常采用流水线的方式实现,空间上并行是目前绝大多数计算问题的并行解决方案,目前常采用区域分解算法[102].

对图12 的全地层模型注CO2问题,采用并行技术进行多场耦合的数值模拟,计算得到注入3.2 年时的竖向地应力如图14 所示.这一计算虽然采用了并行技术(Intel Xeon Platinum CPU 48 核),但模拟耗时仍达40 h.

图14 胜利油田G89-1 区块储层竖向有效应力模拟结果Fig.14 The simulation results of vertical effective stress at the Shengli oilfield G89-1 block

为了大幅度提高CO2地质封存大规模模型数据处理和数值计算的效率,可基于MPI,OpenMP,CUDA 等工具研发高效的并行求解技术.其中,MPI主要适用于分布式存储的可拓展的并行计算机和工作站机群;OpenMP 主要适用于共享内存多处理器系统和多核处理器等体系结构,且与MPI 类似,适用于粗粒度任务并行计算;而CUDA 适用于包含GPU 设备的单个计算节点,或者与MPI/OpenMP 协同应用于包含多个GPU 设备的集群,主要适用于细粒度数据并行计算[104].

5 结论

储层-盖层系统的力学分析和稳定性评估是实现CO2地质封存工业化的关键问题.CO2注入地层后会导致较大的流体压力积聚、有效应力场变化,引发CO2泄漏、地表位移、盖层破坏和地震诱发等地质风险.盖层破坏的常见形式包括由剪切破坏引发的断层活化和由拉伸破坏引发的盖层压裂(裂缝扩展).盖层裂缝的起裂、扩展和断层的活化、滑移,形成了CO2逸出的潜在通道,对盖层完整性构成了重大威胁.CO2泄漏的主要机理是浓度驱动的分子扩散和压力驱动的对流,流动路径主要受裂缝网络几何形状控制,但也受到地层特性和地球化学反应等影响,涉及温度场-渗流场-力学场-化学场的多物理场耦合.

CO2地质封存数值模拟的研究对象是岩石固体骨架与孔隙流体(CO2、水和油等)组成的多相多组分系统,其数学模型主要包括基本控制方程和相平衡计算方法等.基本控制方程包括应力平衡方程、质量守恒方程、能量守恒方程和化学反应等.对离散后的控制方程耦合求解,求解方法大体可以分为全耦合、迭代耦合、弱耦合、显式耦合和拟耦合等算法.

为了适应CO2地质封存储层-盖层体系精细化建模、大时间尺度和大空间尺度的多物理场耦合对计算模型和计算效率的要求,还需要在裂缝性介质数值模拟方法方面进行重点研究.目前,三维地质体精细化建模通常都采用角点网格.因此,需要将角点网格和有限元方法的优点相结合,建立基于角点网格的地应力有限元模拟方法.将离散裂缝模型、多重介质模型和整体嵌入式离散裂缝模型相结合,建立相应的T-H-M-C 多场耦合复合介质数学模型,能够更准确地表征裂缝型储层并提高数值模拟的计算效率.由于CO2地质封存全地层模型的应力计算区域远大于渗流和温度计算区域,采用分区域的有限元-有限体积混合数值离散方法,能够在保证渗流场和温度场计算精度的同时,降低应力场求解的计算量.由于CO2地质封存问题涉及的时间尺度和空间尺度较大,在多物理场耦合方面需要发展高效的并行求解技术.