地质碳封存中井底压力与CO2羽流预测的边界条件保真度

# 1 引言  
来源：https://arxiv.org/html/2606.27515  
地质碳封存中井底压力与CO₂羽流预测的边界条件保真度  

Romal Ramadhan¹, Seyyed A. Hosseini¹,², Larry W. Lake³  

¹ 德克萨斯大学奥斯汀分校杰克逊地球科学学院地球与行星科学系，美国德克萨斯州奥斯汀市，邮编78712  
² 德克萨斯大学奥斯汀分校经济地质局，地址：10611 Exploration Way，美国德克萨斯州奥斯汀市，邮编78758  
³ 德克萨斯大学奥斯汀分校Hildebrand石油与地质系统工程系，美国德克萨斯州奥斯汀市，邮编78712  
*通讯作者：[email protected], [email protected]  

###### 摘要  

准确预测井底压力（BHP）和CO₂羽流运移对于安全的深部地质碳封存至关重要，然而实际模拟常依赖截断域，其中人工边界会扭曲压力扩散和CO₂饱和度足迹。本研究通过将十种降阶边界处理方法与均质及非均质储层中的全域参考模拟进行对比，评估了边界条件保真度对BHP和CO₂羽流预测的影响。我们测试了均匀孔隙体积乘数、传导率修改器、角点校正孔隙体积修正、分层修正及渐变修改器，并采用BHP均方根误差（RMSE）、归一化均方根误差（NRMSE）、峰值压力偏差和羽流交并比（IoU）作为性能指标。结果表明，保持角点孔隙体积是截断域建模的最重要要求。我们发现忽略角点储量的均匀处理方法会产生较大的压力误差：在均质模型中BHP RMSE为362至382 psi，在非均质模型中为250至304 psi，并且羽流IoU值接近0.80至0.84，表明约16%至20%的组合羽流面积被错误表征。角点校正方案显著降低了压力误差，并将羽流IoU提升至0.94以上，但我们观察到传导率修正并非普遍有益。在均质储层中，均匀传导率调整改善了压力保真度；但在非均质储层中，它可能过度限制跨变渗透率边界面的流动，从而增大BHP误差并缩小预测羽流范围。我们发现，带有传导率修正的渐变修改器能提供最一致的性能，在两种储层类型中均实现了BHP NRMSE低于3.7%和羽流IoU超过0.97的结果。  

关键词：地质碳封存；井底压力；CO₂羽流；边界条件；孔隙体积乘数；传导率修改器；  

