陶佳輝，施小清，康學遠，徐紅霞，吳吉春

(表生地球化學教育部重點實驗室/南京大學地球科學與工程學院，江蘇 南京 210023)

隨著社會對石油化工產(chǎn)品需求量的增加，大量以輕非水相液體(light non-aqueous phase liquid，LNAPL)形式存在的有機產(chǎn)品，在貯存及使用過程中發(fā)生泄漏，導致地下水污染[1～2]。一般而言，LNAPL密度比水小且難溶于水，泄漏后在重力與毛細壓力作用下運移[3]。若LNAPL泄漏量足夠大，少量以殘余態(tài)形式滯留在運移路徑中，部分向下運移直至遇到障礙(黏土層等)，形成輕油聚集體，其余LNAPL運移至毛細水帶和地下水位后產(chǎn)生側(cè)向運移[4]。通常將潛水面以上殘留LNAPL的區(qū)域稱為不飽和污染源區(qū)；將飽和帶或潛水面以下殘留LNAPL的區(qū)域稱為飽和污染源區(qū)[5～6]。研究者常采用LNAPL的飽和度分布(saturation distribution)、孔隙中殘余LNAPL體積(volume of possible LNAPL residence)、不連續(xù)離散狀(ganglia)與連續(xù)池狀(pool)LNAPL體積比(ganglia-to-pool ratio，GTP)等指標來定量表征LNAPL的污染源區(qū)結(jié)構(gòu)[7]。殘留在介質(zhì)中的LNAPL在自然條件下難以降解，從而成為長期存在的土壤-地下水污染源。因此，精細刻畫污染源區(qū)結(jié)構(gòu)對于LNAPL污染的風險評估和修復治理極其重要[8]。

國內(nèi)外學者主要通過試驗與數(shù)值模擬方法研究含水量、非均質(zhì)性、毛細壓力等因素對LNAPL運移和分布的影響[9～10]。如Schroth等[11]構(gòu)建了含有毛細障礙柵的砂箱實驗，研究了LNAPL在包氣帶非均質(zhì)介質(zhì)中的運移，結(jié)果表明含水量對LNAPL的運移和分布有較大影響。Wipfler等[12]通過二維砂箱實驗觀察LNAPL在層狀多孔介質(zhì)中的運移，結(jié)果表明毛細壓力是影響LNAPL運移和分布的主要敏感因子。束善治等[4]采用數(shù)值分析方法探究了不同特征的土層結(jié)構(gòu)對LNAPL運移和分布的影響。Van Geel等[13]通過兩相和三相試驗數(shù)據(jù)來預測LNAPL的殘余飽和度，結(jié)果表明LNAPL的殘余飽和度是最大流體飽和度和含水量的函數(shù)。

以往研究中通常視LNAPL從地表泄漏后穿過包氣帶運移至潛水面，前人的分析和探討大多建立在“LNAPL在潛水面漂浮移動”的污染概念模型上。然而，地下介質(zhì)的非均質(zhì)性和包氣帶含水量的空間變異性可能形成復雜的LNAPL污染源區(qū)結(jié)構(gòu)，LNAPL甚至可能無法到達潛水面，而在毛細水帶或局部透鏡體蓄積。因此，探討多種因素對復雜LNAPL污染源區(qū)結(jié)構(gòu)的影響具有重要意義。

本文基于數(shù)值模型研究LNAPL泄漏量、介質(zhì)非均質(zhì)性與含水量空間變異、潛水面周期性變化等因素對LNAPL污染源區(qū)結(jié)構(gòu)的影響。通過分析不同條件下LNAPL的飽和度分布來刻畫LNAPL污染源區(qū)結(jié)構(gòu)，以期為地下水中LNAPL污染的風險評估及修復治理提供理論參考。

1 計算方法

1.1 控制方程

本文采用TOUGH2軟件的T2VOC模塊模擬LNAPL的運移與分布[14]。TOUGH2可以解決不同條件下地下水中水流和熱量運移問題，而T2VOC模塊在模擬飽和帶、非飽和帶非水相液體運移的問題中應用廣泛。

T2VOC程序?qū)APL視為一個單獨的相，采用有限差分法進行空間離散，水、氣及NAPL中每一相的運移可根據(jù)Darcy定律的多相形式由壓力和重力確定[15]。在三相組分組成的非等溫系統(tǒng)中，完整地描述該系統(tǒng)需要三個質(zhì)量守恒方程和一個能量守恒方程。對于任意體積為Vn，表面積為Γn的滲流區(qū)，組分k的平衡方程為：

(1)

