侯岳嵐，夏 源，程亞平,3，周正濤

（1.桂林水文中心，廣西 桂林 541199；2.桂林理工大學(xué)，廣西 桂林 541006；3.巖溶地區(qū)水污染控制與用水安全保障協(xié)同創(chuàng)新中心，廣西 桂林 541006)

0 引言

不與包氣帶及潛水面以下的水發(fā)生混合的液體通常稱為非水相液體（nonaqueous phase liquids，NAPLs），其中密度比水小的叫做輕非水相液體（Light Non Aqueous Phase Liquid，LNAPL）。LNAPL屬于揮發(fā)性有機(jī)污染物（VOCs），不與水混合且很難自然降解，在地下水系統(tǒng)中能停留數(shù)年至數(shù)十年之久，對地下水環(huán)境系統(tǒng)造成嚴(yán)重破壞。其泄漏后，屬于“NAPL-液-氣”多相流的問題，其中：LNAPL自由態(tài)和殘余固態(tài)屬于NAPL 相，溶于水后溶解態(tài)屬于液相，揮發(fā)在空氣中屬于氣相[1-2]。目前大多數(shù)學(xué)者通過“NAPL-液-氣”三相流體系描述LNAPL在地下環(huán)境遷移形式[3]。

LNAPL 污染物遷移過程復(fù)雜，受多種因素影響。施小清、李曄等[4-5]運(yùn)用T2VOC程序、NAPL sim?ulator 模型與土壤相對介電常數(shù)等方法研究出NAPLs 泄漏速率對其在非飽和帶中的運(yùn)移范圍影響。楊明星[6]為研究水位波動帶對不同類型的石油污染物產(chǎn)生成分和分布的變化，利用TMVOC 程序數(shù)值模擬和室內(nèi)試驗(yàn)進(jìn)行驗(yàn)證，不同污染物組分在非均質(zhì)土壤中滲入滯留在界面附近，污染物自由態(tài)的濃度也受水位波動下水流沖刷作用逐漸升高。潘明浩等[7]為揭示透鏡體與地下水波動相互作用下非水相液體的遷移過程，利用TOUGH2 程序中TMOVC 模塊建立LNAPLs 在包氣帶運(yùn)移的多相流數(shù)值模型，結(jié)果表明：地下水位恒定時，粗砂透鏡體有助于LNAPL 垂向遷移。王穎等[8]基于TMVOC 模型模擬分析LNAPL 受水位波動所引起的污染范圍與質(zhì)量變化。

就目前的研究可知，多種因素影響著LNAPL的溶解速率，但對其受因素的變化及影響的敏感度研究還是空白。本文采用TOUGH2 軟件對LNAPL 進(jìn)行數(shù)值模擬，選取地下水流速率、LNAPL 的泄漏強(qiáng)度、地下水位波動等影響因素，對比分析上述幾種因素對LNAPL 溶解速率的變化及其敏感度，為LNAPL的溶解遷移模擬提供源強(qiáng)依據(jù)。

1 計(jì)算方法

1.1 控制方程

由于多相流模型具有高度非線性的特點(diǎn)，本文采用由美國勞倫斯—伯克利國家實(shí)驗(yàn)室（LBNL）開發(fā)的TOUGH2 軟件中T2VOC 模塊求解，T2VOC 模塊求解采用數(shù)值法，適用于多相流模型[9]。在TOUGH（非飽和地下水流及熱流傳輸）數(shù)值模擬程序中，可以分別模擬一維，二維和三維孔隙或裂隙介質(zhì)中多相流、多組分及非等溫的水流及熱量運(yùn)移。T2VOC 程序運(yùn)用積分有限差分法（IFDM）進(jìn)行空間離散，將水、氣及NAPL視為三相且假設(shè)均符合Darcy定律，三相運(yùn)移根據(jù)壓力和重力確定[10]。在包含三個質(zhì)量分量的非等溫系統(tǒng)中，需要三個質(zhì)量平衡方程和一個能量平衡方程來完整描述系統(tǒng)，對于表面積為n的任意流動區(qū)域Vn，以積分的形式寫出：

