《1、 引言》

1、 引言

在《巴黎协定》框架下,到21世纪中叶实现碳中和是应对气候危机的根本解决方案。与电气化、氢能和可持续生物能源一样,碳捕集、利用与封存(CCUS)是实现净零能源系统的关键技术。其中,CCUS是唯一一组既有助于消除难以减少的CO2排放又有助于直接减少关键部门排放的技术[1]。在地质封存中,CO2在发电厂或工业设施中被捕集,经过压缩之后作为超临界流体输送并注入枯竭的油与气储层和含盐深水层,以便长期封存。CCUS技术极具潜力,是因为地质岩层具有的巨大全球存储能力、潜在注入地点的可用性和与发电厂的接近性。然而,这种注入的方法存在重大的技术挑战,在控制成本的情况下,该方法难以确保操作的安全性,也难以最大限度减少较长时间尺度上(几百年,甚至几千年)的地质断层和裂缝可能导致的泄露。在应对这些挑战时,准确预测注入CO2的行为对于碳封存的长期成功非常重要,因为即使是长时间的微小泄漏率也可能导致CO2封存的失败[25]。此外,储层地质力学响应与CO2高注入率对天然裂缝的影响有关,对确定地震效应至关重要[67]。在这里,需要一个用于准确预测注入CO2的命运和地质构造相关响应的物理模型,这一物理模型必须由各化学反应控制,包括多相流、多组分运输、岩石力学、热力学相行为、流体和岩石内的化学反应,并且还要由所有这些现象在多个时间和空间尺度上的耦合控制。该研究对所用的物理模型及其相应的数值近似值具有较高的精度要求。除地质碳封存外,CO2提高采收率(EOR)是通过向油藏注入CO2来提高采收率的另一种关键方法。最近的研究工作[8]表明,封存驱动的CO2 EOR中封存的CO2量可能超过燃烧石油所产生的CO2排放量,导致CO2净排放量为零,甚至为负。

在该研究中,我们使用得克萨斯大学奥斯汀分校地下建模中心开发的隐式并行精确储层模拟器(IPARS)来模拟地质碳封存过程。IPARS提供了一个框架,支持多个物理模型独立的和耦合的模拟。这些模型包括单相流和两相流、黑油模型、状态方程(EOS)组分流模型、热力学模型和地质力学模型等[913]。IPARS已被广泛用于与地质配方中CO2封存相关的实地基准和实验研究[1316]。具体而言,EOS组分流模型被用于模拟含盐深水层中CO2的封存过程。考虑到CO2会溶解到盐水中,CO2和盐水在组分模型中都建模为碳氢化合物组分。Peng-Robinson EOS [17]用于CO2-盐水相行为和性质的计算。随后进行闪蒸计算,以确定两个平衡阶段中CO2和水的摩尔分数。特别地,该模块被增强为可以模拟CO2地质封存应用中的CO2-盐水相行为[1819]。组分流模块包含先进的岩石物理模型,即包括新型三相相对渗透率滞后模型和隐式结构泡沫模型,用于准确描述循环注入过程和泡沫辅助回收过程[5,2021]。这些模型考虑了地质碳储量过程中的残余气体捕获效应。IPARS还被用于耦合流模拟和力学模拟,以研究地质碳储量过程中的孔隙弹性和孔隙塑性效应[5,22]。IPARS中的流动模块是在一般的六面体网格上使用质量保守的多点通量混合有限元方法(MFMFE)制定的[10]。此外,IPARS支持多种预调节和非线性求解器,如稳定双共轭梯度(BCGS)法、代数多重网格(AMG)法和广义最小残差(GMRES)法。该模拟器可以通过消息传递接口(MPI)库处理数百万个网格块,用于高性能计算环境中的并行分布式内存计算[11]。

