水位波动下包气带透镜体影响LNAPL迁移的数值模拟研究

潘明浩 ，时 健 ，左 锐 ，赵 晓 ，刘嘉蔚 ，薛镇坤 ，王金生 ，胡立堂

（1.北京师范大学水科学研究院，北京 100875；2.地下水污染控制与修复教育部工程研究中心，北京 100875）

轻质非水相液体（Light Non-Aqueous Phase Liquid，LNAPL）是土壤及地下水环境中有机污染问题的重要来源，往往具有持续危害性[1-2]。土壤包气带是LNAPL下渗迁移过程中的重要通道[3]，其含水量、氧化还原条件等随土壤深度变化，使得LNAPL在包气带中的多相流迁移行为变得较为复杂[4-5]。一般认为，LNAPL进入包气带后首先以垂向迁移为主，到达地下水位上方的毛细带后将在毛细力作用下驱替孔隙水，发生侧向迁移扩展[6-8]。LNAPL在包气带中的迁移行为受黏度、密度等流体性质的控制[9-10]以及介质岩性、地下水位波动等外部因素的影响[7,11]，其中非均质性（如透镜体、非均质夹层、裂隙等）[12-13]以及水动力条件[14-15]在很大程度上将影响LNAPL的迁移规律和相态分布特征。

国内外学者广泛采用物理模型试验或数值模拟方法探究包气带内的LNAPL迁移过程。非均质条件下的介质含水率差异以及毛细力作用均被认为是控制LNAPL迁移规律的重要因素。Schroth[16]在通过砂箱实验探究LNAPL在非均质包气带中的迁移过程时发现含水量对LNAPL分布存在显著影响。束善治等[7]、陶佳辉等[17]则通过数值模拟探究不同地层结构中LNAPL的迁移过程，发现介质含水率差异决定了非均质结构对LNAPL运动的阻碍或促进作用。此外，Simantiraki[18]与Wipfler[19]分别运用图像分析与数值模拟手段分析LNAPL迁移过程，二者均认为毛细力是污染物迁移过程的重要控制因素，LNAPL垂向运动受阻可能与界面形成的毛细屏障有关。水位波动是土壤及地下水污染问题中的另一关键因素，对于LNAPL污染物而言，水位波动将显著影响气-水-LNAPL三相之间的相互驱替，进而改变污染物迁移及相态分布过程[20]。罗凌云[11]通过实验探究水位波动条件下柴油进入包气带的迁移行为时，发现LNAPL在水位上升时将会优先驱替多孔介质中的气体，在水位下降时则将优先驱替水相；杨明星[15]结合实验及数值模拟手段对石油组分在水位波动带内分异演化机理进行研究，发现污染物在水位波动条件下受“界面效应”影响，趋向于滞留在较低渗透性的介质中，水位波动将加剧LNAPL吸附、溶解等过程并使其对水的驱替过程减弱。因此可以说，地下水位波动将显著影响LNAPL的迁移过程[21]。

为探究土壤包气带内存在透镜体结构这一复杂环境条件下LNAPL的迁移过程及其相态分布，本文基于北京市通州区张家湾包气带试验基地的场地状况及其地下水位波动背景，采用TOUGH2（transport of unsaturated groundwater and heat）程序的TMVOC模块进行算例研究，构建二维非均质剖面中LNAPL入渗迁移的多相流模型，以透镜体介质岩性、地下水位波动作为共同影响因素，模拟非均质包气带内LNAPL的迁移规律、分布特征及挥发过程，以期为复杂条件下的石油污染场地修复提供理论参考。

1 材料与方法

1.1 模型背景

TOUGH2程序可用于模拟多孔介质中的地下水流、非等温热量运移等多相流过程[22]，TMVOC模块是该程序后续版本中用于模拟NAPLs在饱和-非饱和介质中迁移、VOCs（挥发性有机污染物）溶解与挥发等过程的模块[23]，采用积分有限差分法进行模拟区域的空间离散[24]。