式中：Mk為每單位多孔介質(zhì)體積中組分K（K=w，a，c）的質(zhì)量；Fk為組分K進(jìn)入體積Vn的質(zhì)量通量；n為內(nèi)法線方向的單位向量；qk為單位體積內(nèi)組分K的源匯項(xiàng)。

1.2 概念模型

考慮污染物苯泄漏速率、地下水位波動、地下水流速3 種影響因素，采用單因素敏感度分析法對3 種影響因素影響下LNAPL 的溶解速率做敏感度分析，得出各因素對LNAPL溶出質(zhì)量的敏感度。建立以下情景：

（1）地下水流速。建立3種情景地下水流速率模型，分別為情景A-1（低流速）、情景A-2（中流速）、情景A-3（高流速）。

（2）泄漏強(qiáng)度。建立3 種情景泄漏強(qiáng)度模型，分別為情景B-1（弱泄漏）、情景B-2（中泄漏）、情景B-3（強(qiáng)泄漏）。

（3）地下水位波動。建立3 種水位波動模型，分別為情景C-1（無波動）、情景C-2（低波動）、情景C-3（高波動）。

由于模擬的情景多，計(jì)算量大，而敏感度分析的目的是為了評價影響因素對溶出質(zhì)量的影響強(qiáng)弱，對模擬的精度要求相對不高，因此對模型縮小計(jì)算尺度，并減少剖分的網(wǎng)格數(shù)量，以減少計(jì)算量。假設(shè)模擬區(qū)域?yàn)榫|(zhì)中砂，尺寸為10 m×3 m×3 m的矩形砂槽（見圖1），X、Y、Z方向均勻剖分為25×15×15共5625個單元格，其中單位單元格尺寸為0.4 m×0.2 m×0.2 m，頂部邊界為封閉邊界，左右邊界為定水頭邊界，下部邊界為隔水邊界。選取苯（benzene）為本次研究的LNAPL污染物。

為了便于對比分析，不考慮熱物理相關(guān)參數(shù)，苯的理化性質(zhì)參數(shù)選取TOUGH2 程序提供的VOC數(shù)據(jù)集里面的部分參數(shù)（見表1）。

表1 模型中主要參數(shù)

初始模型污染源尺寸為0.4 m×0.2 m 的矩形設(shè)于模型區(qū)域頂部中心軸下游3 m 處，污染物苯以5.55×10-4kg/s速度持續(xù)泄漏6 h，改變初始條件考慮地下水流速、滲透速率、地下水位波動3 種因素，比較他們對LNAPL的溶解速率的影響。

2 結(jié)果與分析

2.1 地下水流速的影響

首先考慮地下水水流速率對LNAPL 運(yùn)移的影響，利用改變左右邊界的水頭差調(diào)整水力梯度，采用達(dá)西定律計(jì)算出地下水流速。情景A-1水流速率為0.38 m/d；情景A-2 水流速率為0.77 m/d；情景A-3水流速率為1.54 m/d。3種情景下壓強(qiáng)分布見圖2。

圖2 3種情景下壓強(qiáng)分布

模型運(yùn)行6 h時，污染源苯停止泄漏，此時液相飽和度分布見圖3。由圖3 可知，3 種情景下XZ剖面液相飽和度的面積分別為5.252、5.566、6.0 615 m2。其中高飽和液相飽和度面積分別為0.1 708、0.2 149、0.2 520 m2。對比情景A-1的液相飽和度面積，情景A-2 中2 倍水流速率增加5.98%，情景A-3中4倍水流速率增加15.41%，相對于情景A-2約2.5倍。對比情景A-1中高飽和液相飽和度面積，情景A-2增加25.82%，情景A-3增加47.54%，相對于情景A-2約2倍左右。結(jié)合質(zhì)量溶出速率與NAPL-水界面面積之間的關(guān)系，可知呈倍數(shù)的增加地下水水流速率促使LNAPL 中NAPL 相溶解到水中數(shù)量呈倍數(shù)的增加，其質(zhì)心面積變異程度受地下水流速率影響較大。