储层井管理中的设计和控制问题具有挑战性,因为前向问题涉及耦合多物理场模拟,这些模拟具有高度非线性、多尺度、计算成本高且地质不确定性的特征。从优化的角度来看,无导数优化是很流行的。例如,这些方法包括粒子群、模拟退火(SA)、遗传算法(GA)和进化策略(ES)[2325]。这些全局优化方法具有使用现有储层模拟软件的灵活性;因此,它们已被广泛用于案例研究。另一组方法利用梯度的随机近似,如同步扰动随机逼近(SPSA)[26]和随机单纯形近似梯度(StoSAG)[27]。伴随方法也已应用于储层井优化[2829]。然而,大规模应用程序对计算的要求仍然很高。

最近,为了加快计算速度,机器学习、降阶建模和其他类型的代理模型也已被广泛集成到地下模拟中[3032]。当目标函数评估成本昂贵或可能难以处理时,贝叶斯优化(BO)成为控制和优化问题的强大解决方案。BO已被证明在超参数调优的机器学习[33]以及材料设计、机器人技术和环境监测[3436]等多个科学领域,特别是针对低维问题取得了成功[37]。与其他无导数方法相比,BO利用了机器学习技术,即高斯过程回归,允许在执行贝叶斯推理时分析可处理性。通过将采集函数与统计推理相结合,在搜索代理解决方案空间时可进行探索-开发权衡。BO已被用于基准储层模型的不确定性量化[38]。本研究采用了BO框架,利用美国密西西比州克兰菲尔德储层的数据,优化地质碳封存中的注入井调度。

本文的组织结构如下:第2节介绍了碳封存过程的数学建模;第3节介绍了BO及其与储层模拟器的耦合;第4节讨论了实验和结果;第5节给出了结论。

《2、 碳封存过程的数学建模》

2、 碳封存过程的数学建模

《2.1 组分流模型》

2.1 组分流模型

在该研究中,我们采用IPARS中的平行、EOS和组分流模块来模拟枯竭储层或含盐深水层中的地质碳储量过程。我们假设等温流、水和Nc烃组分形成三相流系统,即水相、非水液相和气相。ϕ表示当前的流体分数(即孔隙率),α表示α相的饱和度,ρα表示α相的质量密度,表示α相中组分i的注入/生产率,表示α相中组分i的摩尔分数,α表示α相的体积速度,表示α相中组分i的扩散-弥散张量。多相流系统的组分i的质量守恒方程表示如下:

tαϕSαραξiα+αραξiαuα-ϕSαDiαραξiα=αqiα,   i=1,,Nc; α=o,w,g(1)

式中,t是时间变量;o是油相;w是水相;g是气相。

κ表示多孔岩石基质的绝对渗透率张量,krα表示相对渗透率,μα表示黏度,pα表示参考相压,g是重力加速度,z是深度。用达西定律描述α相的体积速度uα

uα=-krαμακpα-ραgz,   (α=o,w,g)(2)

对毛细管压力方程、辅助方程以及初始和边界条件进行了补充,使得系统趋于关闭[3940]。用于求解该系统的主要变量是参考相压pαNc的组分浓度S1,S2,,SNc。水相和烃相分别使用微可压缩型和立方型Peng-Robinson EOS [17]。相位平衡是使用Rachford-Rice方程和等逸性准则计算的。

《2.2 相对渗透率滞后和泡沫模型》

2.2 相对渗透率滞后和泡沫模型

从油藏注气中吸取的经验之一是,由于低气体黏度、密度和地层非均质性引起的重力覆盖,气体体积波及效率本质上较低。因此,提高体积波及效率的方法很有价值。气水交替注入(WAG)可以降低气体的流动性。在注入水段塞的过程中,气体被毛细作用力困在高气体饱和区域,因此气流在随后的气段塞中被转移到气体饱和度较低的区域[41]。然而,WAG的有效性有限,并且WAG在高度非均质或厚的储层中效率不高。因此,可以引入表面活性剂-交替气体(SAG)作为降低气体流动性的第二种补救措施。原位产生的泡沫增加了气体的表观黏度,并阻止了高渗透率条纹,从而将流量从高渗透率转移到低渗透率和未波及区域[42]。

UTHYST模型[20]提供了一种简单的方法来计算任何相的相对渗透率的周期依赖性不可逆滞后行为。它模拟了高捕集数下单调增加的捕集饱和度,包括卡普利尔去饱和效应,并在三相多孔介质流中使用了动态Land系数。有关模型和验证的更多详细信息见参考文献[20]。