式中：Mk——單位體積組分k(k=w,a,n)的質(zhì)能累積量；

Fk——組分k進入體積Vn的質(zhì)能通量；

n——內(nèi)法線方向的單位向量；

qk——單位體積內(nèi)組分k的源匯項。

1.2 概念模型

假設模擬區(qū)域為二維XZ剖面，尺寸為90 cm×5 cm×90 cm，X方向Z方向均勻剖分為30×30共900個網(wǎng)格，網(wǎng)格水平與垂直間距均為3 cm。模型上邊界為大氣邊界，左右邊界為定水頭邊界，下邊界為隔水邊界(圖1)。甲苯密度比水小，是常見LNAPL污染物汽油中的典型組分，對人體與自然環(huán)境有較大危害，故本研究選取甲苯作為LNAPL的特征污染物[16]。污染源設于模擬區(qū)域頂端以下30 cm的中心處，泄漏速率恒為5×10-4kg/s，基于甲苯的運移和分布特征分析LNAPL污染源區(qū)結(jié)構(gòu)的影響因素。假設溫度恒為20 ℃，故未考慮熱物理相關(guān)參數(shù)。

相對滲透率函數(shù)采用Stone模型，Stone模型中Swr、Snr、Sgr分別表示液相、NAPL相以及氣相的殘余飽和度，n為擬合參數(shù)。毛管壓力函數(shù)采用Parker模型，Parker模型中Sm表示殘余液相飽和度,αgn和αnw分別表示氣相和NAPL相、NAPL相和水相之間進氣壓力的倒數(shù)。

圖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的運移和分布受泄漏量以及泄漏條件(事故性泄漏、多年持續(xù)泄漏等)的影響。隨著泄漏量的不斷增加，部分LNAPL蓄積在毛細水帶中，其余LNAPL將運移至潛水位以下[19]?；诓煌腖NAPL泄漏量，本節(jié)分析了泄漏量對LNAPL源區(qū)結(jié)構(gòu)的影響。模型介質(zhì)為均質(zhì)中砂，系統(tǒng)達到水-氣平衡后，以5×10-4kg/s的速度注入甲苯，注入時間分別為500 s、1 000 s和2 000 s，分析停止注入后靜置條件下LNAPL的運移和分布特征。模擬進行240 min后LNAPL的空間分布趨于穩(wěn)定，圖2為不同泄漏量下LNAPL飽和度分布圖。可見泄漏量較小時，LNAPL在毛細水帶中積蓄的壓力較小，幾乎所有的LNAPL都蓄積在毛細水帶中。隨著泄漏量的增大，LNAPL將穿越毛細水帶，少量LNAPL運移至潛水面以下。

2.2 非均質(zhì)性的影響

本節(jié)設計兩種不同的層狀非均質(zhì)模型，并與均質(zhì)中砂模型(圖3a)作對比。兩種層狀非均質(zhì)模型由兩層砂土構(gòu)成：上細下粗(圖3b)、上粗下細(圖3c)。假設甲苯的泄漏量較小，注入時間為500 s，三種模型的邊界條件、甲苯的注入速率和注入方式均相同。

系統(tǒng)在穩(wěn)定水流場中達到水-氣平衡后(不考慮降雨入滲)，介質(zhì)非均質(zhì)性會導致潛水面以上毛細水帶厚度的差異(圖4)。潛水位標高15 cm；潛水面以上，含水量接近飽和的區(qū)域為毛細水帶；毛細水帶上邊緣與大氣邊界之間的區(qū)域為土壤水帶。三種模型毛細水帶厚度分別為10，7，13 cm。上粗下細模型中毛細水帶厚度最大(圖4c)，原因在于：細砂層土顆粒粒徑較小，比表面積相對較大，產(chǎn)生較大的毛細壓力，使得細砂層對水有較強的吸附和滯留能力，通過自吸方式獲得較多水分。

圖2 不同泄漏量下LNAPL飽和度(SO)分布(CF：毛細水帶邊緣)Fig.2 LNAPL saturation(SO)distribution with different leakage (CF: Capillary Fringe)注：為顯示LNAPL主要分布區(qū)域，將SO=0.01作為飽和度閾值，下同。

圖3 均質(zhì)和層狀非均質(zhì)概念模型Fig.3 Conceptual model showing the homogeneous and layered heterogeneity注：實線代表細砂/粗砂、粗砂/細砂分界線，下同。

圖4 不同介質(zhì)條件下水的飽和度(Sw)分布Fig.4 Water saturation(Sw)distribution under different medium situations注：虛線代表地下水位線和毛細水帶上邊緣，下同。

