水动力和介质非均质性影响下DNAPLs运移数值模拟

韩智颖，王锦国，崔孜铭，杨蕴

(河海大学地球科学与工程学院，南京 211100)

随着经济的飞速发展，社会对石油化工产品需求量日益增加，导致每年大量的挥发性液态有机化合物从地表泄漏并渗透进入地下，对土壤环境和地下水资源造成严重污染。这些液态有机化合物一般与水互不相溶，而是以非水相液体(Non-aqueous Phase Liquids，NAPLs)的方式与水相、气相共存于土壤-地下水系统中组成复杂的多相流运移问题。根据NAPLs密度大小与水相密度比较，通常可分为轻非水相液体(Light Non-aqueous Phase Liquids，LNAPLs)和重非水相液体(Dense Non-aqueous Phase Liquids，DNAPLs)[1]。

相较于LNAPLs而言，DNAPLs具有更低的粘滞性和界面张力，其运移规律更为复杂多变；在均匀介质中，DNPALs在非饱和层中的迁移过程与LNAPLs相似，主要受重力和毛细力的共同影响, 同时产生水平向和竖直向的位移；自然条件下重力的影响远远大于毛细力，污染物分布形态总体表现为向下运移；但由于DNAPLs相对于水相密度较大，到达饱和层中后，地下水动力占主导作用，LNAPLs基本沿水平方向运移[2-3]；而DNAPLs则渗过潜水面继续运移，并在地下水动力和重力的协同作用下最终集聚在饱和含水层底部。所以研究DNAPLs在地下介质中的运移分布规律等具有重要的意义和价值。

地下水动力、泄露速率、降雨入渗、多孔介质的非均质性等是影响DNAPLs迁移和空间分布的主要外部因素，其中地下水动力和介质的非均质性影响DNAPL运移规律更为显著[4-9]。Erning等[4]采用TOUGH系列软件TMVOC模块研究了地下水流速变化对DNAPL入渗扩散行为的影响。研究表明，DNAPL体对应用的地下水流速的反应比场景的几何设置更敏感。即使在其位置、尺寸和结构的最低流速下，DNAPL也受到影响。施小清等[7,8]和Yang等[5]研究了介质非均质性、水相饱和度、DNAPL泄露速率等对DNAPL运移分布的影响，结果表明介质非均质性是污染物达到最大饱和度并集聚的主要原因。

非均质性是实际天然土层的空间结构最明显的特征之一，野外的天然地下介质通常存在大孔隙、多裂隙发育、软弱夹层、低渗透性透镜体等结构，有时也夹杂着断裂等地质构造，从而造成地下多孔介质具有高度复杂的非均质性[10]。且由于自然环境的水力梯度一般横跨范围较大且可能有人为因素的影响，例如抽水井附近地区水力梯度一般大于正常水力梯度，因此本文采用不同水力梯度和不同位置和大小的透镜体模拟地下水动力和介质非均质性以及该两种因素协同作用对DNAPLs运移规律和分布特征等影响。

1 多相流运移数学模型

TMVOC是水、土壤气体和挥发性有机化合物(VOCs)多元混合物在非均质多孔介质中三相非等温流动的数值模拟软件，采用FORTRAN语言编写，是劳伦斯伯克利国家实验室开发的TOUGH2通用模拟程序的扩展[11]。TMVOC软件中的多相流计算体系包含水相、NAPL相以及不可压缩气体(NCGs)三部分，多相流在迁移的过程各组分之间会发生扩散并相互转化，并遵循多相、多组分计算体系中流体和热流的质量和能量守衡方程如下：

(1)

式中，q表示质量或热流；q表示汇和源；vn表示单位体积;τn为流体穿过的面积；Mk为组份k在单位体积内的质量；Fk为组份k进入单位体积内的通量；qk为k组份在单位体积内的源汇项；n为单位体积向外的法向量矢量，向内指向vn。

采用积分有限差分方法进行空间和时间的离散化，所有计算网格块的相关主要热力学变量用牛顿-拉夫逊迭代法求解，并假设三相中每相的运动遵循引入相对渗透概念的Darcy定律，如下：

(2)