基于计算机建模小组(CMG)STARS中的隐式结构泡沫模型的改进泡沫模型最近已在IPARS中开发和实施[21]。该模型提出了一种新的泡沫生成函数,该函数能够捕获:①粗泡沫向强泡沫的转变,②泡沫剪切稀化流变性,以及③通过实验观察到的泡沫生成滞后现象。该模型已通过实验数据进行了校准和验证。更多详细信息见参考文献[21]。

《3、 贝叶斯优化》

3、 贝叶斯优化

考虑一个一般的优化问题:

maxxAf(x)(3)

式中,xRd表示d维控制变量;A是可行集合。BO利用高斯过程来表示目标函数f(x)。高斯过程是一个在函数上的分布,由其均值函数m0(x)和协方差函数k(x1,x2)定义:

fx𝒢𝒫m0x,kx1,x2(4)

它对可能的目标函数进行先验信念,并在训练期间通过更新基于观测数据条件的贝叶斯后验,迭代优化模型。协方差函数k(x1,x2)可以采取多种形式;此处我们使用Matérn内核[4344]:

kx1,x2=21-νΓν2ννHν2νx1-x2(5)

式中,Γ()和Hν()分别是ν阶伽马函数和贝塞尔函数。

用高斯过程表示对目标函数f的置信度,用采集函数(表示为a:𝒳R+)选择接下来通过代理优化评估的点:

xn+1=argmaxxax;𝒟n(6)

式中,𝒟n={(xi,yi)}i=1n表示使用yif(xi)进行第n次迭代中的观测值集。

预期改进(EI)是广泛使用的采集函数之一[4546]:

aEIx;𝒟nmnx-fxbest-εΦmnx-fxbest-εσnx+σnxϕmnx-fxbest-εσnx(7)

式中,xbestx1, x2,...,xn中的最佳候选;σnxn次迭代的后验标准差;ϕ()是标准正态分布的概率密度函数;Φ()是标准正态分布的累积分布函数。EI鼓励探索和开发,直观地说,当替代均值高(开发)和替代方差高(探索)时,EI的值更高。超参数ε>0,使探索-开发之间的权衡具有更大的灵活性[46]。

参考文献[47]中提出了情境改进(CI)的概念,它是一种简单而有效的启发式方法,用cvσn¯xfxbest替换了EI中的超参数ε

aCIx;𝒟n=mnx-fxbest-cvΦmnx-fxbest-cvσnx+σnxϕmnx-fxbest-cvσnx(8)

参考文献[47]表明,CI对较差的先验更为稳健。

综上所述,在算法1中给出了BO的伪代码[48]。

《算法1》

算法1 BO的伪代码[]

1:

Initialize the probabilistic surrogate with n0 points, D0=x1,y1,x2,y2,,xn0,yn0;

2:

for n=0 to N do

3:

 Select new query point xn+1 by optimizing acquisition function a, xn+1=argmaxxax;Dn;

4:

 Evaluate yn+1=fxn+1;

5:

 Augment data Dn+1=Dnxn+1,yn+1;

6:

Update the probabilistic mode.

《3.1 并行和分布式贝叶斯优化》

3.1 并行和分布式贝叶斯优化

在研究中,我们采用了基于并行汤普森采样(TS)的BO批处理实现[49]。TS的获取函数表示如下[50]:

aTSx;𝒟nfnx(9)

式中,fθfn𝒢𝒫(θ|𝒟n)参数化。

该算法描述见算法2 [49]。它可以以完全并行和分布式的方式实现;因此,它可以充分利用超级计算机上的多个处理器来处理[49]。

《算法2》

算法2 基于TS的并行和分布式BO []

1:

Initialize the probabilistic surrogate with n0 points, D0={(x1,y1),(x2,y2),,(xn0,yn0)}. Batch size B.

2:

for n=0 to N do

3:

 Compute current posterior p(θ|Dn).

4:

for b=1 to B do