圖5為三種不同介質(zhì)中不同時刻LNAPL飽和度分布圖。均質(zhì)介質(zhì)中(圖5a)，滲流初始階段(T=10 min)，LNAPL受重力作用影響較大，以垂向入滲為主。T=60 min時，LNAPL受毛細壓力影響較大，開始以橫向遷移為主，并在毛細水帶邊緣達到最大飽和度0.38。模擬結(jié)束時大部分LNAPL進入到毛細水帶中，但由于LNAPL密度比水小且泄漏量較小，最終滯留在毛細水帶中。

上細下粗模型中(圖5b)，T=30 min時，LNAPL開始在細砂層下部蓄積，T=60 min時，細砂層中蓄積的LNAPL達到最大飽和度0.38。模擬結(jié)束時大部分LNAPL穿越細砂層進入到毛細水帶中，但仍有部分LNAPL滯留在細砂層無法下滲。

上粗下細模型中(圖5c)，滲流初始階段由于粗砂層較大的滲透率和顆粒粒徑，LNAPL在粗砂較快地運移至毛細水帶邊緣，隨后在毛細水帶邊緣蓄積并橫向遷移。模擬結(jié)束時LNAPL均勻水平地展布在毛細水帶上邊緣，幾乎沒有LNAPL進入到毛細水帶中。

圖5 三種介質(zhì)中不同時刻LNAPL飽和度(SO)分布Fig.5 LNAPL saturation(SO)distribution at different times in three kinds of media

圖6為模擬結(jié)束時三種介質(zhì)中水和LNAPL的飽和度分布隨深度變化的剖面圖。上細下粗模型中(圖6b)細砂層阻礙LNAPL的垂向入滲，最終仍有部分LNAPL(最大飽和度達0.22)滯留于細砂層無法向下運移，成為長期污染源。原因是細砂層固有滲透率低，LNAPL在細砂中運移較為困難；細砂層部分孔隙體積被水占有，含水量較大，LNAPL的相對滲透率較低[20]；細砂層土顆粒粒徑較小，較大的毛細壓力驅(qū)使LNAPL產(chǎn)生較多的橫向遷移，阻礙LNAPL的垂向入滲。

上粗下細模型中(圖6c)，LNAPL最終蓄積在細砂層毛細水帶邊緣。LNAPL無法進入到細砂層中的毛細飽和帶，原因是細砂層土顆粒粒徑較小，LNAPL相對滲透率較低，不容易發(fā)生運移；細砂層毛細水帶含水量較高，對LNAPL產(chǎn)生向上的浮力作用，阻礙LNAPL的垂向入滲；細砂層的進氣壓力較大，LNAPL與水之間的壓力差不足以克服細砂顆粒對水的毛細束縛，最終在毛細水帶邊緣蓄積[21]。

一般認為LNAPL在包氣帶中可直接運移至潛水面，然而圖6(b)、6(c)顯示，上細下粗非均質(zhì)條件下，部分LNAPL滯留在細砂層中無法下滲，形成雙重LNAPL高飽和度區(qū)；上粗下細條件下幾乎所有的LNAPL均滯留在毛細水帶邊緣，不會運移至潛水面。兩者產(chǎn)生差異的原因在于：包氣帶中LNAPL的運移和分布是油-水-氣三相相互驅(qū)替的過程，水相對于空氣具有更大的密度和黏度，LNAPL到達毛細水帶后更容易驅(qū)替空氣向四周遷移，并在毛細水帶處發(fā)生蓄積[22]。由于LNAPL相對空氣而言是浸潤相，而相對水而言是非浸潤相，且細砂層固有滲透率較低，進氣壓力較大，阻礙LNAPL的下滲，因此最終其無法穿越細砂層的毛細水帶。介質(zhì)非均質(zhì)性將使LNAPL污染源區(qū)結(jié)構(gòu)趨于復雜，增大污染場地修復治理的難度[23]。

圖6 水的飽和度(Sw)和LNAPL飽和度(SO)分布對比Fig.6 Comparison of water saturation(Sw)and LNAPL saturation(SO)distribution

2.3 含水量與非均質(zhì)性的共同影響

圖7 含透鏡體概念模型及水的飽和度(Sw)分布Fig.7 Conceptual model showing the clay lens and water saturation(Sw)distribution

