轻非水相液体污染源区结构的影响因素数值分析

陶佳辉，施小清，康学远，徐红霞，吴吉春

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

随着社会对石油化工产品需求量的增加，大量以轻非水相液体(light non-aqueous phase liquid，LNAPL)形式存在的有机产品，在贮存及使用过程中发生泄漏，导致地下水污染[1～2]。一般而言，LNAPL密度比水小且难溶于水，泄漏后在重力与毛细压力作用下运移[3]。若LNAPL泄漏量足够大，少量以残余态形式滞留在运移路径中，部分向下运移直至遇到障碍(黏土层等)，形成轻油聚集体，其余LNAPL运移至毛细水带和地下水位后产生侧向运移[4]。通常将潜水面以上残留LNAPL的区域称为不饱和污染源区；将饱和带或潜水面以下残留LNAPL的区域称为饱和污染源区[5～6]。研究者常采用LNAPL的饱和度分布(saturation distribution)、孔隙中残余LNAPL体积(volume of possible LNAPL residence)、不连续离散状(ganglia)与连续池状(pool)LNAPL体积比(ganglia-to-pool ratio，GTP)等指标来定量表征LNAPL的污染源区结构[7]。残留在介质中的LNAPL在自然条件下难以降解，从而成为长期存在的土壤-地下水污染源。因此，精细刻画污染源区结构对于LNAPL污染的风险评估和修复治理极其重要[8]。

国内外学者主要通过试验与数值模拟方法研究含水量、非均质性、毛细压力等因素对LNAPL运移和分布的影响[9～10]。如Schroth等[11]构建了含有毛细障碍栅的砂箱实验，研究了LNAPL在包气带非均质介质中的运移，结果表明含水量对LNAPL的运移和分布有较大影响。Wipfler等[12]通过二维砂箱实验观察LNAPL在层状多孔介质中的运移，结果表明毛细压力是影响LNAPL运移和分布的主要敏感因子。束善治等[4]采用数值分析方法探究了不同特征的土层结构对LNAPL运移和分布的影响。Van Geel等[13]通过两相和三相试验数据来预测LNAPL的残余饱和度，结果表明LNAPL的残余饱和度是最大流体饱和度和含水量的函数。

以往研究中通常视LNAPL从地表泄漏后穿过包气带运移至潜水面，前人的分析和探讨大多建立在“LNAPL在潜水面漂浮移动”的污染概念模型上。然而，地下介质的非均质性和包气带含水量的空间变异性可能形成复杂的LNAPL污染源区结构，LNAPL甚至可能无法到达潜水面，而在毛细水带或局部透镜体蓄积。因此，探讨多种因素对复杂LNAPL污染源区结构的影响具有重要意义。

本文基于数值模型研究LNAPL泄漏量、介质非均质性与含水量空间变异、潜水面周期性变化等因素对LNAPL污染源区结构的影响。通过分析不同条件下LNAPL的饱和度分布来刻画LNAPL污染源区结构，以期为地下水中LNAPL污染的风险评估及修复治理提供理论参考。

1 计算方法

1.1 控制方程

本文采用TOUGH2软件的T2VOC模块模拟LNAPL的运移与分布[14]。TOUGH2可以解决不同条件下地下水中水流和热量运移问题，而T2VOC模块在模拟饱和带、非饱和带非水相液体运移的问题中应用广泛。

T2VOC程序将NAPL视为一个单独的相，采用有限差分法进行空间离散，水、气及NAPL中每一相的运移可根据Darcy定律的多相形式由压力和重力确定[15]。在三相组分组成的非等温系统中，完整地描述该系统需要三个质量守恒方程和一个能量守恒方程。对于任意体积为Vn，表面积为Γn的渗流区，组分k的平衡方程为：

(1)

式中：Mk——单位体积组分k(k=w,a,n)的质能累积量；

Fk——组分k进入体积Vn的质能通量；

n——内法线方向的单位向量；

qk——单位体积内组分k的源汇项。

1.2 概念模型