通州区为北京城市副中心，近年来该地区地下水位总体呈上升趋势，但伴有水位的持续波动。在水位波动背景下，其土壤地下水中有机物污染问题愈发受到关注。通州区张家湾包气带试验基地位于永定河与潮白河冲洪积扇中下部的交汇地带，属于第四系沉积地层，其上部地层岩性主要包括细砂、粉细砂、中砂等，场地介质非均质性较强且局部存在厚度为几厘米至几十厘米不等的微观非均质透镜体结构。在此基础上，本文设计包气带内存在透镜体的非均质地层结构，假定LNAPL污染物由近地表泄漏，采用TMOVC模块建立XZ方向的二维剖面数值模型。如图1所示，模型尺寸为90 cm×10 cm×110 cm，其中X方向剖分为18列网格，Z方向剖分为22层网格，网格水平及垂直间距均为5 cm。透镜体区域坐标范围为X=30～60 cm；Z=-70～-75 cm，LNAPL泄漏源坐标范围为X=40～50 cm；Z=0～-5 cm。此外，模型共设置6个观察点用于观察透镜体周围的LNAPL运移状况。

图1 概念模型示意图Fig.1 Schematic diagram of the conceptual model

1.2 模型参数

本文采用中砂介质作为包气带内的主要介质，分别采用细砂及粗砂作为透镜体介质。采用Stone模型作为相对渗透率函数，其中Swr为水相残余饱和度；Snr为NAPL相残余饱和度；Sgr为气相残余饱和度；n为拟合参数。采用van Genuchten模型作为毛细压力函数，其中m=1-1/n；Slr为液相残余饱和度；Sls为饱和液相的饱和度；P0为进气值；Pmax为最大进气压力值。介质相关各参数如表1所示。

表1 介质相关参数[17,25-26]Table 1 Medium-relevant parameters

苯系物（BTEX）是较为常见的LNAPL污染物[27]，本文采用甲苯作为算例中的LNAPL污染物，考虑甲苯的挥发及溶解等相态传质过程，相关参数如表2所示。模型在等温模式下运行，温度恒定为25 ℃。

表2 甲苯相关参数[24]Table 2 Toluene-relevant parameters

1.3 边界条件及初始条件

模型顶部单元格设定为大气边界（Dirichlet边界），垂向范围为Z=0～5 cm，压力值恒定为大气压1.013×105Pa，气体饱和度为0.99。模型底部一层单元格为定压Dirichlet边界，其深度范围Z=-100～-105 cm，初始压力值按照地下水位-100 cm计算并设定为1.017 9×105Pa。模型构建完成后首先计算得到初始水分分布，图2所示分别为中砂介质包气带内存在细砂、粗砂透镜体时模型中的初始水分分布情况。

图2 各透镜体模型中的含水率变化规律Fig.2 Variation pattern of water content in each model

2 结果与讨论

2.1 不同岩性透镜体条件下自由相迁移分布

本节设计包气带内存在细砂、粗砂透镜体两种模型算例，分别以LNAPL饱和度SNAPL、LNAPL在水相中的质量分数x(aq)作为表征自由相、溶解相分布的定量指标，分析水位恒定时LNAPL在包气带内的迁移分布。LNAPL释放时间为0～180 min，释放速率为4.483 mL/min。图3及图4分别代表包气带内存在细砂、粗砂透镜体时LNAPL的迁移过程。

图3 细砂透镜体中不同时刻LNAPL饱和度分布Fig.3 Distribution of LNAPL saturation at different moments in the fine-sand lens model

图4 粗砂透镜体中不同时刻LNAPL饱和度分布Fig.4 Distribution of LNAPL saturation at different moments in the coarse-sand lens model