本節(jié)在均質(zhì)中砂背景中設置黏土透鏡體(圖7a)，分別考慮飽和與非飽和兩種不同含水量情形。為避免毛細水帶對結(jié)果的影響，將模擬范圍擴大為原來的100倍，同時將地下水位降為10 m，這樣可使得LNAPL運移區(qū)域遠離毛細水帶。飽和透鏡體模型考慮降雨，系統(tǒng)首先在持續(xù)降水入滲條件下達到平衡，然后停止降雨入滲，以5×10-4kg/s的速度注入甲苯，分析持續(xù)注入1 a過程中甲苯的運移和分布情況。系統(tǒng)在降雨條件下達到水-氣平衡后，背景中砂的含水飽和度(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均在黏土透鏡體上方蓄積并側(cè)向運移形成污染池，使得水平擴散范圍顯著增大。兩種情形僅含水量存在差異，但二者LNAPL污染源區(qū)結(jié)構(gòu)存在較大差異。污染土體中一般同時存在氣相、水相和NAPL相，當LNAPL蓄積的壓力超過LNAPL/水界面的“進氣壓力”時，LNAPL才能驅(qū)替孔隙中的水相或氣相[4,25～26]。由于飽和透鏡體含水量較高，毛細壓力較大，從模擬結(jié)果來看，LNAPL無法驅(qū)替飽和透鏡體中的水分和空氣，飽和黏土透鏡體對于LNAPL來說是不可滲的(圖8a)。而在非飽和透鏡體情形下，LNAPL可以驅(qū)替非飽和透鏡體中的水分和空氣，并穿透非飽和透鏡體繼續(xù)向下入滲(圖8b)。由此可見包氣帶中黏土透鏡體并非都是LNAPL運移的阻礙，當黏土透鏡體中含水量低時，LNAPL可以穿透黏土透鏡體；而黏土透鏡體含水量較高時，將阻礙LNAPL的入滲。

2.4 潛水面周期性變化的影響

地下介質(zhì)中的LNAPL主要分為自由相、殘留相與溶解相三種[27～28]。由于LNAPL多集中分布于毛細水帶及地下水位波動帶，不同相態(tài)之間的分配轉(zhuǎn)化主要是由潛水位波動造成的[29～30]。鑒于此，本節(jié)在均質(zhì)中砂模型的基礎上，將最下層30個網(wǎng)格作為注水和抽水的水位波動控制單元格，通過改變其源匯項條件來升高/降低地下水位(h)，水位先上升至0.3 m，再下降至0 m，分析地下水位一個升降周期內(nèi)LNAPL源區(qū)結(jié)構(gòu)特征的改變。

圖9為地下水位波動時LNAPL飽和度分布圖，可見水位波動使LNAPL污染范圍變大。波動初期水位上升至0.3 m(圖9b)，由于LNAPL密度比水小，水位上升過程中油-水-氣三相相互驅(qū)替，水位上升產(chǎn)生的浮力使LNAPL一同抬升。與水位波動前相比，被攜帶上升的LNAPL的橫向分布范圍明顯擴大。水位下降至0 m時(圖9c)，隨著水位下降LNAPL受重力作用垂向入滲。水位下降過程中，重力使LNAPL向下運移使得垂向分布范圍變大，又由于LNAPL向下驅(qū)替水比較困難，而向四周驅(qū)替空氣相對容易，故LNAPL橫向分布范圍變大，最終LNAPL的分布范圍逐漸擴大至整個水位波動帶[31]。

地下水位波動導致土壤含水量以及其他環(huán)境因素的改變，水位波動帶內(nèi)油-水-氣三相相互作用促進了地下NAPL相的重新分配[32]。水位上升會為LNAPL相提供浮力，驅(qū)替土顆粒內(nèi)的LNAPL相從孔隙中漂浮出來；水位下降時LNAPL受重力作用向下層運移，使得污染范圍明顯增大。附著于土壤的LNAPL隨地下水位波動而不斷地被溶解和裹挾，使污染源區(qū)結(jié)構(gòu)趨于復雜。

圖9 地下水位波動時LNAPL飽和度(SO)分布Fig.9 LNAPL saturation(SO)distribution under the water table fluctuations

3 結(jié)論

(1)當泄漏量較大時，LNAPL可運移至潛水面。

(2)當泄漏量較小時，均質(zhì)條件下，LNAPL主要在潛水面上方的毛細水帶滯留；上細下粗條件下，LNAPL在土壤水帶及毛細水帶均會滯留；上粗下細條件下，LNAPL主要在毛細水帶上邊緣發(fā)生蓄積，無法到達潛水面。

(3)包氣帶中黏土透鏡體并非都是LNAPL運移的阻礙，LNAPL可以穿透低含水量的黏土透鏡體，只有高含水量的黏土透鏡體才對LNAPL的入滲有阻礙作用。

(4)潛水面周期性變化導致地下LNAPL分布范圍變大，擴大污染范圍。