基于T2VOC對比分析LNAPL多種因素運移的數值模擬

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

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

0 引言

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

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

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

1 計算方法

1.1 控制方程

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

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

1.2 概念模型

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

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

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

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

由于模擬的情景多，計算量大，而敏感度分析的目的是為了評價影響因素對溶出質量的影響強弱，對模擬的精度要求相對不高，因此對模型縮小計算尺度，并減少剖分的網格數量，以減少計算量。假設模擬區域為均質中砂，尺寸為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污染物。

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

表1 模型中主要參數

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

2 結果與分析

2.1 地下水流速的影響

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

圖2 3種情景下壓強分布

模型運行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倍左右。結合質量溶出速率與NAPL-水界面面積之間的關系，可知呈倍數的增加地下水水流速率促使LNAPL 中NAPL 相溶解到水中數量呈倍數的增加，其質心面積變異程度受地下水流速率影響較大。

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

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

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

2.2 泄漏強度的影響

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

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

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

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

2.3 地下水位波動的影響

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

表2 不同降雨條件下模擬過程設計

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

2.4 敏感度分析

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

式中：ΔF為不確定因素F的變化率；ΔA為不確定因素F發生ΔF變化率時，評價指標A的相應變化率；E為評價指標A對于不確定因素F的敏感度系數。

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

表3 地下水流速、泄漏強度、地下水位波動的敏感度系數擬合預測模型

圖8 地下水流速率敏感度系數分布

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

圖9 LNAPL滲透速率敏感度系數分布

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

圖10 地下水位波動敏感度系數分布

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

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

3 結語

（1）LNAPL 的NAPL 相達到穩定狀態時，地下水流速率與溶解速率呈現正相關趨勢；泄漏速率與NAPL相飽和度呈現反相關，NAPL相飽和度與溶解速率呈現正相關，NAPL 相飽和度對研究溶解速率變化起關鍵作用。

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

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