对于细砂透镜体模型，LNAPL迁移到达透镜体上方时存在明显的蓄积和横向扩展现象，180 min前透镜体两侧SNAPL几乎为0，且透镜体上方（观察点2）及下方（观察点5）SNAPL相对较大，表明部分LNAPL已优先以“穿过透镜体”形式向下迁移。对于粗砂透镜体模型，LNAPL更多以垂向迁移为主，表现为150 min时LNAPL分布深度超过细砂模型约17%，105～180 min时LNAPL在透镜体下方发生横向扩展。污染物释放阶段的模拟结果表明，粗砂透镜体成为LNAPL垂向迁移的“优势通道”，细砂透镜体则表现出对LNAPL垂向运动的阻碍，蓄积于透镜体上方的LNAPL，在其横向扩展超过透镜体宽度之前以“穿过透镜体”为主。上述不同的迁移过程与介质参数及含水率差异存在较大联系：由于水-气-NAPL三相体系中，NAPL驱替水相、气相发生运移的难易程度与NAPL/水界面的进气压力值密切相关[7,28]，故粗砂介质相对较小的进气压力值Pmax（表1）决定了LNAPL更容易驱替粗砂介质中的水相、气相，因此LNAPL到达粗砂透镜体附近时具有更强的迁移能力。其次，细砂介质较强的持水性显著改变包气带内水分分布，水分占据孔隙空间使细砂透镜体对LNAPL垂向迁移表现出一定阻碍作用；而粗砂透镜体则使-70～-90 cm深度范围内的含水率变化更加显著，促使LNAPL在透镜体及其下方发生显著的横向扩展。

污染物停止释放（180 min）后，粗砂透镜体模型中LNAPL迁移过程与释放阶段一致，而对于细砂透镜体模型，180～200 min内透镜体两侧4个观察点中SNAPL均呈现“迅速上升-逐渐降低”规律（图5），表明LNAPL横向扩展超过细砂透镜体宽度后将转而主要以“绕流”形式运移，释放停止约5.8 h后透镜体上方污染物横向扩展范围逐渐减小，“绕流”减弱。

图5 细砂透镜体周围各观察点SNAPL对比（0～360 min）Fig.5 Comparison of SNAPL at each observation point around the fine-sand lens (0～360 min)

就污染物分布状况而言，180 min后两模型-85 cm深度以下的LNAPL横向扩展情况趋于一致，且两透镜体模型在该深度范围内的含水率分布规律一致（图2），故可推断包气带内的含水率分布规律是控制LNAPL分布形态的主要因素。透镜体岩性差异亦将影响两模型中的LNAPL分布：由于细砂介质孔隙度较大且渗透性较差，10～24 h部分LNAPL在透镜体内部产生滞留现象；而粗砂介质较好的渗透性则决定了其对污染物的滞留能力较差，大部分LNAPL集中分布于透镜体下方。模拟至24 h时，粗砂透镜体模型中污染物分布深度范围相比细砂透镜体模型小28.5%，包气带内LNAPL分布面积（约1 950 cm2）相比细砂透镜体模型（约2 300 cm2）小15.2%。

2.2 水位波动条件下自由相迁移分布

各种原因造成的地下水补排不平衡将导致水位短时波动（如地表水-地下水相互转化、降水入渗等[29]），使包气带内LNAPL迁移分布状况复杂化[15]。本节以细砂、粗砂透镜体模型作为水位波动算例基础，保持LNAPL瞬时释放（4.483 mL/min，释放180 min），通过下层18个网格作为控制水位波动的源汇项，在250 min时控制水位由-100 cm升至-80 cm，并于12 h时恢复至-100 cm。图6(a)、(b)分别代表水位波动条件下细砂、粗砂透镜体模型中LNAPL迁移分布。

图6 地下水位波动情景中不同时刻LNAPL饱和度分布Fig.6 Distribution of LNAPL saturation at different moments in the water table fluctuation scenario

对于细砂透镜体模型，水位上升阶段污染物在浮力作用下向上抬升，至12 h时LNAPL主要分布于透镜体上方并呈连续池状（Pool）[30]，此时透镜体上方（观察点2）SNAPL相比水位恒定时上升0.15。同时，水位升高使LNAPL迁移过程发生显著变化，图7所示为100 min至24 h内透镜体上方（观察点2）、下方（观察点5）以及透镜体内部SNAPL变化状况，对比可知，水位升高时细砂透镜体内部SNAPL相比水位恒定情景减小0.13～0.20，透镜体两侧SNAPL平均值相比水位恒定时增大257%～386%，表明水位升高时LNAPL将主要以“绕流”形式在包气带内运动，其原因为水位升高时细砂透镜体含水率增大至35.8%，更加显著的透镜体高含水率特征导致其对LNAPL迁移的阻碍作用增强。

