基于T2VOC對(duì)比分析LNAPL多種因素運(yùn)移的數(shù)值模擬

侯岳嵐，夏 源，程亞平,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ù)十年之久，對(duì)地下水環(huán)境系統(tǒng)造成嚴(yán)重破壞。其泄漏后，屬于“NAPL-液-氣”多相流的問(wèn)題，其中：LNAPL自由態(tài)和殘余固態(tài)屬于NAPL 相，溶于水后溶解態(tài)屬于液相，揮發(fā)在空氣中屬于氣相[1-2]。目前大多數(shù)學(xué)者通過(guò)“NAPL-液-氣”三相流體系描述LNAPL在地下環(huán)境遷移形式[3]。

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

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

1 計(jì)算方法

1.1 控制方程

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

式中：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 概念模型

考慮污染物苯泄漏速率、地下水位波動(dòng)、地下水流速3 種影響因素，采用單因素敏感度分析法對(duì)3 種影響因素影響下LNAPL 的溶解速率做敏感度分析，得出各因素對(duì)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）地下水位波動(dòng)。建立3 種水位波動(dòng)模型，分別為情景C-1（無(wú)波動(dòng)）、情景C-2（低波動(dòng)）、情景C-3（高波動(dòng)）。

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

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

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

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

2 結(jié)果與分析

2.1 地下水流速的影響

首先考慮地下水水流速率對(duì)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)分布見(jiàn)圖2。

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

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

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

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

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

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

泄漏強(qiáng)度對(duì)LNAPL 的污染范圍和濃度分布有明顯的影響，故考慮其對(duì)泄漏在地下水環(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 相溶解，故對(duì)NAPL 相苯進(jìn)行監(jiān)測(cè)，結(jié)果見(jiàn)圖5。由圖5 可知，3 種情景中高NAPL相飽和度分別可達(dá)0.342，0.241，0.236。

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

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

圖6 3種情景下各時(shí)刻溶解速率圖

2.3 地下水位波動(dòng)的影響

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

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

圖7 3種情景下各時(shí)刻溶解速率圖

2.4 敏感度分析

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

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

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

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

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

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

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

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

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

綜合圖8～圖10 可知，敏感度系數(shù)E1、E2＞E6、E5＞E7、E8＞E3、E4，故對(duì)NAPL相質(zhì)量溶出速率影響的因素對(duì)比如下：每增大一倍地下水流速率＞地下水位每上升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é)語(yǔ)

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

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

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