准确预测井底压力（BHP）和二氧化碳（CO₂）羽流运移对于地质碳封存（GCS）的安全与经济部署至关重要[8 (https://arxiv.org/html/2606.27515#bib.bib4)]。这两个响应变量涵盖了GCS风险中相互耦合的压力和空间维度：BHP控制着注入速率设计以及将注入压力维持在低于地层破裂梯度的水平，而羽流范围则定义了位移、溶解和捕获过程发生的地下足迹。两者均具有直接的监管权重[4 (https://arxiv.org/html/2606.27515#bib.bib5),10 (https://arxiv.org/html/2606.27515#bib.bib6)]。BHP是表明作业保持在安全压力限制内并保护地下饮用水源（USDW）的主要指标，而羽流范围则界定了美国环保署（EPA）VI类许可下的审查区域（AoR），在该区域内作业者必须识别并修复潜在的泄漏路径[9 (https://arxiv.org/html/2606.27515#bib.bib7),7 (https://arxiv.org/html/2606.27515#bib.bib9)]。在作业过程中，模拟与实测的BHP或羽流行为之间的偏差可能触发注入速率调整、额外的监测义务或监管审查，因此这两个维度的模拟保真度直接关系到项目的经济性和时间进度[3 (https://arxiv.org/html/2606.27515#bib.bib8)]。  

实现这一保真度因GCS模拟中固有的尺度不匹配而变得复杂[5 (https://arxiv.org/html/2606.27515#bib.bib10)]。CO₂注入产生的压力扰动在数十年的时间尺度上从井筒向外传播数十至数百公里，远远超出高分辨率数值模型实际可捕获的空间范围，而羽流本身根据渗透率、非均质性和浮力作用可能仅运移数百米到几公里[1 (https://arxiv.org/html/2606.27515#bib.bib11)]。因此，从早期筛选到详细表征和作业监测，在GCS工作流程的几乎每个阶段，域截断都是实践上的必要手段[6 (https://arxiv.org/html/2606.27515#bib.bib3)]。然而，截断的代价是引入人工边界，这些边界将压力反射回注入井，而非允许其向远场消散。这种反射使模拟的BHP高于无限域应有的值，并通过人为约束驱动侧向运移的压力梯度来扭曲羽流轨迹。  

已有几种基于修改器的边界处理方法被提出以减轻这些影响。[2 (https://arxiv.org/html/2606.27515#bib.bib1)] 表明，在边界单元上应用孔隙体积乘数（PVM）和传导率修改器（TM）方案会显著影响BHP和审查区域预测，并且渐变修改器方案能改善与全域参考行为的一致性。[11 (https://arxiv.org/html/2606.27515#bib.bib2)] 也证明了将边界单元扩展与等效的PVM和TM相结合，比简单的定压或过大孔隙体积简化能提供更多增益，特别是对于区域性的开阔含水层，对羽流运移也有较小但仍然相关的影响。在我们之前的工作中，我们为储量守恒的截断模型引入了一个渐变修改器框架，并通过在边界单元上分布孔隙体积和传导率修正，展示了更优的AoR预测效率[6 (https://arxiv.org/html/2606.27515#bib.bib3)]。尽管取得了这些进展，现有研究通常孤立地评估压力或羽流行为，且尚未有研究确定是否存在一种单一的边界处理方法能够同时在从均质到非均质再到全尺度的储层类型中保持两者。  

这里还存在另一个我们之前未能明确的误差来源。角点单元位于两个域边界相遇之处，它们最终代表了均匀修改方案简单忽略的、出乎意料大的储层体积部分。当我们失去这些角点体积时，我们也失去了注入井所依赖的部分压力缓冲，从而推高BHP并在边界附近引入梯度误差。这些误差并非固定不变；它们会渗透到模型预测的羽流运移路径中。一个提出的修正方案是对修改器进行空间分级，将修正作为距域中心距离的幂律函数进行分布。其思路是恢复缺失的角点体积，同时平滑均匀方案所产生的尖锐属性对比。我们尚未证明的是，这种几何修正是否能在我们预期实践中遇到的各种储层结构范围内，为BHP和羽流行为提供更好的预测。  

我们通过在全域参考模拟的基础上，对均质和非均质储层模型中的边界条件配置进行评估来填补这些空白。使用BHP时间序列和CO₂羽流范围作为同等重要的验证标准，我们在整个注入和注入后监测期内量化了峰值压力偏差、归一化均方根误差（NRMSE）和羽流面积指标的性能。结果为GCS项目各阶段的边界条件选择提供了实用指导，对注入设计、审查区域划定和监管合规证明具有直接影响。  

## 2 方法  

### 2.1 数值模型描述  

我们构建了地质复杂性递增的储层模型，以评估边界条件如何影响BHP和CO₂羽流预测。均质模型作为不受地质变异性影响的受控基准，而非均质模型则引入了空间复杂性。两个模型共享相同的域几何形状和注入计划，以便将BHP与羽流行为的差异完全归因于边界处理和储层非均质性。我们使用CMG-GEM组分模拟器（Computer Modelling Group Ltd.）在等温条件下进行所有模拟。排除了地球化学反应和地质力学耦合，以使主导物理机制仅限于压力扩散和浮力驱动的CO₂运移。表1 (https://arxiv.org/html/2606.27515#S2.T1) 汇总了关键模型参数。  

#### 2.1.1 均质模型  

我们将全域参考案例（以下称为“真实模型”）离散化为一个501×501×1的网格，单元尺寸为113.65 ft × 113.65 ft × 328.1 ft，覆盖侧向面积为56,939 ft × 56,939 ft。地层顶面深度为3,281 ft，均匀厚度为328.1 ft。我们将储层初始压力设为1,600 psi，并在整个网格中保持孔隙度和渗透率恒定，分别为0.20和10 mD。一个中心注入井以地表当量速率2,500吨/天注入CO₂，持续10年，之后关井并继续监测系统40年。  

我们围绕注入井同心提取了两个截断变体：251×251×1和101×101×1个单元。所有内部单元保持与真实模型相同的属性；我们仅根据边界处理方法修改最外层单元环。  

#### 2.1.2 非均质模型  

为了在全域501×501×1网格上生成空间相关的一维孔隙度和渗透率场，我们实现了序贯高斯模拟。然后对实现进行重新缩放，使其算术平均值与均质值（孔隙度0.20，渗透率10 mD）一致。图1 (https://arxiv.org/html/2606.27515#S2.F1) 展示了真实模型以及两个降阶域提取结果中的孔隙度和渗透率分布。降阶子模型（251×251×1和101×101×1个单元）通过使用NULL关键字从真实网格中切割得到，该关键字会去激活感兴趣区域之外的单元，同时保留活动域内全分辨率属性场。我们应用与均质情况相同的边界层修改；所有其他模拟设置与真实模型保持一致。  

表1：全域模拟的储层模型属性（根据[6 (https://arxiv.org/html/2606.27515#bib.bib3)]）  
[表内容]  

图1：非均质储层的孔隙度（上行）和渗透率（下行，对数尺度）空间分布。左侧为全域真实模型（501×501），中间和右侧面板分别为251×251和101×101的降阶域提取结果。黑色轮廓和连线表示嵌套的子域边界（修改自[6 (https://arxiv.org/html/2606.27515#bib.bib3)]）。  

### 2.2 边界条件公式化  

我们对每个边界单元的孔隙体积进行缩放，乘数考虑了超出网格边缘被移除的储层体积。对于251×251和101×101网格，我们将均匀PVM值分别设为126和201（后文称为cc），以使边界单元的总孔隙体积等于全域被切除部分的体积[6 (https://arxiv.org/html/2606.27515#bib.bib3)]。然后，我们根据相邻两个单元各自半单元传导率的调和平均值计算它们之间的传导率。对于外边界上j方向的面对面传导率[2 (https://arxiv.org/html/2606.27515#bib.bib1)]：  

T_{j+1/2} = \frac{K_a (\Delta x_a \cdot \Delta z_a)}{\Delta y_a / 2}   (1)  
T_{j-1/2} = \frac{K_b (\Delta x_b \cdot \Delta z_b)}{\Delta y_b / 2}   (2)  

其中Ka和Kb是单元渗透率，Δx、Δy、Δz是单元在x、y、z方向的尺寸。界面处的调和平均传导率为：  

T_j = \frac{1}{\frac{1}{T_{j+1/2}} + \frac{1}{T_{j-1/2}}}   (3)  

当对边界单元施加PVM时，必须调整传导率以反映修改后的单元几何形状。在i和k方向，由于共享界面的两个单元都承受相同的体积缩放，调整后的半单元传导率变为：  

T_{i+1/2}^{\mathrm{adj}} = \frac{K_b (\mathrm{PVM} \cdot \Delta y \cdot \Delta z)}{\Delta x / 2}   (4)  
T_{i-1/2}^{\mathrm{adj}} = \frac{K_c (\mathrm{PVM} \cdot \Delta y \cdot \Delta z)}{\Delta x / 2}   (5)  

由于PVM对称地作用于两项，i和k方向的传导率修改器简化为：  

TM_{i,k} = \frac{T_i^{\mathrm{adj}}}{T_i} = \mathrm{PVM}   (6)  

在j方向，界面连接一个修改后的边界单元（应用了PVM）和一个未修改的内部单元（PVM=1），形成不对称配置。边界侧的调整后半单元传导率为：  

T_{j-1/2}^{\mathrm{adj}} = \frac{K_b (\Delta x_b \cdot \Delta z_b)}{\mathrm{PVM} \cdot \Delta y_b / 2} = \frac{T_{j-1/2}}{\mathrm{PVM}}   (7)  

在 dissimilar 单元之间进行调和平均，得到传导率修改器：  

TM_j = \frac{T_j^{\mathrm{adj}}}{T_j} = \frac{2}{\mathrm{PVM} + 1}   (8)  

位于域角点的单元同时受到两个截断边缘的影响。因此，它们占缺失储层体积的比例比沿单一边界的单元大得多，均匀的PVM校正无法完全恢复这一损失。图2 (https://arxiv.org/html/2606.27515#S2.F2)a显示，仅对边缘进行缩放仍会留下四个角点区域补偿不足，而图2b则标识了需要调整的传导率面。为了修正这一剩余缺陷，角点校正方案（VA和VAT）对每个角点单元的PVM额外施加因子c²，其中c是基础乘数。此调整确保降阶模型保留真实模型的总孔隙体积。  

图2：传统边界单元修改中角点孔隙体积亏缺示意图：(a) 应用于最外层单元环（红色）的均匀PVM方案未考虑四个角点区域（标记为“缺失PV”）；(b) 相应的传导率调整