5:

 Sample a random realization of the objective function from its posterior by sampling θ from p(θ|Dn)

6:

 Select new query point xn+1,bargmaxxDnfθ(x)

7:

 Evaluate yn+1,b=f(xn+1,b)

8:

 Augment data Dn+1=Dn(xn+1,yn+1)b=1B

《3.2 与储层模拟器耦合》

3.2 与储层模拟器耦合

国际商业机器(IBM)公司的贝叶斯优化加速器(BOA)托管在IBM的服务器上。它通过接口函数功能耦合到储层模拟器IPARS,IPARS在内部集群上运行。图1展示了BOA-IPARS耦合的原理图。对于每个批次,储层模拟在本地计算机集群上并行运行;模拟完成后对结果进行后处理,以计算目标函数值。然后将目标函数值传递给IBM云以进行BOA计算;接下来,将BOA建议的新评估点返回到集群,并为下一批次启动模拟。

《图1》

图1 BOA-IPARS耦合原理图。

《4、 数值实验》

4、 数值实验

《4.1 克兰菲尔德储层模型》

4.1 克兰菲尔德储层模型

克兰菲尔德站点位于美国密西西比州,是一个枯竭的储层,具有盐芯圆顶地质结构。2009年12月,得克萨斯大学经济地质局(BEG)研究团队在枯竭的克兰菲尔德站点启动了实地规模的CO2封存试点项目[15]。到2015年,克兰菲尔德的场地已经注入了约50万吨的CO2。储层模型包含80 ft × 7200 ft × 7200 ft(1 ft = 0.3048 m)的区域,深度为9901 ft,该区域被离散化为由20×144×144个单元组成的扭曲六面体网格,块尺寸为4 ft × 50 ft × 50 ft。该场地模拟了5个注入井。由于流动模型需要无流动边界条件,因此将8口拟生产井分配到储层边界,以模拟开放边界条件(图2)。拟生产井被限制在等于初始储层压力的恒定井底压力下,即4650 psi (1 psi = 6.894757 kPa);初始储层温度为125 ℃。假设该储层最初已被盐水饱和。考虑到CO2溶解到盐水中的影响,CO2和盐水在组分模型中都建模为碳氢化合物组分。Peng-Robinson EOS [17]用于CO2-盐水相行为和属性计算。接下来进行闪蒸计算,以确定两个平衡阶段中CO2和水的摩尔分数。

《图2》

图2 储层模型中的油井分布,5个注入井利用了现场的真实注水井建模,伪边界井模拟了储层模型三个边界上的开放流动边界条件。

表1总结了压力、体积和温度(PVT)数据,这些数据来自参考文献[15]。根据已发表的CO2在盐水中的溶解度和盐水[1819,51]密度的相关性,对EOS中使用的二元相互作用系数和体积位移参数进行了修正,该相关性是根据实验数据拟合得到的。使用了历史匹配的全张量渗透率、孔隙率、相对渗透率和毛细管压力曲线[15](图3图5)。克兰菲尔德储层数据已被用于许多地质碳储量的数值模拟研究[5,13,15,22]。

《表1》

表1 克兰菲尔德CO封存的EOS数据[]

ComponentCritical temperature (K)Critical pressure (psi)Compressibility factorAcentric factorMolecular weight (mol-1)Volume shift parameterBinary interaction coefficientDensity (kg·m-3)Viscosity (×10-3 Pa·s)CO2 mole fraction in brine
CO2304.131070.40.2550.22444.01-0.20000.09576.220.0440.013
Brine647.093540.90.2000.22418.010.29600.091033.290.440

《图3》

图3 储层初始孔隙度。

《图4》

图4 对数尺度下的储层垂直渗透率(logpermx)。md:毫达西。

《图5》

图5 相对渗透率和毛细压力曲线。krg:气相相对渗透率;krw:水相相对渗透率;Pc:毛细压力。

《4.2 基线模拟》

4.2 基线模拟