式中，β表示相位;kγβ为相的相对渗透率(介于0和1之间);pβ=p+pcβ是相中的流体压力;pcβ是参考相(通常是气相)中的压力和毛细管压力的总和(毛细管压力为负)。

每个相在压力和重力的作用下流动，其中包括相之间的相对渗透率和毛细管压力的影响。本文选用相对渗透率Stone模型和毛管压力Parker模型，表达式如下：

相对渗透率Stone模型：

(3)

(4)

(5)

式中，kwr、kgr、knr分别为水相、气相以及NAPL相的相对渗透率;Sw、Sg、Sn分别为水相、气相以及NAPL相的饱和度(Sw+Sg+Sn=1);Swr、Sgr、Snr分别为水相、气相以及NAPL相的残余饱和度;n为拟合参数。

毛管压力Parker模型：

(6)

(7)

(8)

(9)

式中，Sm为液相的残余饱和度;Pcgn为NAPL相-气相之间的毛管压力;Pcgw为水相-NAPL相之间的毛管压力;Swe和Sle分别为水相和液相的有效饱和度;αng、αnw分别为NAPL相-气相、水相-NAPL相之间的进气压力的倒数;n为拟合参数[12]。

2 理想模型概况与基本参数

本文采用TMVOC模拟不同水力梯度和不同位置和大小的透镜体影响下DNAPLs运移规律和分布特征。算例概化如下：理想模型为二维尺度，长(X)为100 m、宽(Y)为1 m和高(Z)为15 m的长方体(图1)，共分为680(17×1×40)个网格，其水平方向和垂向的网格剖分如下(表1、表2)。

图1 概念模型示意图

表1 X、Y方向网格剖分

表2 Z方向网格剖分

模型的顶部设定为具有降雨入渗的大气边界，为通量边界；底部以及前后两个面(XZ剖面)为零通量边界，左右边界(YZ剖面)为定水头边界，左边界水头高于右边界。模型在水-气两相系统中处于平衡状态，水相饱和度为1.0，模型温度恒为20 ℃，不考虑热物理相关过程运算。氯苯作为DNAPLs代表，以点源方式在模型顶部以下0.5 m处均匀注入，泄露速率为 1×10-5kg/s，模拟其在连续注入3 a 内的运移和分布结果。将地质体整体概化为粘土岩相，用SAND表示；大气边界概化为ATOMS岩相；低渗透性层状透镜体模型概化为粘土岩相，用IMP表示。模型的主要参数值见表3。

表3 模型参数的设置

3 结果与分析

3.1 水力梯度

首先考虑不同水力梯度对DNAPL运移分布的影响，水力梯度分别设置为0、0.01和0.02，分别对应小、中、大3种流速，并选取0.5 a、1 a、2 a 三个时间节点XVOCM(液相VOC的总质量分数)Y=0横剖面图作为实验结果展示(图2)。

图2 不同水力梯度下0.5 a、1 a、2 a XVOCM Y=0横剖面图

模拟结果表明：i=0时，氯苯在重力的影响下垂向迁移，污染范围基本呈对称结构；i=0.01时，随着水力梯度的增加，在0～0.5 a内，垂向入渗占据绝对的主导地位，氯苯沿注入方向快速向下运移。这主要是由于增大的水流流速减小了氯苯驱替水相的阻力[15]。随着时间的增加，污染羽受水力梯度的影响逐渐向右运移；在1 a时，污染羽垂向运移的深度和向下游侧运移的长度基本一致；i=0.02时，当接近潜水面时，由于水力梯度较大，氯苯便开始随水流的方向进行弥散。随着时间的增加，氯苯基本向右运移，污染羽向下运移的深度和向右运移的长度约为1∶5[13]。

3.2 非均质性

大部分的学者主要通过数值模拟一般采用渗透率随机场实现介质的非均质性[5]或室内实验一般用具有不同性质或位置大小的夹层或透镜体来实现介质的层状非均质性[16]。本模拟水力梯度设置为0，在饱水带设置不同位置和大小的层状低渗透性透镜体，并选取6月、12月、18月、24月、30月、36月时间节点XVOCM(液相VOC的总质量分数)Y=0横剖面图作为实验结果展示(图3)。为了定量分析透镜体上方饱和度的变化，分别选择网格251、249、254作为3个监测单元，监测透镜体上方饱和度在36个月内的变化(图4(a))。