假设模拟区域为二维XZ剖面，尺寸为90 cm×5 cm×90 cm，X方向Z方向均匀剖分为30×30共900个网格，网格水平与垂直间距均为3 cm。模型上边界为大气边界，左右边界为定水头边界，下边界为隔水边界(图1)。甲苯密度比水小，是常见LNAPL污染物汽油中的典型组分，对人体与自然环境有较大危害，故本研究选取甲苯作为LNAPL的特征污染物[16]。污染源设于模拟区域顶端以下30 cm的中心处，泄漏速率恒为5×10-4kg/s，基于甲苯的运移和分布特征分析LNAPL污染源区结构的影响因素。假设温度恒为20 ℃，故未考虑热物理相关参数。

相对渗透率函数采用Stone模型，Stone模型中Swr、Snr、Sgr分别表示液相、NAPL相以及气相的残余饱和度，n为拟合参数。毛管压力函数采用Parker模型，Parker模型中Sm表示残余液相饱和度,αgn和αnw分别表示气相和NAPL相、NAPL相和水相之间进气压力的倒数。

图1 概念模型(WT：潜水面)Fig.1 Conceptual model(WT: Water Table)

粗砂中砂细砂黏土透镜体孔隙率0.340.380.410.5渗透率/10-10 m22.260.1930.080.001颗粒比重/(kg·m-3)1 7491 6501 5642 600Swr0.050.100.120.60Stone相对Snr0.030.030.040.05渗透率模型Sgr0000.001n3333Sm0.050.060.100.60Paker毛管压力模型n1.841.841.841.84αgn10.810.011.0500.0αnw11.005.001.580.05

2 影响因素分析

2.1 LNAPL泄漏量的影响

LNAPL的运移和分布受泄漏量以及泄漏条件(事故性泄漏、多年持续泄漏等)的影响。随着泄漏量的不断增加，部分LNAPL蓄积在毛细水带中，其余LNAPL将运移至潜水位以下[19]。基于不同的LNAPL泄漏量，本节分析了泄漏量对LNAPL源区结构的影响。模型介质为均质中砂，系统达到水-气平衡后，以5×10-4kg/s的速度注入甲苯，注入时间分别为500 s、1 000 s和2 000 s，分析停止注入后静置条件下LNAPL的运移和分布特征。模拟进行240 min后LNAPL的空间分布趋于稳定，图2为不同泄漏量下LNAPL饱和度分布图。可见泄漏量较小时，LNAPL在毛细水带中积蓄的压力较小，几乎所有的LNAPL都蓄积在毛细水带中。随着泄漏量的增大，LNAPL将穿越毛细水带，少量LNAPL运移至潜水面以下。

2.2 非均质性的影响

本节设计两种不同的层状非均质模型，并与均质中砂模型(图3a)作对比。两种层状非均质模型由两层砂土构成：上细下粗(图3b)、上粗下细(图3c)。假设甲苯的泄漏量较小，注入时间为500 s，三种模型的边界条件、甲苯的注入速率和注入方式均相同。

系统在稳定水流场中达到水-气平衡后(不考虑降雨入渗)，介质非均质性会导致潜水面以上毛细水带厚度的差异(图4)。潜水位标高15 cm；潜水面以上，含水量接近饱和的区域为毛细水带；毛细水带上边缘与大气边界之间的区域为土壤水带。三种模型毛细水带厚度分别为10，7，13 cm。上粗下细模型中毛细水带厚度最大(图4c)，原因在于：细砂层土颗粒粒径较小，比表面积相对较大，产生较大的毛细压力，使得细砂层对水有较强的吸附和滞留能力，通过自吸方式获得较多水分。

图2 不同泄漏量下LNAPL饱和度(SO)分布(CF：毛细水带边缘)Fig.2 LNAPL saturation(SO)distribution with different leakage (CF: Capillary Fringe)注：为显示LNAPL主要分布区域，将SO=0.01作为饱和度阈值，下同。

图3 均质和层状非均质概念模型Fig.3 Conceptual model showing the homogeneous and layered heterogeneity注：实线代表细砂/粗砂、粗砂/细砂分界线，下同。

图4 不同介质条件下水的饱和度(Sw)分布Fig.4 Water saturation(Sw)distribution under different medium situations注：虚线代表地下水位线和毛细水带上边缘，下同。