通过模拟表2中总结的4种注入方案:连续注入CO2、无相对渗透率滞后建模的WAG、具有相对渗透率滞后建模的WAG和具有相对渗透率滞后建模的SAG,来验证气体流动性控制技术(包括WAG和SAG [25])的效果。在4个模拟研究中,气体注入是由速率控制的,使得4种情况下注入的CO2总量相等。该设计对每种进样和模拟策略的CO2封存量进行无偏比较。图6显示了4种情况下进样结束时CO2分布的模拟结果(顶部和底部的视图)。正如预期的那样,在连续注入CO2期间[图6(a)],由于气体密度和黏度低,气体羽流向上迁移,直到到达封闭盖层;然后它会水平迁移,导致储层底部的波及效率差。从封存的角度来看,这种情况是不利的:经过16年的注入后,大部分注入的气体会通过边界井溢出。通过比较案例[图6(b)、(c)]展现出的性能,说明了相对渗透率滞后模型的显著效果。两种情况都模拟了20年的WAG流程;图6(c)中的案例支持渗透率滞后模型,而图6(b)中的案例则不能。与图6(b)中的案例相比,图6(c)中的案例的预测显示,其储层底部的CO2波及面积更广。相对渗透率滞后解释了气相的毛细管捕获现象,导致预测更多的CO2封存量。如表3所示,在不考虑相对渗透率效应的情况下,与连续注入CO2相比,WAG过程导致累积CO2封存量减少了9%;当考虑相对渗透率效应时,WAG过程导致累积CO2封存量增加55%。

《表2》

表2 CO注入方案总结

DescriptionInjection scheme
Gas (year)Water (year)Gas (year)Water (year)Gas (year)Water (year)Gas (year)Water (year)Gas (year)
CO2 injection16
WAG without hysteresis413131313
WAG with hysteresis413131313
SAG with hysteresis413131313

《图6》

图6 每种注入情况下,进样结束时CO2气体饱和度分布(Sgas)的模拟结果。(a)连续注入CO2;(b)无相对渗透率滞后的WAG;(c)具有相对渗透率滞后的WAG;(d)具有相对渗透率滞后的SAG。

《表3》

表3 不同注入方案的CO封存量总结

DescriptionCO2 storage volume (Mscf)Percentage of change (%)
CO2 injection0.64 × 108
WAG without hysteresis0.58 × 108-9
WAG with hysteresis0.99 × 108+55
SAG with hysteresis1.13 × 108+77

SAG在气体流动性控制中的突出效果如图6(d)所示。由于原位产生的泡沫降低了液/气迁移率,气相可以侵入低渗透区域并克服重力分层。与连续注入CO2和WAG相比,垂直气体波及效率显著提高。表3汇总了各注入案例中的CO2封存量,结果表明,与连续注入CO2相比,SAG的CO2封存量显著提高了77%。在下一小节中,表2中注入方案的SAG模拟被设置为注入周期优化的基线。

《4.3 用于气体/表面活性剂注入计划的贝叶斯优化》

4.3 用于气体/表面活性剂注入计划的贝叶斯优化

目标是优化气体/表面活性剂注入周期,以便在特定注入周期内最大化上述储层模型中的累积碳储量。T表示总注入时间,rinj,i表示注入井i的CO2注入速率,Ninj表示注入井的总数,rprod,i表示生产井i的CO2生产率,Nprod表示生产井的总数。则最优化问题表述如下:

argmaxC1,C2,, CNS1,S2,, SN k=0Ki=1NinjSkCk+1rinj,idt-i=1Nprod0Trprod,idt  

s.t.  0(10)

式中,控制变量C1,C2,,CN+1是每个CO2注入周期的进样停止时间(其中,CN+1=TK是循环数);S0,S1,S2,,SN是每个表面活性剂注入周期的进样停止时间(其中,S0=0)。最多4个(Ninj=4)的注入周期和20年的总注入时间(T=20)是固定的。

在BO实验中,我们使用配备了Matérn 52协方差核的高斯过程和CI获取函数。我们分配了15个函数评估来初始化高斯过程。除了初始化之外,计算预算被设置为最多200次函数评估。对于批量BO,我们使用的批量大小为5。我们还在同一优化问题上实验了GAs和协方差矩阵自适应(CMA)-进化策略(ES)[25];为了实现这些算法,我们使用了Python(DEAP)[52]中的开源ES软件包的分布式进化算法。

