介质润湿性和非均质性对重非水相运移影响孔隙尺度机制研究

彭锦潇，蒋建国，吴吉春

南京大学 地球科学与工程学院，表生地球化学教育部重点实验室，南京 210023

重非水相（DNAPLs）是一大类有机污染物，例如全氯乙烯（PCE）、煤焦油，其密度大于水且不与水互溶。它们很容易到达含水层底板，形成污染羽从而加大修复难度（Rathfelder et al., 2001；Lorenzo et al., 2015）。大量研究表明DNAPL 的运移主要受介质润湿性，介质结构，DNAPL 质量和粘度影响（Ryder et al., 2008；Molnar et al., 2010；Koch et al., 2015；Zhang et al., 2012）。最近，程洲等利用砂箱实验表明当石英砂多孔介质由水相润湿变为油相润湿后，DNAPL 的下落速度减小，甚至可能会悬停在介质中。实验中，将PCE 注入扁平的砂箱，背景液是水或表面活性剂Tween 80 溶液，增加Tween 80 的浓度能够使介质接触角增大从而改变润湿性。结果表明接触角增加PCE 下落速度减小且残余饱和度升高。尽管文章对宏观机理进行了简单解释，但是缺少对迁移细节的详细分析（Cheng et al., 2016）。

DNAPL 在地下水饱和带迁移过程是典型的多相流现象。由于多孔介质多相流在许多领域，如污染物修复、石油储藏、碳截存等有重要作用（Ran et al., 2018），许多研究使用连续介质模型对饱和多孔介质的DNAPL 位移进行模拟（Zhang et al., 2001；Nayagum et al., 2004；Kokkinaki et al.,2013）。但这些模拟难以刻画流动时可能出现的大量指流，指流相对于稳定流动有更快的迁移速度，对DNAPL 在含水层的最终分布有重要影响（Chatzis et al., 1983；Lee et al., 2012；Mason et al.,2013）。因此前人提出许多微观模拟方法（Ran et al., 2015；Trojer et al., 2015），格子玻尔兹曼方法（LBM）能够轻松刻画复杂孔隙边界的物理流动而被广泛采用（Pan et al., 2007）；孔隙网络模型（PNM）通过将复杂的三维孔隙简化为孔喉和孔体，大大降低了计算量，但精度相对较低（Qiang et al., 2016；Bondino et al., 2009）；计算流体力学（CFD）方法，例如水平集和相场等能够更加准确捕捉界面的压力和形状。由于本文侧重机理研究，而相场方法收敛性好对计算能力要求相对较低，故选择相场方法（Amiri et al., 2013）。

笔者选择在较小的空间尺度上进行模拟，以揭示流体在颗粒之间的流动细节。尽管前人在孔隙尺度上进行了一些研究，但他们侧重于单个界面的驱替过程。DNAPL 在重力作用下落具有自身特点，这个过程受到界面处压力差驱动，与前后界面高度差相关，前后界面均会影响下落过程，从模拟结果来看也有所印证。因此单界面模型不足以简单明了解释实验现象，笔者将提出一个双界面模型。

本文首先介绍模型和模拟方法，然后给出模拟结果并与实验结果对比。随后结合DNAPL 前后界面对其流动细节进行分析。同时绘制压力分布图，以直观地解释润湿性对DNAPL 迁移的影响。此外通过具有完全相同尺寸颗粒的均质模型佐证分析的正确性，也表明润湿性的影响取决于多孔介质的非均质性。最后，总结文章主要结论。

1 模型和模拟方法

1.1 模型介绍

将平板砂箱简化为二维模型，如图1 所示，其中圆盘代表石英颗粒，实验中砂箱为40～50 目，孔径为0.27～0.38 mm （Cheng et al., 2016）。模拟中为了使得邦德数大于1，若使用实验孔径的圆盘将使得非水相高度非常高从而无法计算，所以笔者将颗粒半径放大处理。