图5为三种不同介质中不同时刻LNAPL饱和度分布图。均质介质中(图5a)，渗流初始阶段(T=10 min)，LNAPL受重力作用影响较大，以垂向入渗为主。T=60 min时，LNAPL受毛细压力影响较大，开始以横向迁移为主，并在毛细水带边缘达到最大饱和度0.38。模拟结束时大部分LNAPL进入到毛细水带中，但由于LNAPL密度比水小且泄漏量较小，最终滞留在毛细水带中。

上细下粗模型中(图5b)，T=30 min时，LNAPL开始在细砂层下部蓄积，T=60 min时，细砂层中蓄积的LNAPL达到最大饱和度0.38。模拟结束时大部分LNAPL穿越细砂层进入到毛细水带中，但仍有部分LNAPL滞留在细砂层无法下渗。

上粗下细模型中(图5c)，渗流初始阶段由于粗砂层较大的渗透率和颗粒粒径，LNAPL在粗砂较快地运移至毛细水带边缘，随后在毛细水带边缘蓄积并横向迁移。模拟结束时LNAPL均匀水平地展布在毛细水带上边缘，几乎没有LNAPL进入到毛细水带中。

图5 三种介质中不同时刻LNAPL饱和度(SO)分布Fig.5 LNAPL saturation(SO)distribution at different times in three kinds of media

图6为模拟结束时三种介质中水和LNAPL的饱和度分布随深度变化的剖面图。上细下粗模型中(图6b)细砂层阻碍LNAPL的垂向入渗，最终仍有部分LNAPL(最大饱和度达0.22)滞留于细砂层无法向下运移，成为长期污染源。原因是细砂层固有渗透率低，LNAPL在细砂中运移较为困难；细砂层部分孔隙体积被水占有，含水量较大，LNAPL的相对渗透率较低[20]；细砂层土颗粒粒径较小，较大的毛细压力驱使LNAPL产生较多的横向迁移，阻碍LNAPL的垂向入渗。

上粗下细模型中(图6c)，LNAPL最终蓄积在细砂层毛细水带边缘。LNAPL无法进入到细砂层中的毛细饱和带，原因是细砂层土颗粒粒径较小，LNAPL相对渗透率较低，不容易发生运移；细砂层毛细水带含水量较高，对LNAPL产生向上的浮力作用，阻碍LNAPL的垂向入渗；细砂层的进气压力较大，LNAPL与水之间的压力差不足以克服细砂颗粒对水的毛细束缚，最终在毛细水带边缘蓄积[21]。

一般认为LNAPL在包气带中可直接运移至潜水面，然而图6(b)、6(c)显示，上细下粗非均质条件下，部分LNAPL滞留在细砂层中无法下渗，形成双重LNAPL高饱和度区；上粗下细条件下几乎所有的LNAPL均滞留在毛细水带边缘，不会运移至潜水面。两者产生差异的原因在于：包气带中LNAPL的运移和分布是油-水-气三相相互驱替的过程，水相对于空气具有更大的密度和黏度，LNAPL到达毛细水带后更容易驱替空气向四周迁移，并在毛细水带处发生蓄积[22]。由于LNAPL相对空气而言是浸润相，而相对水而言是非浸润相，且细砂层固有渗透率较低，进气压力较大，阻碍LNAPL的下渗，因此最终其无法穿越细砂层的毛细水带。介质非均质性将使LNAPL污染源区结构趋于复杂，增大污染场地修复治理的难度[23]。

图6 水的饱和度(Sw)和LNAPL饱和度(SO)分布对比Fig.6 Comparison of water saturation(Sw)and LNAPL saturation(SO)distribution

2.3 含水量与非均质性的共同影响

图7 含透镜体概念模型及水的饱和度(Sw)分布Fig.7 Conceptual model showing the clay lens and water saturation(Sw)distribution

本节在均质中砂背景中设置黏土透镜体(图7a)，分别考虑饱和与非饱和两种不同含水量情形。为避免毛细水带对结果的影响，将模拟范围扩大为原来的100倍，同时将地下水位降为10 m，这样可使得LNAPL运移区域远离毛细水带。饱和透镜体模型考虑降雨，系统首先在持续降水入渗条件下达到平衡，然后停止降雨入渗，以5×10-4kg/s的速度注入甲苯，分析持续注入1 a过程中甲苯的运移和分布情况。系统在降雨条件下达到水-气平衡后，背景中砂的含水饱和度(Sw=15%～19%)仅略高于残余饱和度而黏土透镜体基本达到饱和(图7b，Sw>98%)[24]。非饱和透镜体模型不考虑降雨条件，黏土透镜体中含水量为其残余饱和度(图7c)。