图7 水位恒定与波动条件下细砂透镜体附近SNAPL对比Fig.7 Comparison of SNAPL near the fine-sand lens under the stationary and fluctuating water table conditions

水位降低时（12～24 h）包气带内含水率分布恢复为水位波动前状况，已大范围横向扩展的LNAPL在重力作用下向下迁移，故18 h后透镜体内部及其周围SNAPL变为与水位恒定情景一致（图7）。经水位波动，-60～-80 cm深度内LNAPL分布范围大幅增加，导致包气带内LNAPL分布面积相比水位恒定情景中增大约51%。此外，24 h时部分LNAPL蓄积于地下水位以上呈透镜体状，其横向分布范围、SNAPL最大值相比水位恒定时分别增大12.5%、5.6%，表明水位开始降低时部分LNAPL仍以“绕流”形式运动，由此造成更多污染物向下部区域迁移。

对于粗砂透镜体模型，水位升高时大部分LNAPL分布于透镜体内部及上方，结合图8可知，250 min后粗砂透镜体内部及其上方SNAPL均出现上升，透镜体内部SNAPL相比水位恒定时平均增大0.19，其原因为：水位升高使透镜体进一步成为LNAPL迁移分布过程的相对“优势空间”。此外，水位升高时透镜体两侧SNAPL相比水位恒定时平均增大0.13～0.16，表明水位抬升使LNAPL转为在透镜体区域发生横向迁移，透镜体所处的-70～-75 cm深度内污染物横向扩展范围相对增大约64.3 cm。

图8 水位恒定与波动条件下粗砂透镜体附近SNAPL对比Fig.8 Comparison of SNAPL near the coarse-sand lensunder the stationary and fluctuating water table conditions

水位降低过程中，粗砂透镜体内部及其上方聚集的LNAPL向透镜体下方迁移，至24 h时同样在地下水位附近蓄积形成LNAPL污染晕。经历水位波动后，该模型包气带内污染物分布面积相较水位恒定时增大约63%，对比细砂透镜体模型可知，水位波动对粗砂透镜体模型中污染物分布的影响更加显著。

2.3 溶解相分布特征

以LNAPL在水相中的质量分数x(aq)作为表征LNAPL溶解相分布特征的定量指标，图9（a）（b）分别代表水位恒定时，两种透镜体模型中溶解相LNAPL的分布变化情况。LNAPL释放阶段两模型中溶解相分布差异主要体现在透镜体附近，细砂透镜体模型观察点1、3中x(aq)平均值均为6.98×10-5，相比粗砂透镜体模型高81%。其次，由于粗砂透镜体的“优势通道”作用以及粗砂较好的渗透性，至180 min时粗砂透镜体下方的溶解相分布面积相比细砂透镜体大46%左右。释放结束后，随自由相的迁移，两模型内溶解相均出现较大扩展，至24 h时细砂透镜体模型包气带内溶解相分布面积（约2 600 cm2）相比粗砂透镜体模型（约2 150 cm2）大20.9%。上述结果表明，两模型包气带内溶解相的分布规律均与自由相一致，但其分布范围的差异相比自由相表现得更为显著。

图9 细砂、粗砂透镜体中不同时刻LNAPL在水相中的质量分数x(aq)Fig.9 Distribution of x(aq) at different moments in the fine sand-lens model

地下水位波动条件下，溶解相分布特征亦与自由相LNAPL密切相关，水位升高时两模型中自由相均集中分布于-70～-85 cm，由此使该深度成为水位波动条件下溶解相分布最广泛的区域。水位波动使溶解相分布状况趋于一致：12 h时两模型在-60 cm以下的溶解相分布面积仅相差约100 cm2，至24 h时两模型的溶解相分布差异仅体现在包气带上部。

2.4 LNAPL挥发过程

LNAPL通过质量交换进入气相，将以挥发方式向上穿过包气带到达地表[31]，本节通过模型顶部单元格计算得到的平均气相VOC通量数据表征LNAPL挥发，图10及图11所示分别为水位恒定、水位波动条件下两模型中挥发通量的变化情况。

