非均质润湿性会改变浅层含水层中甲烷的迁移和泄漏情况

《Environmental Science & Technology》：Heterogeneous Wettability Alters Methane Migration and Leakage in Shallow Aquifers

　　高分辨率图像
下载MS PowerPoint幻灯片

毛细管异质性越来越被认为是控制气体羽流在含水层和储层中迁移和捕获的主要因素。我们发现，由矿物学和盐度决定的水-甲烷接触角的空间变异性会重塑毛细管进入压力，进而影响迁移路径。通过分子动力学模拟，我们估算了石英和高岭石在淡水和咸水条件下的接触角，并将这些结果通过基于接触角的Leverett J函数嵌入到连续介质多相流模拟中，将润湿性直接映射到连续介质尺度上的流动特性上。考虑到接触角的异质性，甲烷的行为发生了变化：含水层中可移动和残留的甲烷量减少了多达10%，而向大气的泄漏量增加了多达20%。这种效应的幅度与渗透率对比、泄漏率、盐度和岩相比例有关。通过将分子尺度上的润湿性与连续介质尺度上的流动和传输过程结合起来，这种跨尺度框架为地下水保护和风险评估提供了更坚实的物理基础，并产生了更可靠的排放量估计。该方法可以推广到其他地下气体传输问题，包括氢气和二氧化碳的储存以及永久冻土融化产生的甲烷自然释放。

1. 引言
甲烷（CH4）从地下源泄漏到浅层含水层和大气中是一个重要的环境问题。(1?4) 逃逸的甲烷可能来源于受损的井、天然断层、裂缝或地表活动，并且还受到天然生成的生物甲烷的复杂影响。(5?11) 作为一种强效的温室气体，甲烷对地下水质量、气候和公共安全构成威胁，因此预测其迁移路径至关重要。(12?15) 最近的研究强调了小尺度地质特征的影响，其中毛细管异质性被认为是控制羽流迁移和捕获的主要因素。(16,17) 忽视毛细管压力的变异性会导致预测偏差，通常会高估大气排放量并低估含水层的滞留能力。(16,18) 在此基础上，我们研究了依赖于矿物和盐度的润湿性（表现为水-甲烷接触角的异质性）如何在分子尺度到含水层尺度之间进一步重塑泄漏路径。

实地研究和数值模型已经证明，小尺度结构特征和毛细管性质的变化可以改变羽流路径、改变捕获效率并影响泄漏结果。(16,17,19) 毛细管异质性已被认为是控制气体迁移的关键因素。(16,20?26) 然而，实地研究表明，忽略这一因素可能会在气体羽流迁移预测中引入系统偏差。(20?22,26,27) 除了基于孔隙度和渗透率的毛细管压力缩放方法外，润湿性代表了另一个未被充分探索的异质性维度，对羽流行为同样具有重要影响。

传统的甲烷泄漏研究及相关问题（如CO2和H2储存）中的毛细管压力缩放通常依赖于孔隙度和渗透率的空间变化，同时假设界面性质保持不变。(16,28?33) 这种方法忽略了润湿性取决于矿物学和水化学（如盐度）的事实，而这些因素会改变局部毛细管进入压力，从而影响迁移路径，但在泄漏模型中很少被考虑。

直接测量接触角（例如静止/悬挂滴）很有价值，但面临诸如表面污染、不必要的粗糙度以及岩心制备过程中的平滑处理等已知挑战。因此，测得的接触角可能无法真实反映自然表面条件。这些挑战在未固结的浅层含水层中尤为突出。(34,35) 分子动力学（MD）模拟通过解析矿物-流体-流体界面处的原子尺度相互作用提供了一种补充方法。先前的研究表明，MD计算出的接触角与实验室测量的接触角在多种岩石-流体系统中具有良好的一致性。(36?40) 通过这种方式，MD模拟可以提供特定岩相的润湿性输入，这些输入可以通过Leverett J函数传递到连续介质尺度模型中，从而在纳米尺度表面物理和连续介质多相流之间建立明确的联系。

在这里，我们通过将基于MD模拟计算的接触角（作为矿物类型和盐度的函数）整合到浅层含水层井筒中甲烷泄漏的连续介质多相流模拟中，探索了微观尺度润湿性与连续介质尺度泄漏行为之间的联系。这种整合使用了基于接触角的Leverett J函数，并评估了明确表示异质润湿性对甲烷迁移、捕获和泄漏的影响。我们进一步通过改变盐度、渗透率对比、泄漏率和地质上合理的岩相比例，确定了这种表示方法最为重要的条件。