图8 含透镜体模型不同时刻LNAPL饱和度(SO)分布Fig.8 LNAPL saturation(SO)distribution at different times in the clay lens model

图8为含透镜体模型不同时刻LNAPL饱和度分布图。在饱和与非饱和两种不同含水量情形下，LNAPL均在黏土透镜体上方蓄积并侧向运移形成污染池，使得水平扩散范围显著增大。两种情形仅含水量存在差异，但二者LNAPL污染源区结构存在较大差异。污染土体中一般同时存在气相、水相和NAPL相，当LNAPL蓄积的压力超过LNAPL/水界面的“进气压力”时，LNAPL才能驱替孔隙中的水相或气相[4,25～26]。由于饱和透镜体含水量较高，毛细压力较大，从模拟结果来看，LNAPL无法驱替饱和透镜体中的水分和空气，饱和黏土透镜体对于LNAPL来说是不可渗的(图8a)。而在非饱和透镜体情形下，LNAPL可以驱替非饱和透镜体中的水分和空气，并穿透非饱和透镜体继续向下入渗(图8b)。由此可见包气带中黏土透镜体并非都是LNAPL运移的阻碍，当黏土透镜体中含水量低时，LNAPL可以穿透黏土透镜体；而黏土透镜体含水量较高时，将阻碍LNAPL的入渗。

2.4 潜水面周期性变化的影响

地下介质中的LNAPL主要分为自由相、残留相与溶解相三种[27～28]。由于LNAPL多集中分布于毛细水带及地下水位波动带，不同相态之间的分配转化主要是由潜水位波动造成的[29～30]。鉴于此，本节在均质中砂模型的基础上，将最下层30个网格作为注水和抽水的水位波动控制单元格，通过改变其源汇项条件来升高/降低地下水位(h)，水位先上升至0.3 m，再下降至0 m，分析地下水位一个升降周期内LNAPL源区结构特征的改变。

图9为地下水位波动时LNAPL饱和度分布图，可见水位波动使LNAPL污染范围变大。波动初期水位上升至0.3 m(图9b)，由于LNAPL密度比水小，水位上升过程中油-水-气三相相互驱替，水位上升产生的浮力使LNAPL一同抬升。与水位波动前相比，被携带上升的LNAPL的横向分布范围明显扩大。水位下降至0 m时(图9c)，随着水位下降LNAPL受重力作用垂向入渗。水位下降过程中，重力使LNAPL向下运移使得垂向分布范围变大，又由于LNAPL向下驱替水比较困难，而向四周驱替空气相对容易，故LNAPL横向分布范围变大，最终LNAPL的分布范围逐渐扩大至整个水位波动带[31]。

地下水位波动导致土壤含水量以及其他环境因素的改变，水位波动带内油-水-气三相相互作用促进了地下NAPL相的重新分配[32]。水位上升会为LNAPL相提供浮力，驱替土颗粒内的LNAPL相从孔隙中漂浮出来；水位下降时LNAPL受重力作用向下层运移，使得污染范围明显增大。附着于土壤的LNAPL随地下水位波动而不断地被溶解和裹挟，使污染源区结构趋于复杂。

图9 地下水位波动时LNAPL饱和度(SO)分布Fig.9 LNAPL saturation(SO)distribution under the water table fluctuations

3 结论

(1)当泄漏量较大时，LNAPL可运移至潜水面。

(2)当泄漏量较小时，均质条件下，LNAPL主要在潜水面上方的毛细水带滞留；上细下粗条件下，LNAPL在土壤水带及毛细水带均会滞留；上粗下细条件下，LNAPL主要在毛细水带上边缘发生蓄积，无法到达潜水面。

(3)包气带中黏土透镜体并非都是LNAPL运移的阻碍，LNAPL可以穿透低含水量的黏土透镜体，只有高含水量的黏土透镜体才对LNAPL的入渗有阻碍作用。

(4)潜水面周期性变化导致地下LNAPL分布范围变大，扩大污染范围。