微观模拟时，构造二维孔隙通常有两种方式：将圆盘圆心以正三角方式均匀摆放，然后改变圆盘半径长度（Ran et al., 2018）；或使用两组一大一小的圆盘，大圆盘仍然均匀摆放，小圆盘随机插入大圆盘孔隙（Di Palma et al., 2017）。第一种方式收敛性要求高，且无法突出小孔隙和大孔隙特点，结合模拟结果发现残留几乎都发生在有小圆盘的小孔隙中，故笔者选择第二种方式模拟。如图1a，白色圆盘代表固体颗粒，在y轴（重力方向）上交错排列。如图1b，定义大圆盘半径为R，相邻圆心距离为d=R+Lt，其中Lt为大孔喉特征直径。中心随机插入半径为Rmin的小圆盘。这样小孔喉特征直径为Ltmin=d/ √ 3 -R-Rmin，模型中两孔喉直径比值为Lt/Ltmin= 1.5。

图 1 二维孔隙模型Fig. 1 2D pore model

在实验中，DNAPL 下落分为两个阶段，分别为4 min 的连续注入和后续自由下落。由于注入时间相比于下落时间短很多，笔者更关注后阶段的行为，所以直接模拟下落过程。如图1 中所示，圆盘空隙的蓝色部分表示非水相的初始位置，灰色部分充满水。DNAPL 的前界面初始高度为45.5 mm，后界面高63.6 mm。上边界为水相入口，压强为0；下边界为出口，压强为纯水压强653.7 Pa。左右两侧边界均为对称边界，多孔颗粒壁设为无滑移边界和恒定水相接触角θ。通过简化的推导可知：在完全均质的孔隙中，如果接触角互补，那么前后界面的表面张力的合力是相等的。为排除接触角本身大小变化带来的差异，表明润湿性的影响机制。对于水相润湿，笔者选择与实验中水相和石英接触角一致为40°，而对于油相润湿，控制接触角与40°互补为140°（与实验中124°均为油相润湿）。需要说明的是，这两个接触角下的模拟结果均未观察到薄膜流动，事实上，只有当接触角小于20°或大于160°时，薄膜流动才显著出现。但由于和实验参数偏差较大所以没有考虑。模拟除了接触角，流体和介质其他参数均相同，根据程洲实验设置（Cheng et al., 2016；见表1）。

1.2 模拟方法

相场方法通过耦合Cahn-Hilliard 方程，Navier-Stokes 方程和连续性方程来求解多相流问题。首先定义各流体体积分数V与相场变量φ的关系：Vw= (1+φ)/2 ，Vnw= (1-φ)/2 。其中，φ=±1 代表流体内部；- 1＜φ＜1代表两相界面，同时使用φ定义边界的其他物理属性：

表 1 模型参数Table 1 Model parameters

其中ϑ表示每个相的属性，例如粘度μ 和密度ρ。将不可压缩的Navier-Stokes方程和连续性方程耦合相场方程，得到方程组如下：

式中Ψ是相场助变量，λ是混合能量密度， ∈是界面厚度参数，∈设置为最小剖分长度的一半，γ是迁移率，这个参数必须适中，以保持界面厚度的同时防止过度抑制对流，模型中设置为γ=ε2(Amiri et al., 2013)。G是化学式，，而表面张力系数可以关联λ和∈得到，

模型中固体颗粒使用无滑移边界条件，方程组如下：

模型中两侧边界使用对称边界，边界垂直方向流速为0，且没有剪切力损失，方程组如下：

由于模型构建在垂直方向上，所以bond数也很重要，其定义为：

其中rt表示特征孔喉半径，笔者设置DNAPL的初始高度满足邦德数大于1（约为1.3），使得DNAPL可以自发迁移。

模型使用三角形网格进行剖分，为保证多相界面足够狭长，剖分尺寸设置为Lmin/12，其中Ltmin表示最小孔喉长度（Amiri et al., 2013）。同时使用自适应对网格进行局部加密，加密使用的误差指示函数为相场变量φ梯度的L2范数

2 模拟结果与讨论

2.1 模拟结果