圖3 3種情景下苯泄漏6 h液相飽和度分布

持續(xù)監(jiān)測受此影響因素下LNAPL 溶解速率變化（見圖4）。由圖4 可知，模型運(yùn)行前13 h 溶解速率數(shù)值波動較大，情景A-3 基本處于3 種情景中的較大值，13 h后LNAPL的NAPL相達(dá)到穩(wěn)定狀態(tài)，溶解速率趨于穩(wěn)定，可以看出：地下水流速率與溶解速率呈現(xiàn)正相關(guān)趨勢。

圖4 3種情景下各時刻溶解速率變化

2.2 泄漏強(qiáng)度的影響

泄漏強(qiáng)度對LNAPL 的污染范圍和濃度分布有明顯的影響，故考慮其對泄漏在地下水環(huán)境中LNAPL 的溶解速率的變化，保持泄漏總量不變，改變泄漏速率設(shè)置3 種情景。情景B-1 模型與情景A-1模型一致，苯泄漏速率為5.55×10-4kg/s，持續(xù)泄漏6 h；情景B-2 中苯泄漏速率增加5 倍，為2.775×10-3kg/s，持續(xù)泄漏72 min；情景3 中苯泄漏速率增加10 倍，為5.55×10-3kg/s，持續(xù)泄漏36 min。模型運(yùn)行6 h 后，3 種情景下污染物苯都已泄漏完成，因LNAPL 的液相主要由NAPL 相溶解，故對NAPL 相苯進(jìn)行監(jiān)測，結(jié)果見圖5。由圖5 可知，3 種情景中高NAPL相飽和度分別可達(dá)0.342，0.241，0.236。

圖5 3種情景下模型運(yùn)行6 h后NAPL相飽和度分布

3 種情景下各時刻溶解速率圖見圖6。由圖6可知，前15 h曲線波動劇烈，情景1中溶解速率曲線值忽高忽低，這是由于其泄漏速率低，持續(xù)泄漏量小，受地下水流速率影響造成一定波動；情景B-2中前2 h由于隨著污染物的持續(xù)泄漏導(dǎo)致溶解速率曲線值有一個變大的波動值，待泄漏停止，溶解速率降低，后隨著模型運(yùn)行溶解速率趨向穩(wěn)定。對20 h后曲線值進(jìn)行放大，可以看出模型穩(wěn)定后，泄漏速率越大，溶解速率越低，LNAPL中NAPL相對越少溶解于水中，造成更多殘留相滯留于地下水環(huán)境中。結(jié)合圖5 可知，泄漏速率與NAPL 相飽和度呈現(xiàn)反相關(guān)，NAPL相飽和度與溶解速率呈現(xiàn)正相關(guān)。

圖6 3種情景下各時刻溶解速率圖

2.3 地下水位波動的影響

LNAPL 泄漏在地下水環(huán)境中主要分布在毛細(xì)水帶與地下水位波動帶附近，潛水面波動能夠?qū)ζ湎鄳B(tài)分配轉(zhuǎn)化產(chǎn)生重大影響。鑒于此，模型中通過上邊界注水模擬降雨條件，建立3 種情景（見表2）對地下水位波動因素引起的LNAPL 溶解速率的變化進(jìn)行對比分析。情景C-1 與情景A-1、情景B-1模型初始條件一致；情景C-2 中考慮日降雨量50 mm（大雨），降雨速率為4.6 296×10-5kg/s，待情景C-1中污染物苯泄漏完成，通過降雨1d使得LNAPL在水位波動狀態(tài)下自由擴(kuò)散以致達(dá)到穩(wěn)定狀態(tài)；情景C-3模擬情景與情景C-2相似，降雨1 d考慮日降雨量為100 mm（大暴雨），降雨速率為9.2 592×10-5kg/s。3種情景下各時刻溶解速率圖見圖7。由圖7可知，3種情景中情景C-3降雨量最大，其處于水位上升狀態(tài)時，LNAPL 中NAPL 相的溶解速率處于最高值，而當(dāng)其處于水位下降直至穩(wěn)定過程中，溶解速率處于3種情景中的最小值。水位上升引起的溶解速率增加有兩點(diǎn)原因：第一點(diǎn)是由于降雨引起的淋濾作用下毛細(xì)水帶中污染物L(fēng)NAPL浸涂區(qū)里NAPL相被沖刷出來，導(dǎo)致LNAPL 的NAPL 相飽和度增加，故增大苯自由相的溶解速率；第二點(diǎn)是由于水位上升為自由態(tài)LNAPL提供浮力，驅(qū)替土壤孔隙中自由態(tài)漂浮出來，也對LNAPL 的NAPL 相飽和度增加起關(guān)鍵作用。情景C-3 中水位下降攜帶更多的NAPL相苯使得溶解速率有一定的增幅，待水位穩(wěn)定后，由于先前的一系列波動，NAPL相飽和度已經(jīng)處于3種情景中的最低值，故溶解速率會處于最低值。