图3 饱水层透镜体9～24月XVOCM Y=0横剖面图

模拟结果表明：结合图3和图4(a)看出经过6个月左右后，DNAPL运移到透镜体位置，DNAPL开始在透镜体正上方集聚，并形成小型污染池；在24个月后，DNAPL饱和度达到一定时，但并未完全渗透穿过透镜体，便开始向左右两侧水平运移，并在重力影响下，形成明显的污染指。可以看出透镜体对于DNAPL在饱水带的垂向渗透不仅具有明显的阻滞作用，并加速DNAPL水平运移，使其运移速度明显大于垂向速度；绕过透镜体后向下运移，图4(b)下游侧入渗深度明显大于上游侧，并在上下游两侧形成较为明显的污染指。透镜体主要发挥的是阻滞作用，明显影响了DNAPL运移的分布、结构特征[14]。

图4 饱水层透镜体对DNAPL运移速度和分布范围的影响

3.3 不同水力梯度下透镜体对DNAPL运移的影响

通过改变透镜体在模型中的大小和分布位置(图5)，并选取21月、24月、27月3个时间节点XVOCM(液相VOC的总质量分数)Y=0横剖面图作为实验结果展示(图5)。探讨水力梯度和介质非均质性协同作用对DNAPL运移分布规律、结构特征的影响。

图5 不同位置透镜体21～27月XVOCM Y=0横剖面图

模拟结果表明：

(1) 有透镜体存在时，在相同的时间内污染物的运移明显迟缓，27月时由于透镜体的存在，DNAPL污染物比无透镜体时少垂向运移了近7 m，证明透镜体对DNAPL污染物的运移具有明显的阻滞作用。在21月时，DNAPL污染物已经在透镜体上积蓄形成污染池，饱和度逐渐增加；当饱和度增加到一定时，DNAPL开始绕过透镜体垂向运移，形成了明显的污染指；当设置了从左至右水力梯度后，上游污染指运移速度明显大于下游侧，且垂向运移的速度也明显加快。

(2) 当水力梯度和透镜体同时存在时，水力梯度的存在使得介质非均质性影响进一步增强，进一步增大了优势通道的出现概率；造成DNAPL结构几何效应的原因主要有两个方面：一是流动水的输运效应，将DNAPL推向下游方向；二是地下水流速越高，传质速率越高，溶解渗滤路径的速度比DNAPL池快[11]。在无水力梯度的均质介质中 DNAPL污染物分布的范围基本以对称结构分布，当地下水动力存在时增强了地下多孔介质的非均质性，导致 DNAPL迁移路径的空间变异性和不规则性增强产生DNAPL污染物在透镜体的影响下更容易产生蓄积和复杂的绕流，形成形态多变的污染池和污染指交替组合。

4 结论

本文采用TMVOC软件通过改变水力梯度和透镜体位置和大小模拟水力梯度、介质不均匀性以及该两种因素共同对DNAPL污染物运移的具体影响。主要结论归纳如下：

(1) 水力梯度对DNAPL在入渗过程中和入渗后的行为有显著影响。随着水力梯度的增加，其渗流路径以及稀释和置换的DNAPL池逐渐朝水流方向倾斜。对于地下水力梯度较高的现场，重力的影响会被削弱。

(2) 在透镜体的影响下，DNAPL在低渗透性透镜体上方滞留集聚并加强了其水平迁移的速度逐渐形成污染池，污染的范围显著增大；污染池达到一定饱和度时，DNAPL污染物绕过透镜体，朝着上下游两侧形成的优势通道继续入渗。

(3) 水力梯度和透镜体同时存在时，水力梯度的存在使得介质非均质性影响进一步增强，进一步增大了优势通道的出现概率；地下多孔介质的非均质性增强造成 DNAPL运移路径的不规则性增强，使其在饱水含水层产生集聚和绕流，形成了多种形态的污染指和污染池交替组合。