2. 材料与方法
将MD模拟计算的、依赖于矿物和盐度的接触角与连续介质多相流模拟相结合分为两个阶段进行。首先，MD模拟在不同盐度条件下提供了特定岩相的接触角。然后，将这些计算出的接触角纳入连续介质模型中以完善毛细管压力缩放。所得模拟用于量化甲烷向大气的泄漏，并确定量化接触角变异性对准确预测最为关键的条件。MD模拟和连续介质模拟的计算成本详见支持信息中的S1部分。

所有MD模拟均使用LAMMPS软件包（41）进行，采用MPI并行化和GPU加速。(42?47) 力场、平衡协议和分析程序的完整细节以及与先前实验室结果的比较详见支持信息中的S2部分。尽管浅层含水层包含多种矿物相，但先前的研究表明，小尺度矿物学异质性对溶质传输的影响有限，羽流的迁移和扩散主要由更大尺度的水文地质结构控制。(25) 基于这些证据，沉积物变异性被简化为两种主要单矿物相：粗粒（CG）相和细粒（FG）相（即沉积物体）。(48) 石英被选为CG相的代表矿物，而高岭石涂层表面代表FG相。

石英存在多种多形体，每种多形体都具有独特的原子排列，这取决于多孔介质的压力和温度条件。(49) 在这些多形体中，选择α-石英（具有(0001)表面），因为它在目标模拟条件下的稳定性较好。(49,50) 由于α-石英呈六方结构，而模拟盒子需要矩形几何形状，因此需要对晶体晶格进行修改。为此，将α-石英单元格(51)沿x和y方向复制，然后沿特定晶体学平面修剪，以确保二维周期性和基底的中性电荷。然后构建一个正交单元格，并沿x、y和z轴复制，保持整个结构的原子化学计量比和键合连接性。(52) 从生成的α-石英表面提取Q2切割面并完全羟基化，以模拟实际的水环境条件。对于Q2，硅原子与两个桥接氧原子和两个羟基(OH)团相连。这是一种不太稳定、高反应性的硅醇基团，其结构为(≡Si–O?)2Si(OH)2。相比之下，高岭石(Al2Si2O5(OH)4是一种不膨胀的1:1型粘土矿物，由一层四面体二氧化硅(T)和一层八面体氧化铝(O)组成。(53) 与蒙脱石不同，高岭石的层电荷可以忽略不计，其层间通过氧化铝层上的羟基与二氧化硅层上的基态氧原子之间的氢键稳定，形成刚性的TO–TO堆叠，层间膨胀最小。高岭石单元格首先沿x和y方向扩展，然后沿z轴复制三次以生成模拟域。

使用已建立的MD模拟协议估算这些矿物上的水-甲烷接触角，这些协议已被证明能够可靠地再现实验室规模的测量结果。(38?40) 在本研究中，将半圆柱形水滴放置在含有CH4的矿物基底上，并使其达到稳态几何形状（支持信息，图S1）。在模拟盒子中，x和y方向施加周期性边界条件以模拟无限系统，确保粒子从一个边界退出后可以从相反侧无缝重新进入。而在z方向施加非周期性边界条件，原子可能被反射、吸收或在模拟盒子内受限。为了最小化边界效应，在模拟盒子顶部放置了反射壁，并选择足够大的盒子高度以防止墙壁对接触角测量产生影响。模拟盒子的大小约为120纳米（包括基底）。接触角是在水滴、矿物表面和CH4相之间的三相边界处测量的，用于连续介质模拟。

储层基本体积（REV）的概念对于连续介质尺度储层模拟至关重要，因为它定义了润湿性和毛细管压力等性质变得具有代表性并适合进行尺度放大的最小尺度。(54) 在本研究中，通过使用MD模拟量化与实验室规模测量一致的接触角，明确解决了REV假设。MD衍生的接触角与已发表的实验数据之间的一致性表明，模拟可靠地捕捉了从孔隙尺度到实验室尺度的润湿性行为。由于连续模型中使用的毛细管压力函数本质上依赖于REV平均的润湿性参数，MD衍生的接触角处于适当的代表性尺度范围内，可以直接作为物理意义上的REV一致输入纳入连续介质多相流模拟中。

连续介质模拟使用Computer Modeling Group（CMG，版本2024.30）的GEM模块进行。连续介质模拟中使用的性质范围基于先前的现场研究和CH4在无限制浅层含水层中的数值模拟。(17,24) 模型域被划分为三个部分，代表一个无限制的浅层含水层及其上方的非饱和带和大气层（支持信息，S3部分）。初始饱和度指定为含水层中100%的水分、非饱和带中40%的水分和60%的空气以及大气中100%的空气。甲烷泄漏作为点源引入含水层，以模拟向上和横向迁移。空气被建模为假纯组分，假设在浅层地下条件下不溶解。(55) 使用CMG的Winprop模块校准流体性质和相行为（支持信息，S3.1部分）。为表示多相流动行为和潜在的向大气泄漏，分配了特定部分的相对渗透率和毛细管压力函数。

尽管毛细管压力曲线通常在岩相尺度上分配，(56) 但在单个岩相内部仍存在矿物组成、颗粒涂层、胶结作用和粘土含量的小尺度变化。与这种亚岩相矿物异质性相关的毛细管异质性可以通过我们的方法明确纳入。

每个网格块的毛细管压力缩放通过基于接触角的(Leverett J函数)纳入润湿性变异性：
$$
P_c^{\text{grid}} = \cos(\theta_{\text{grid}} \cos(\theta_{\text{ref}}) k_{\text{ref}} \phi_{\text{grid}} \phi_{\text{ref}}
$$
其中下标ref和grid分别表示参考条件和模型网格块值。k、?、θ和PC分别代表绝对渗透率、孔隙度、接触角和毛细管压力。

第一项（$\cos(\theta_{\text{grid}} \cos(\theta_{\text{ref}} \cos(\theta_{\text{ref}})$）考虑了岩相之间的润湿性变化，而第二项（$k_{\text{ref}} \phi_{\text{grid}} \phi_{\text{ref}}$）考虑了孔径相关长度的缩放。在没有接触角（NCA）的情况下，润湿性假设为常数，因此第一项等于1，缩放简化为大多数气体传输研究中使用的传统孔隙度-渗透率公式。

每次模拟都从30天的预泄漏阶段开始，在此期间，模型边界处的注入井和生产井以约0.6 m3/天的速率运行，以建立稳定的地下水流动。然后以0.14 rm3/天的恒定速率引入甲烷泄漏，持续60天，并将模拟延长至100天，以确保泄漏路径的稳定。

使用T-PROGS中实现的基于转移概率的方法对含水层和非饱和带内的CG和FG岩相的空间分布进行了建模(58)（见支持信息，S3.2部分）。该方法生成了25个随机实现的岩相结构，每个实现都代表了CG和FG岩相的不同空间配置，同时保持相同的体积比例。这些实现之间的差异在于岩相的几何形状、空间排列和连接性，从而捕捉地质不确定性并允许评估模型结果的随机变异性。

为了表示不同的水文地球化学条件，对零盐度（0 ppm NaCl）和咸水（12,000 ppm NaCl）地下水进行了MD和连续介质模拟。淡水情况代表主要通过降水补给的内陆含水层，而咸水情况代表沿海含水层、蒸发盆地或受农业补给影响的区域。(59?62) 对于基准情景，CG岩相的渗透率、甲烷泄漏率和CG岩相的体积比例分别设置为5,000 mD、0.14 m3/天和40%。从这个基线出发，通过系统地改变关键参数，开发了额外的情景：(i) 将CG渗透率从5,000增加到20,000 mD，同时将对数尺度方差从0.3增加到1.0；(ii) 将CH4泄漏率和总泄漏体积增加2.5倍；(iii) 将盐度从0增加到12,000 ppm，以评估地下水化学成分的影响；(iv) 将CG相的体积比例从40%增加到60%，以研究沉积变异性和连通性的影响。选择60%的CG相作为上限，因为在二维系统中，连通性阈值接近50%。然后比较了包含接触角（CA）和不包含接触角（NCA）情景下CH4在移动相、被困相和溶解相中的分布，以评估接触角对泄漏预测的影响。

3. 结果与讨论
3.1. 接触角估计
基于分子动力学（MD）模拟得到的接触角估计值列在图1中。为了确认模拟时间足够，跟踪了每种矿物在两种盐度水平（0和12,000 ppm NaCl）下水滴质心的标准化垂直位置。如支持信息中的图S2所示，所有情况都达到了稳定状态（石英和高岭石分别为6–10纳秒和2纳秒）。MD计算出的接触角落在之前实验室测量结果的范围内（支持信息中的图S3），这支持了该方法的可靠性。实验值之间的差异可能反映了样品制备、表面条件、测量技术和操作者处理方式的差异。我们的MD结果与这一实验范围的一致性进一步支持了建模方法的有效性；然而，也应承认实验室测量中的不确定性，并将其纳入分析中。

图1
图1. 使用分子动力学模拟基于生产阶段最后2纳秒内水滴形状的时间平均值估计的接触角：(a) 零盐度地下水–CH4–石英，(b) 盐水地下水–CH4–石英，(c) 零盐度地下水–CH4–高岭石，以及(d) 盐水地下水–CH4–高岭石。浅绿色线条表示用于确定接触角的水滴边界。

高分辨率图像
下载MS PowerPoint幻灯片

石英的亲水性始终高于高岭石。盐度的增加降低了两种矿物的亲水性，使石英的接触角（θ）增加了约14°（从约20°增加到约34°），高岭石的接触角增加了约3°（从约63°增加到约66°）。石英对盐度的更强响应在物理上是合理的（63–68），因为溶解的盐分在矿物-水界面屏蔽了静电电荷，减弱了水对表面的粘附力并增强了气体的润湿性。在零盐度和盐水条件下，使用MD模拟最后1纳秒内的10帧数据估计了两种矿物的接触角不确定性。对于石英，在零盐度地下水中接触角的标准差约为1.7°，在盐水中约为2.6°；而对于高岭石，分别为约3.7°和4.0°。

3.2. 连续介质尺度模拟
完成所有模拟后，我们首先检查了随着实现次数的增加，平均CH4（%）分布如何收敛。这项收敛性分析是在基础情景下进行190天模拟后进行的，考虑了不同储层部分和流体相的变化。目的是验证模拟结果的统计稳定性和代表性，确保观察到的CH4迁移和捕获行为不受矿物分布随机变异性的影响（支持信息，第S4.1节）。一旦确认了收敛性，就比较了CA和NCA情景下的羽流传播。

如图2所示，将接触角变异性纳入毛细压力缩放（CA）从根本上重塑了储层中毛细进入压力的空间分布。通过减少高压区和低压区之间的强烈对比，这种方法产生的毛细异质性模式比不包含接触角变异性（NCA）的传统缩放更加平滑和空间分布更均匀。在NCA情况下，粘土覆盖的区域保持高进入压力，成为毛细屏障，迫使CH4通过有限的优先路径迁移，导致羽流前沿尖锐。当包含接触角变异性时，这些区域的部分地方cos(θ)的降低降低了进入压力，使CH4能够侵入以前无法进入的孔隙空间，并激活了额外的流动路径。这种进入压力的重新分布改变了羽流形态，使得羽流前沿更加平滑，并修改了横向羽流的扩散，如图2b所示。

图2
图2. 比较了接触角包含对基础情景中毛细进入压力异质性和泄漏CH4羽流传播的影响：(a) 毛细进入压力分布，(b) 泄漏停止后的CH4羽流传播（90天后），以及(c) 190天后的CH4羽流传播。CA和NCA分别代表包含和不包含接触角的情景。

高分辨率图像
下载MS PowerPoint幻灯片

图2b和c的比较进一步显示，一旦泄漏源关闭，横向CH4迁移变得最小。此时，羽流演化主要由向上迁移驱动，由浮力或通过毛细捕获固定。在这个阶段，由异质进入压力塑造的毛细力与CH4和地下水之间的密度差产生的重力一起，超过了原本会促进水平扩散的粘性力。总体而言，这些结果表明接触角变异性削弱了有效的毛细屏障，重新分配了甲烷流量，并从根本上改变了迁移路径，对甲烷泄漏行为有重要影响。

视频S1提供了代表性情景下CH4羽流通过含水层和非饱和带的动画。结果显示，包含接触角变异性产生了额外的流动路径，减少了仅考虑传统毛细异质性时通常观察到的尖锐羽流前沿。与NCA相比，CA还导致更多的甲烷泄漏到大气中。此外，动画展示了注入停止后含水层中的移动CH4相逐渐分裂成不连续的簇，表明了水力连通性的丧失。这些孤立的气体簇随着时间的推移由于毛细捕获而逐渐固定，阻止了进一步的向上迁移，并有助于CH4在地下长期稳定。

如图所示，从MD模拟计算出的接触角本质上包含一定程度的不确定性。与实验测量结果的比较表明，MD得出的接触角落在实验室结果范围内；然而，某些未考虑的因素仍可能影响润湿行为。为了考虑这些潜在影响，通过将每种相的接触角变化±10°进行了敏感性分析。为了捕捉可能的结果范围，对代表最大和最小相特异性接触角对比的两种极端情景进行了连续介质尺度模拟。这些情景用于评估接触角不确定性对CH4泄漏预测的影响。更多细节见支持信息，第S4.2节。

图3展示了25次随机实现中CH4在含水层、非饱和带和大气中的时间演变的不确定性，比较了NCA和CA情况。这些相配置的变化导致了羽流迁移和捕获行为的差异，这些差异在图中有所体现。在含水层和非饱和带中，CH4被分为移动相、断裂捕获相和溶解相。在模拟器中，溶解的CH4是根据亨利定律计算得出的，而断裂捕获的甲烷则使用Land捕获模型进行量化（详见支持信息，第S3节）。然后通过从总CH4中减去溶解的和剩余的捕获CH4来获得每个网格块中的移动CH4量。由于非饱和带中的CH4量及其相关变异性相对较小，这里仅报告了总值（有关非饱和带中捕获CH4的详细信息，请参见支持信息，第S4.3节）。在大气中，仅显示总CH4，因为该区域中的CH4仅以气态存在。图中每一行对应于基础情景下的特定敏感性情景，例如渗透率增加、盐度升高、CH4泄漏增加或CG相的体积比例增加，而每一列显示了不同部分或相中CH4的时间分布。所得到的CH4百分比表示每个相或部分中CH4摩尔数相对于每次模拟中总泄漏CH4的比例。变异性通过25次实现中的标准差来量化，反映了地质异质性对羽流传播和相分配的影响。图中，实线表示平均值，阴影区域表示平均值上下一个标准差的范围，灰色阴影区域突出显示了CA和NCA情况之间的重叠区域，说明了它们的结果在哪里收敛。

图3
图3. 25次随机实现中泄漏CH4体积百分比的变异性。模拟在含水层、非饱和带和大气中进行，比较了包含和不包含接触角的情况。这些相配置的变化导致了羽流迁移和捕获行为的差异，这些差异在图中有所体现。在含水层和非饱和带中，CH4被分为移动相、断裂捕获相和溶解相。在模拟器中，溶解的CH4是根据亨利定律计算得出的，而断裂捕获的甲烷则使用Land捕获模型进行量化（详见支持信息，第S3节）。然后通过从总CH4中减去溶解的和剩余的捕获CH4来获得每个网格块中的移动CH4量。由于非饱和带中的CH4量及其相关变异性相对较小，这里仅报告了总值（有关非饱和带中捕获CH4的详细信息，请参见支持信息，第S4.3节）。在大气中，仅显示总CH4，因为该区域中的CH4仅以气态存在。图中每一行对应于基础情景下的特定敏感性情景，如渗透率增加、盐度升高、CH4泄漏增加或CG相的体积比例增加，而每一列显示了不同部分或相中CH4的时间分布。所得到的CH4百分比表示每次模拟中每个相或部分中CH4摩尔数相对于总泄漏CH4的比例。变异性通过25次实现中的标准差来量化，反映了地质异质性对羽流传播和相分配的影响。图中，实线表示平均值，阴影区域表示平均值上下一个标准差的范围，灰色阴影区域突出显示了CA和NCA情况之间的重叠区域。

图3
图3. 25次随机实现中矿物异质性下泄漏CH4体积百分比的变异性。模拟在含水层、非饱和带和大气中进行，持续190天，比较了包含和不包含接触角的情况。情景包括(a) 基础情景，(b) 渗透率增加（4 × 基础），(c) 泄漏率增加（2.5 × 基础），(d) 含水层盐度增加（从0增加到约12,000 ppm），以及(e) 体积比例增加（粗粒（CG）体积比例从40%增加到60%）。在含水层中，结果显示了移动CH4气体、断裂捕获CH4和溶解CH4；在非饱和带和大气中，仅显示总CH4气体。缩写：CA = 在毛细压力缩放中包含接触角；NCA = 不包含接触角。

高分辨率图像
下载MS PowerPoint幻灯片

在含水层中，移动CH4浓度在泄漏阶段增加，然后在泄漏停止后随着CH4被捕获、溶解或向上迁移而减少。在各相中，溶解CH4的变异性最小（约10%）。包含接触角异质性使含水层中的CH4减少了多达10%，而大气泄漏增加了多达20%。相比之下，非饱和带中的CH4仅略有增加，CA和NCA平均值之间的差异可以忽略不计。

在更高的泄漏率下，CA和NCA情景之间的差异减小，实现之间的变异性也减小。这种减少是因为以粘性为主的流动抑制了毛细效应，导致CH4无论润湿性假设如何都能更好地进入大气。虽然图3使用标准差对称地展示了变异性，但这种方法可能掩盖了实现之间的全部结果范围。为了更明确地捕捉CH4分布的扩散，我们使用箱形图（图4）评估了190天后的结果，这提供了每个部分中情景之间变异性的更清晰表示。每个箱子代表四分位数范围（第25–75百分位数），中心线表示中位数。须状线显示了非异常数据的分布，单个标记表示异常值。这种统计表示突出了有无接触角（CA与NCA）效应的情景之间的中心趋势和变异性。

图4
图4. 190天后含水层、非饱和带和大气中泄漏CH4体积百分比的比较，当CH4羽流在不同情景中稳定时：基础情景、渗透率增加、泄漏率增加和体积比例增加。结果显示了包含和不包含接触角的情况。(a) 含水层中的移动CH4气体，(b) 含水层中的断裂捕获CH4，(c) 含水层中的溶解CH4；(d) 非饱和带中的总CH4气体；(e) 大气中的总CH4气体。CA = 在毛细压力缩放中包含接触角；NCA = 不包含接触角。

高分辨率图像
下载MS PowerPoint幻灯片

图4的比较进一步显示，一旦泄漏源关闭，横向CH4迁移变得最小。此时，羽流演化主要由向上迁移驱动，由浮力或通过毛细捕获固定。在这个阶段，由异质进入压力塑造的毛细力与CH4和地下水之间的密度差产生的重力一起，超过了原本会促进水平扩散的粘性力。总体而言，这些结果表明接触角变异性削弱了有效的毛细屏障，重新分配了甲烷流量，并从根本上改变了迁移路径，对甲烷泄漏行为有重要影响。

视频S1提供了代表性情景下CH4羽流通过含水层和非饱和带的动画。结果显示，包含接触角变异性产生了额外的流动路径，减少了通常仅考虑传统毛细异质性时观察到的尖锐羽流前沿。CA还导致比NCA更多的甲烷泄漏到大气中。此外，动画展示了注入停止后含水层中的移动CH4相逐渐分裂成不连续的簇，表明了水力连通性的丧失。这些孤立的气体簇随着时间的推移由于毛细捕获而逐渐固定，阻止了进一步的向上迁移，并有助于CH4在地下长期稳定。

如所示，从MD模拟计算出的接触角本质上包含一定程度的不确定性。与实验测量结果的比较表明，MD得出的接触角落在实验室结果范围内；然而，某些未考虑的因素仍可能影响润湿行为。为了考虑这些潜在影响，通过将每种相的接触角变化±10°进行了敏感性分析。为了捕捉可能的结果范围，对代表相特异性接触角最大和最小对比的两种极端情景进行了连续介质尺度模拟。这些情景用于评估接触角不确定性对CH4泄漏预测的影响。更多细节见支持信息，第S4.2节。

图3展示了25次随机实现中CH4在含水层、非饱和带和大气中的时间演变的不确定性，比较了NCA和CA情况。这些相配置的变化导致了羽流迁移和捕获行为的差异，这些差异在图中有所体现。在含水层和非饱和带中，CH4被分为移动相、断裂捕获相和溶解相。在模拟器中，溶解的CH4是根据亨利定律计算得出的，而断裂捕获的甲烷则使用Land捕获模型进行量化（详见支持信息，第S3节）。然后通过从总CH4中减去溶解的和剩余的捕获CH4来获得每个网格块中的移动CH4量。由于非饱和带中的CH4量及其相关变异性相对较小，这里仅报告了总值（有关非饱和带中捕获CH4的详细信息，请参见支持信息，第S4.3节）。在大气中，仅显示总CH4，因为该区域中的CH4仅以气态存在。图中，每一行对应于基础情景下的特定敏感性情景，如渗透率增加、盐度升高、CH4泄漏增加或CG相的体积比例增加，而每一列显示了不同部分或相中CH4的时间分布。所得到的CH4百分比表示每次模拟中每个相或部分中CH4摩尔数相对于总泄漏CH4的比例。变异性通过25次实现中的标准差来量化，反映了地质异质性对羽流传播和相分配的影响。图中，实线表示平均值，阴影区域表示平均值上下一个标准差的范围，灰色阴影区域突出显示了CA和NCA情况之间的重叠区域。

图3
图3. 25次随机实现中泄漏CH4体积百分比的变异性。模拟在含水层、非饱和带和大气中进行，持续190天，比较了包含和不包含接触角的情况。情景包括(a) 基础情景，(b) 渗透率增加（4 × 基础），(c) 泄漏率增加（2.5 × 基础），(d) 含水层盐度增加（从0增加到约12,000 ppm），以及(e) 体积比例增加（粗粒（CG）体积比例从40%增加到60%）。在含水层中，结果显示了移动CH4气体、断裂捕获CH4和溶解CH4；在非饱和带和大气中，仅显示总CH4气体。缩写：CA = 在毛细压力缩放中包含接触角；NCA = 不包含接触角。

高分辨率图像
下载MS PowerPoint幻灯片

在含水层中，移动CH4浓度在泄漏阶段增加，然后在泄漏停止后随着CH4被捕获、溶解或向上迁移而减少。在各相中，溶解CH4的变异性最小（约10%）。包含接触角异质性使含水层中的CH4减少了多达10%，而大气泄漏增加了多达20%。相比之下，非饱和带中的CH4仅略有增加，CA和NCA平均值之间的差异可以忽略不计。

在更高的泄漏率下，CA和NCA情景之间的差异减小，实现之间的变异性也减小。这种减少是因为以粘性为主的流动抑制了毛细效应，导致CH4无论润湿性假设如何都能更好地进入大气。虽然图3使用标准差对称地展示了变异性，但这种方法可能掩盖了实现之间的全部结果范围。为了更明确地捕捉CH4分布的扩散，我们使用箱形图（图4）评估了190天后的结果，这提供了每个部分中情景之间变异性的更清晰表示。每个箱子代表四分位数范围（第25–75百分位数），中心线表示中位数。须状线显示了非异常数据的分布，单个标记表示异常值。这种统计表示突出了有无接触角（CA与NCA）效应的情景之间的中心趋势和变异性。

图4
图4. 190天后含水层、非饱和带和大气中泄漏CH4体积百分比的比较，当CH4羽流在不同情景中稳定时：基础情景、渗透率增加、泄漏率增加和体积比例增加。结果显示了包含和不包含接触角的情况。(a) 含水层中的移动CH4气体，(b) 含水层中的断裂捕获CH4，(c) 含水层中的溶解CH4；(d) 非饱和带中的总CH4气体；(e) 大气中的总CH4气体。CA = 在毛细压力缩放中包含接触角；NCA = 不包含接触角。

高分辨率图像
下载MS PowerPoint幻灯片

CA和NCA情况之间的四分位数范围比较表明，包含接触角增强了与矿物异质性相关的变异性，从而增加了随机实现的影响。在所有情景中，渗透率增加时移动CH4的变异性最大，其中四分位数范围最宽，CA和NCA之间的中位数差异接近10%。相比之下，涉及泄漏率或体积比例增加的情景在变异性和中间值方面显示出最小的差异。

对于断裂捕获CH4，当包含接触角时，四分位数范围通常会增加，除了渗透率增加的情况。在中位数方面，基础情景和盐度增加情景中的CA–NCA差异最大（约20%），尽管绝对值有所不同。这表明含水层的盐度会随着环境条件的变化而变化，从而影响CH4在不同相态之间的分配。值得注意的是，尽管存在这些变化，但在被捕获的相态中，CA（接触角）和NCA（非接触角）之间的差异仍然保持一致，这表明必须考虑接触角效应才能准确量化CH4的捕获情况。溶解态CH4的趋势与被捕获的CH4相似。在非饱和带中，CA情景下的CH4迁移分数始终略高于NCA情景，中位数差异低于3%。唯一的例外是泄漏率增加的情景，在这种情况下，中位数和变异性几乎相同。该区域的CH4百分比在CA情景下为9-12%，反映了由于沉积物松散而缺乏毛细管捕获作用，以及非饱和带在含水层和大气之间的过渡作用。

大气泄漏也是另一个关键因素。正如预期的那样，泄漏率和体积增加的情景会导致更多的CH4释放到大气中，CA情景下的中位数约为75%，NCA情景下为70%。这些发现强调了将接触角变异性纳入建模工作的重要性，因为它会改变CH4在不同相态之间的分配，并直接影响进入大气的排放量。

4. 局限性、意义和未来研究
本研究开发的建模框架旨在考察非均匀润湿性和毛细管异质性对CH4迁移和捕获的相对影响。因此，应对结果中的一些假设和局限性有所认识。首先，模拟基于由达西定律控制的连续介质多相流，并假设每个网格块内都处于局部毛细管平衡状态。诸如膜流、团簇动态和动态接触角行为等孔隙尺度过程并未得到明确解决。润湿性是通过分子动力学（MD）模拟得出的静态接触角来表示的。尽管这些值与实验室测量结果一致，但未考虑由地球化学反应或生物活动引起的潜在时间依赖性润湿性变化。其次，毛细管压力-饱和度和相对渗透率-饱和度关系是使用依赖于经验参数化的本构模型来定义的。相对渗透率滞后现象被包括在内以捕捉排水和吸水行为；然而，由于缺乏相关条件下CH4-水系统的实验滞后数据，因此采用了已建立的经验模型。因此，对残余捕获量和截止值的估计应被视为一般趋势的指示，而不是精确的定量预测。第三，气体和水相之间的质量传递是在平衡假设下处理的，而忽略了反应性传输、微生物对CH4的消耗以及地质力学变形。尽管这些过程在特定现场条件下可能很重要，但在这里被排除在外，以便单独研究润湿性驱动的毛细管异质性效应。在这些限制条件下，结果提供了关于非均匀润湿性如何改变甲烷迁移性、影响毛细管捕获效率以及改变泄漏路径的机制性见解。因此，该框架最好被视为一种基于敏感性和过程导向的工具，它补充了实验室和现场观测，并为未来的更全面建模工作提供了信息。

我们的结果为之前未充分探索的非均匀润湿性在产生毛细管异质性方面的作用提供了概念验证，扩展了早期基于现场校准和模拟的研究，这些研究展示了毛细管异质性对浅层含水层中CH4羽流几何形状和命运的影响。通过将润湿性作为一个独立机制进行隔离，这项工作增强了我们对甲烷迁移和捕获过程的基本理解，并有助于指导未来旨在直接评估润湿性驱动效应的实验、现场和监测研究。同时，简化或忽略毛细管异质性会导致大气中CH4的过度估计以及含水层中可迁移CH4的低估。我们证明了非均匀润湿性（通过矿物学和盐度变化的空间可变水-CH4接触角来表达）通过直接改变局部毛细管入口压力，即使渗透率和孔隙率固定，也会产生独立且可量化的效应。纳入接触角变异性（CA）后，含水层中可迁移的CH4和被捕获的CH4减少了多达10%，而大气泄漏增加了多达20%。在高泄漏率下，CA与NCA之间的差异减小，这与向粘性主导流动的转变一致。

我们的发现对泄漏风险评估具有重要的实际意义。传统模型假设润湿性恒定（NCA），在矿物学和地下水盐度共变的条件下（如富含石英的沿海或蒸发岩含水层），会低估大气泄漏量。这种偏差是因为盐度对石英的接触角影响远大于高岭石，从而降低了以石英为主的含盐系统中的毛细管入口压力，促进了气体更早的突破和更高的泄漏通量。相反，降低地下水盐度可以减小接触角并提高毛细管入口压力，从而增强气体捕获并限制泄漏。为了量化这些动态，应使用基于接触角的Leverett J函数（方程1）来替代仅考虑孔隙率和渗透率的传统缩放方法，特别是在矿物学对比明显或盐度变化的混合相含水层中。这种方法补充了Ershadnia等人（2021年）的研究结果，他们强调了米级相结构在控制羽流迁移中的主导作用。他们表明，非均匀润湿性在相同的尺度上起作用，因此必须对其进行联合表征、绘制和建模，以改进对地下气体传输和泄漏风险的预测。

最近在微CT成像方面的进展使得可以直接从孔隙尺度图像中提取接触角场，证明润湿性本质上是空间异质的，而不是单一的均匀值。基于这一证据，一种在野外应用CT结果的实际方法是使用CT数据作为相级校准器。在这种工作流程中，可以对每个相的代表性岩芯样本进行成像，并计算其接触角分布。然后可以将得到的统计结果作为整个储层模型中的相条件接触角场进行映射。我们基于接触角的J函数旨在直接纳入这些场，从而将孔隙尺度异质性传播到连续介质尺度流动模拟中。当CT数据不可用时，可以使用MD模拟计算的接触角作为物理上可靠的替代方案，随着更多有限的CT测量数据的可用性，可以进一步改进这一方法。这种混合策略规避了直接从CT到野外的放大挑战，将润湿性估计与测量的岩石纹理联系起来，并自然地通过集合采样促进了不确定性量化，同时考虑了渗透率对比、泄漏率和盐度的变化。

在这项研究中，我们专注于清洁的石英和高岭石表面，并假设接触角是静态的；然而，天然沉积物通常表现出额外的复杂性，如表面粗糙度（凹凸不平）、有机涂层、氧化铁沉淀物和接触角滞后（前进角与后退角），所有这些都会显著影响润湿性和羽流迁移。未来的实验和建模工作应考虑这些因素的影响，以更好地模拟未受限浅层含水层中的实际泄漏情况。此外，我们的模拟中盐度保持不变，而在自然界中，由于边界条件的变化，地下水化学成分会随时间演变，这反过来又会影响接触角和毛细管行为。因此，连续介质尺度上接触角变异性和泄漏CH4羽流在空间和时间上的测量对于验证和细化预测模型以及捕捉CH4泄漏动态至关重要。这些数据结合改进的模型将提高对地下水质量和逃逸CH4排放风险评估的准确性。

我们的发现的意义不仅限于浅层含水层。在深层地质系统中，非均匀润湿性仍然是一个关键未知因素，包括CO2和H2储存或有毒化学物质（如NAPL）的储存，其中高迁移率使得毛细管异质性成为控制羽流传播和短期捕获的关键因素。此外，注入流可能含有或与残留CH4相互作用，这是由于共注入、遗留碳氢化合物或捕获和传输过程中的杂质携带所致。因此，在评估更大空间尺度上的流动和捕获行为时，也应考虑涉及CH4的相互作用，特别是考虑到先前的现场校准证据表明毛细管异质性对CH4羽流迁移和泄漏具有第一级控制作用。通过整合实验室和现场数据与多尺度建模研究，将改善对浅层和深层地下流体迁移的预测能力，并支持更安全和更有效的能源和气候缓解技术的发展。