表2 不同降雨條件下模擬過程設(shè)計(jì)

圖7 3種情景下各時刻溶解速率圖

2.4 敏感度分析

綜上，地下水流速率、泄漏速率、地下水位波動均對LNAPL 的溶解速率產(chǎn)生一定的影響。單因素敏感性分析可用敏感度系數(shù)表示項(xiàng)目評價指標(biāo)對不確定因素的敏感程度[11-13]，其表達(dá)式為：

式中：ΔF為不確定因素F的變化率；ΔA為不確定因素F發(fā)生ΔF變化率時，評價指標(biāo)A的相應(yīng)變化率；E為評價指標(biāo)A對于不確定因素F的敏感度系數(shù)。

本文E為無量綱數(shù)，如計(jì)算出正負(fù)值，則做絕對值處理（方便觀測比較）。對地下水流速率因素，敏感度系數(shù)E表示NAPL相質(zhì)量溶出速率變化倍數(shù)比地下水流速增大倍數(shù)，其中敏感度系數(shù)E1、E2中地下水流速增大倍數(shù)分別為2 倍、4 倍，對E1、E2進(jìn)行回歸分析，結(jié)果見圖8 和表3。由圖8 和表3 可知，擬合預(yù)測模型1、模型2的R2分別為0.57、0.85，均擬合效果較好。從擬合結(jié)果來看，地下水流速率每增大1 倍，其敏感度系數(shù)E1、E2 變化在0.1～0.3 之間；當(dāng)?shù)叵滤魉俾试龃? 倍時，擬合出敏感度系數(shù)在區(qū)間變化中相對較小。

表3 地下水流速、泄漏強(qiáng)度、地下水位波動的敏感度系數(shù)擬合預(yù)測模型

圖8 地下水流速率敏感度系數(shù)分布

對LNAPL 泄漏強(qiáng)度因素，敏感度系數(shù)E表示NAPL相質(zhì)量溶出速率倍數(shù)比其LNAPL泄漏強(qiáng)度增大倍數(shù)，其中敏感度系數(shù)E3、E4中LNAPL泄漏強(qiáng)度增大倍數(shù)分別為5 倍、10 倍。對E3、E4 所得的敏感度數(shù)據(jù)結(jié)果進(jìn)行觀測，發(fā)現(xiàn)在22～44 h 之間E值波動較為劇烈，44 h 后，E3、E4 值均趨向于0.0 085，故對在22～44 h之間的E值采用3項(xiàng)式線性擬合，結(jié)果見圖9和表3。由圖9和表3可知，回歸模型3、模型4 的預(yù)測值與敏感度實(shí)際值之間的R2分別為0.88、0.93，擬合效果較好。回歸模型3 敏感度系數(shù)在0.004～0.012之間，模型4敏感度系數(shù)在0.002～0.011之間，綜合來看，E3＞E4，說明采用10倍LNAPL泄漏強(qiáng)度計(jì)算敏感度比用其5倍泄漏強(qiáng)度所得的敏感度系數(shù)低，模型運(yùn)行后期，泄漏強(qiáng)度對LNAPL 的敏感度系數(shù)較低，趨向于0.0 085。