BO的最佳实验结果如表4所示。值得注意的是,基线SAG方案的目标是5个注入周期,第一个注入周期从气体开始。优化后的SAG方案只有两个注入周期,第一个注入周期从表面活性剂开始。研究结果与原位泡沫生成的物理机制一致——表面活性剂注入的初始周期将有助于在整个注入方案的早期生成泡沫。与基线SAG方案相比,优化后的方案在20年的注入中使CO2封存量增加了16%,减少了56%的水和表面活性剂使用量。

《表4》

表4 通过BO优化的气体/表面活性剂方案与连续气体注入和基线注入方案的比较

DescriptionInjection schemeCO2 storage volume (Mscf)
Gas (year)Surfactant (year)Gas (year)Surfactant (year)Gas (year)Surfactant (year)Gas (year)Surfactant (year)Gas (year)
CO2 injection160.64 × 108
Baseline SAG41.03131.00313.001.13 × 108
Optimized SAG1.520.2516.251.31 × 108

为测试BO算法的稳健性,我们在重复实验中比较了BO算法和其他两种算法的性能。每种算法的实验都重复了三次。图7比较了三种算法在重复实验中的最佳运行性能。实验表明,BO得到了比GA和CMA-ES更好的解,且其函数评估次数更少。更具体地说,BO在73次功能评估后达到最佳观测值1.31×108 Mscf (1 Mscf = 28.317 m3),而CMA-ES在183次功能评估后才达到1.28×108 Mscf。BO以减少60%的函数评估实现了具有竞争力的目标函数值。这些结果表明BO算法的高效性。图8通过重复实验绘制了三种算法的整体性能。实线表示函数值的平均值,阴影区域表示标准差。这些重复实验证明了批量BO算法具有稳健性优势。

《图7》

图7 BO、GA和CMA-ES的最佳运行性能比较。

《图8》

图8 BO、GA和CMA-ES的整体性能比较。

《5、 总结和讨论》

5、 总结和讨论

本文提出的框架首次尝试将高保真物理模型与机器学习技术相结合,这些技术可以扩展到地质碳封存的现场应用。我们将高保真组分储层模拟器IPARS与IBM BOA相结合,优化注入井调度,以便使用原位生成的泡沫控制气体流动性,从而最大限度地提高累积CO2封存量。在此框架下,使用IPARS的正向模拟严格遵循了地质碳封存过程中的复杂物理特性。组分流模型包括相对渗透率滞后模型和泡沫模型,可以模拟CO2封存的三种主要捕获机理,即结构捕获、残余捕获和溶解度捕获。由于目标函数的计算代价很高,因而BO在设计和控制方面被证明是成功的。在这里BO由两个组件组成:一个是为目标函数构建贝叶斯代理项的高斯过程机器学习模型;另一个是决定采样点序列的采集函数。高斯过程使用贝叶斯统计进行分析,并对无法获得评估的未知目标函数值中的不确定性给予可处理性,采集函数用它来平衡探索-开发权衡。这两个组件的组合可以有效地搜索代理解决方案空间。

我们将该框架应用于实地碳储量,利用来自克兰菲尔德储层的现场数据,设计泡沫辅助碳封存的注入周期。数值实验结果表明,与基线情况相比,优化后进样方案的储气量增加了16%,水/表面活性剂的使用量减少了56%。我们进一步进行了比较研究,以检查BO的性能和稳健性,并将其与其他常用算法(包括GA和CMA-ES)进行比较。结果表明,BO算法具有较高的效率和稳健性,因为它以更少的函数评估次数实现了具有竞争力的目标函数值。结果表明,BO有望进一步应用于CCUS项目中的其他控制和井管理应用,如经济约束下的优化、井位布置优化等。由于IPARS包括与组分流模型有效耦合的地质力学模块,该框架可以轻松地扩展到包括碳存储过程中的地质力学效应。此外,该框架还可以扩展到地质不确定性情况下BO算法的应用。