图10 不同透镜体模型中LNAPL挥发通量对比Fig.10 Comparison of LNAPL volatile flux in each lens structure model

图11 不同水位情景中LNAPL挥发通量对比Fig.11 Comparison of LNAPL volatile flux in different water table scenarios

水位恒定时，由于污染物横向扩展阻碍了气体-LNAPL质量交换，因此两模型中挥发通量在此阶段均呈降低趋势（图10）。此阶段两模型挥发过程存在一定区别：结合图3、图4对比挥发通量可知，由于粗砂透镜体对LNAPL垂向运动呈“优势通道”作用，透镜体下方污染物分布范围更大，故粗砂透镜体模型中溶解过程更显著，其挥发通量平均值（1.40×10-6mol/(s·m2)）相比细砂模型（1.53×10-6mol/(s·m2)）小8.5%。LNAPL释放停止后，两模型中挥发过程均呈现出“先减弱-后增强”趋势，并在23 h左右达到峰值（1.3×10-6mol/(s·m2)）。结合气体质点方向分析，释放停止后包气带上部的LNAPL对水、气驱替减缓，部分气体快速向原LNAPL下渗路径上运动，填补多孔介质中LNAPL减少所产生的“空缺”，且由于180 min后包气带下部的LNAPL分布范围扩大，故挥发通量上升。

地下水位波动时，由于水位抬升打破饱和带内原有的三相平衡状态，且水位升高时LNAPL将优先驱替多孔介质内的气体[11]，故水位开始抬升时细砂、粗砂透镜体模型中挥发通量分别迅速上升至9.64×10-6mol/(s·m2)、1.09×10-5mol/(s·m2)，此时分别于细砂透镜体两侧、粗砂透镜体上部观察到挥发通量快速增大，部分LNAPL迅速以气相形态向上逸散。挥发通量快速到达峰值后，由于两模型中LNAPL均在透镜体区域横向扩展，阻碍了气体-LNAPL质量交换过程，故此后模型中挥发通量呈逐渐降低趋势。两透镜体模型在水位上升阶段（250 min至12 h）的平均挥发通量相比水位恒定时分别增大124%、126%，表明水位上升将大幅增强LNAPL的挥发过程。水位降低时两模型的挥发通量均出现小幅度波动，结合气体质点运动方向分析，随着水位降低，模型下部介质孔隙内LNAPL及水分均向下迁移，模型上部挥发相随之转为向下部运动，故模型上部挥发通量均有所降低。

3 结论

（1）地下水位恒定时，介质含水率差异导致粗砂透镜体表现出对LNAPL垂向迁移的“优势通道”作用，细砂透镜体则起不完全阻碍作用，使LNAPL呈现“蓄积穿透-横向扩展-绕流”迁移过程。细砂、粗砂透镜体模型中LNAPL主要分布于透镜体内部、下方，两模型包气带内污染物分布面积相差约15.2%。

（2）地下水位波动将使粗砂透镜体进一步成为LNAPL迁移分布的“优势空间”，而其导致的细砂透镜体含水率变化则使“绕流”显著增强。水位波动使-60～-80 cm深度内LNAPL分布范围大幅增加，细砂、粗砂模型包气带内LNAPL分布面积相比水位恒定情景分别增大51%、63%。

（3）溶解相分布规律与自由相迁移密切相关，水位恒定时两模型中溶解相分布差异主要体现在透镜体附近，至24 h时溶解相分布面积相差20.9%，地下水位波动则使该差距显著减小至约100 cm2。

（4）气体-LNAPL接触条件、LNAPL分布状况决定LNAPL挥发强度变化，两模型中挥发强度变化均呈“先减弱-后增强”趋势，前180 min粗砂透镜体模型平均挥发量相比细砂透镜体模型小8.5%。水位波动使包气带内三相平衡状态受扰动，主要表现在水位抬升阶段LNAPL挥发增强，该阶段两模型中平均挥发量相比水位恒定时分别增大124%、126%。