图2 分别展示了介质为不同润湿性时，DNAPL 下落1 s 后的分布图，其中左图为水相润湿，右图为油相润湿。圆盘间隙中红色部分表示DNAPL，其余部分表示水。DNAPL 的质心位置用蓝实心圆标注。

DNAPL 的下落速度由质心位移和时间关系计算得来。如图3 所示，在开始的0.18 s 内,不同润湿性条件下落速度并没有差异。当速度达到26 mm/s 后，曲线出现分离。水相润湿时速度继续增加，而油相润湿速度开始逐渐减小。这是因为在流动前期，DNAPL 形状相同，且接触角互补，前后界面贡献的表面张力总和也相等，所以下落速度一致。但0.18 s 后，水相润湿条件中，DNAPL 前界面开始出现指流，油相润湿DNAPL 后界面开始出现残留，导致下落速度发生变化。因此，分析DNAPL 的形状是润湿性对运移影响关键因素。

2.2 速度分析

DNAPL 的下落由压力差驱动，压力差与其高度成正比，微观模型可以准确地模拟压力变化。观察图2 可见，DNAPL 的形状差异明显，水相润湿时DNAPL 前界面发育指流，而油相润湿DNAPL后界面出现大量残留。如图3 所示，指流会增大DNAPL 的高度从而增大DNAPL 前后界面的压力差，这使得水相润湿下落速度增高。而残留意味着压力差减小从而下落速度降低。

为研究形状差异机制，笔者截取不同时刻的两相界面位置细节，如图4 所示。其中图4a 表示非湿润相驱替湿润相，4b 表示湿润相驱替非湿润相。图中彩色弯液面表示不同时刻的两相界面，例如从蓝至红，分别表示0.01 s、0.02 s、0.03 s和0.04 s 时刻的两相界面，黑色箭头表示结束时刻的流场。

图4a 中的弯液面a总是朝向防御流体，并在0.03 s 后加速移至c，然后触碰到下一个颗粒形成新的弯液面d1和d2。新形成的弯液面d2毛细管压力与上部小孔隙弯液面平衡，阻止其下落从而形成指流和残留。因此，当非湿润相置换湿润相时，侵入相绕过被侵入相并形成指流，同时被侵入相阻塞而形成残留。

图 2 DNAPL下落1s后分布图Fig. 2 DNAPL distribution after falling for 1 second

图 3 DNAPL中心下落速度Fig. 3 The falling speed of the center of DNAPL

图 4 不同时刻的两相边界变化图Fig. 4 Two-phase interface variation diagram at different times

而当湿润相换置换非湿润相时，图4b 所示，边界则趋于平滑。若彩色弯液面在相邻两孔中速度差异不大，先通过孔喉的驱替界面会在毛细管压力为0 处附近停留，等待相邻孔隙的驱替界面。这是由于弯液面在此处水平很难触碰到下一个颗粒，同时此处作为拉力的毛细管压力很小。此时另一侧弯液面a快速通过孔喉达到d，随后e1和e2汇合形成一个新的液面f之后再触碰到下一个颗粒。这种交替前进，协同填充机制抑制了指流的发育和被侵入相的残留。协同填充的能力主要和孔隙度有关，孔隙度决定临界接触角θc，当低于θc时形成稳定位移。根据胡冉等人的补充相图，在低速且接触角接近0 时，流体形成薄膜流动沿着孔隙壁前进，增加了单孔隙snap-off 模式的捕获使残留量增大，但笔者没有在40°的接触角下观察到薄膜流（Ran et al., 2018）。

此次研究的模拟结果与Holtzman 的论文的结论一致，该论文表明，增加侵入流体的接触角有利于协同填充，从而形成稳定驱替（Holtzman et al., 2015）。

图 5 不同润湿性的压力分布图Fig. 5 The pressure distribution of fluids