圖9 LNAPL滲透速率敏感度系數(shù)分布

對于地下水位波動因素的敏感度分析，由圖7可知，降雨造成的水位波動下使得LNAPL的質(zhì)量溶出速率受水位上升下降有不同的變化，故其敏感度按水位上升和水位下降進(jìn)行分析（見圖10）。假設(shè)水位上升0.05 m為水位上升波動1倍，水位上升0.1 m為波動2倍，圖10（a）中E表示NAPL相質(zhì)量溶出速率倍數(shù)比水位上升下降波動倍數(shù)，其中敏感度系數(shù)E5、E6 中水位上升波動倍數(shù)分別為1 倍、2 倍，敏感度系數(shù)E7、E8中水位下降波動倍數(shù)分別為1倍、2倍。觀察E5、E6敏感度系數(shù)分布，選用一元線性模型進(jìn)行擬合，結(jié)果顯示R2分別為0.62、0.71，擬合結(jié)果較好，其E值為0～0.15，E6 斜率比E5 較大，說明水位上升較高時計(jì)算所得E值較大。圖10（b）中E7、E8選用一元線性模型進(jìn)行擬合，結(jié)果顯示R2分別為0.91、0.93，擬合效果良好，其E值為0.02～0.08，可以看出水位下降對NAPL相質(zhì)量溶出速率影響最高可達(dá)0.08，后隨時間變化逐漸降低。

圖10 地下水位波動敏感度系數(shù)分布

綜合圖8～圖10 可知，敏感度系數(shù)E1、E2＞E6、E5＞E7、E8＞E3、E4，故對NAPL相質(zhì)量溶出速率影響的因素對比如下：每增大一倍地下水流速率＞地下水位每上升0.05 m＞地下水位每下降0.05 m＞LNAPL的泄漏強(qiáng)度每增大一倍。

從表3 中的擬合模型表達(dá)式可以看出，地下水位下降的敏感度最符合線性規(guī)律，即地下水位下降的程度與溶出質(zhì)量近似線性負(fù)相關(guān)關(guān)系；其次是地下水位上升的敏感度，也比較符合線性規(guī)律，即地下水位上升的程度與溶出質(zhì)量近似線性正相關(guān)關(guān)系；而地下水流速和泄露強(qiáng)度的敏感度的非線性規(guī)律較強(qiáng)，都是先上升后下降。

3 結(jié)語

（1）LNAPL 的NAPL 相達(dá)到穩(wěn)定狀態(tài)時，地下水流速率與溶解速率呈現(xiàn)正相關(guān)趨勢；泄漏速率與NAPL相飽和度呈現(xiàn)反相關(guān)，NAPL相飽和度與溶解速率呈現(xiàn)正相關(guān)，NAPL 相飽和度對研究溶解速率變化起關(guān)鍵作用。

（2）影響LNAPL 溶解速率因素比較結(jié)果：地下水流速率每增加1 m/s 的溶解速率變化值＞地下水位每上升1 m 的溶解速率變化值＞LNAPL 的滲透速率每增加1 kg/s 的溶解速率變化值＞地下水位每下降1 m的溶解速率變化值。

（3）對于地下水流速率較快與經(jīng)常造成地下水位上升的污染源區(qū)，質(zhì)量溶出效率較高，造成溶解相擴(kuò)散面積大，污染羽范圍廣。而對于滲透速率較快的污染源區(qū)，LNAPL 的NAPL 相進(jìn)入潛水面附近較多，污染源區(qū)范圍大，后期也會造成溶解相的增多。因此，對于地下水流速率較快與經(jīng)常造成地下水位上升的污染源區(qū)，應(yīng)當(dāng)格外注意擴(kuò)大液相的監(jiān)測，對于滲透速率較快的污染源區(qū)，應(yīng)當(dāng)注意排查LNAPL的NAPL相。