可以结合双界面模型来解释DNAPL 下落速度差异原因。在水相润湿条件下，前界面按照图4a的模式驱替，发育指流；后界面按照图4b的模式驱替，协同填充。两个界面都对下落速度有贡献，导致下落速度增加。相反在油相润湿中，前界面对应图4b模式，后界面则对应图4a模式，几乎没有指流增强下落速度且残留很多，最终下落速度减缓。

双界面模型可以清楚的解释程洲等的实验结果，在他们的实验中，油相润湿中的PCE 保持着平整的前进形状且形成大量残留，所以下落速度相比于水相润湿中低很多。推测宏观试验效果更显著可能是因为非均质性更强（Cheng et al.,2016），笔者尝试了从模拟角度证明推测，但是非均质性更强意味着收敛难度更大，将在后续的研究中探索更有效的剖分方式。

2.3 残留分析

残余饱和度是污染物修复的重要指标。模拟结果与实验一致，残余饱和度也受介质润湿性影响（Cheng et al., 2016）。如图2 和表2 所示，非均匀介质中，油相润湿时残留量比水相润湿大得多。残余饱和度主要由后界面驱替特征决定，在油相湿润条件下，水相按照图4a 的方式置换DNAPL形成指流。这些指流绕过小孔侵入DNAPL，包围小孔中的DNAPL 形成隔离，最终残留在多孔介质中。相反水相润湿时，根据图4b 中的模式，后界面平滑，残留很少。这样笔者提供了一种直观的方式来解释残余饱和度差异原因。

同时，下落初期垂直方向的压力分布图也能对残余饱和度差异原因给出启发性解释。如图5所示，红线表示油相润湿，蓝线表示水相润湿。油相润湿中，后方边界的液滴内部压力小于周围水相压力，导致水相容易进入负压区域而形成残留。相反，在水相润湿中，后界面的内部压力大于周围水相压力，导致残留较少。

表 2 模拟结束时下落速度和残余饱和度表Table 2 The statistics of falling speed and residual saturation

图 6 完全均匀介质中DNAPL下落速度图Fig. 6 The falling speed of DNAPL in homogeneous medium

2.4 均质模型

上述分析的合理性可以用均匀介质佐证。DNAPL 的下落速度受前界面的指流和后界面的残留影响很大，如图2 所示，指流在水相润湿中大孔中形成，而残留则发生在油相润湿中的小孔中。如果孔径大小一致，指流很难发育，DNAPL 也难以残留，使得接触角改变对迁移的影响会小很多。为了验证这一假设，通过去除模型（图1）中的所有小圆盘来构造相同孔径的均匀介质，其余参数与上述模型一致。图6 所示，不同润湿性下DNAPL 下落速度基本相等，残余饱和度分别为0.00和0.01，具体数据如表2 所示。因此，润湿性的影响是基于多孔介质非均质性强弱的。笔者的模拟时间只有1 s，如果时间更长，下降速度的差异可能还是会显现。

3 结论

DNAPL 在地下水中的迁移很大程度受多孔介质润湿性的影响。程洲等实验表明：当水和固体颗粒之间的接触角变大时，PCE 的下落速度降低、残余饱和度增加（Cheng et al., 2016）。通常微观多相流理论研究仅考虑单个界面的流动，这不足以解释实验现象。在本文中，使用双界面模型模拟了DNAPL 在孔隙尺度的重力下落过程，在孔隙尺度重现与实验结果相似的现象。通过对DNAPL 微观流动的模拟和定性分析，得到如下结果：

（1）水相润湿时，DNAPL 前界面发育指流，后界面平整，均有利于下落。油相润湿情况相反，DNAPL 前界面平整，后界面DNAPL 残留，均抑制下落。

（2）残余饱和度主要由后界面性质决定。水相润湿后界面平整，DNAPL 残留较少。而油相润湿后界面易发育指流，DNAPL 易在小孔隙处被截断，形成残留，因此残余饱和度较大。

（3）此外，构建相同孔径的均质多孔介质，在模拟时间尺度内，润湿性的改变对DNAPL 下落速度和残余饱和度影响不大。因此，润湿性对DNAPL 迁移的影响需要通过多孔介质的非均质性起作用。