輕非水相液體污染源區(qū)結(jié)構(gòu)的影響因素?cái)?shù)值分析

陶佳輝，施小清，康學(xué)遠(yuǎn)，徐紅霞，吳吉春

(表生地球化學(xué)教育部重點(diǎn)實(shí)驗(yàn)室/南京大學(xué)地球科學(xué)與工程學(xué)院，江蘇 南京 210023)

隨著社會(huì)對(duì)石油化工產(chǎn)品需求量的增加，大量以輕非水相液體(light non-aqueous phase liquid，LNAPL)形式存在的有機(jī)產(chǎn)品，在貯存及使用過(guò)程中發(fā)生泄漏，導(dǎo)致地下水污染[1～2]。一般而言，LNAPL密度比水小且難溶于水，泄漏后在重力與毛細(xì)壓力作用下運(yùn)移[3]。若LNAPL泄漏量足夠大，少量以殘余態(tài)形式滯留在運(yùn)移路徑中，部分向下運(yùn)移直至遇到障礙(黏土層等)，形成輕油聚集體，其余LNAPL運(yùn)移至毛細(xì)水帶和地下水位后產(chǎn)生側(cè)向運(yùn)移[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)等指標(biāo)來(lái)定量表征LNAPL的污染源區(qū)結(jié)構(gòu)[7]。殘留在介質(zhì)中的LNAPL在自然條件下難以降解，從而成為長(zhǎng)期存在的土壤-地下水污染源。因此，精細(xì)刻畫污染源區(qū)結(jié)構(gòu)對(duì)于LNAPL污染的風(fēng)險(xiǎn)評(píng)估和修復(fù)治理極其重要[8]。

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

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

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

1 計(jì)算方法

1.1 控制方程

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

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

(1)

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

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

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

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

1.2 概念模型

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

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

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

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

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

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

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

圖2 不同泄漏量下LNAPL飽和度(SO)分布(CF：毛細(xì)水帶邊緣)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注：實(shí)線代表細(xì)砂/粗砂、粗砂/細(xì)砂分界線，下同。

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

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

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

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

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

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

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

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

圖6 水的飽和度(Sw)和LNAPL飽和度(SO)分布對(duì)比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ì)中砂背景中設(shè)置黏土透鏡體(圖7a)，分別考慮飽和與非飽和兩種不同含水量情形。為避免毛細(xì)水帶對(duì)結(jié)果的影響，將模擬范圍擴(kuò)大為原來(lái)的100倍，同時(shí)將地下水位降為10 m，這樣可使得LNAPL運(yùn)移區(qū)域遠(yuǎn)離毛細(xì)水帶。飽和透鏡體模型考慮降雨，系統(tǒng)首先在持續(xù)降水入滲條件下達(dá)到平衡，然后停止降雨入滲，以5×10-4kg/s的速度注入甲苯，分析持續(xù)注入1 a過(guò)程中甲苯的運(yùn)移和分布情況。系統(tǒng)在降雨條件下達(dá)到水-氣平衡后，背景中砂的含水飽和度(Sw=15%～19%)僅略高于殘余飽和度而黏土透鏡體基本達(dá)到飽和(圖7b，Sw>98%)[24]。非飽和透鏡體模型不考慮降雨條件，黏土透鏡體中含水量為其殘余飽和度(圖7c)。

圖8 含透鏡體模型不同時(shí)刻LNAPL飽和度(SO)分布Fig.8 LNAPL saturation(SO)distribution at different times in the clay lens model

圖8為含透鏡體模型不同時(shí)刻LNAPL飽和度分布圖。在飽和與非飽和兩種不同含水量情形下，LNAPL均在黏土透鏡體上方蓄積并側(cè)向運(yùn)移形成污染池，使得水平擴(kuò)散范圍顯著增大。兩種情形僅含水量存在差異，但二者LNAPL污染源區(qū)結(jié)構(gòu)存在較大差異。污染土體中一般同時(shí)存在氣相、水相和NAPL相，當(dāng)LNAPL蓄積的壓力超過(guò)LNAPL/水界面的“進(jìn)氣壓力”時(shí)，LNAPL才能驅(qū)替孔隙中的水相或氣相[4,25～26]。由于飽和透鏡體含水量較高，毛細(xì)壓力較大，從模擬結(jié)果來(lái)看，LNAPL無(wú)法驅(qū)替飽和透鏡體中的水分和空氣，飽和黏土透鏡體對(duì)于LNAPL來(lái)說(shuō)是不可滲的(圖8a)。而在非飽和透鏡體情形下，LNAPL可以驅(qū)替非飽和透鏡體中的水分和空氣，并穿透非飽和透鏡體繼續(xù)向下入滲(圖8b)。由此可見(jiàn)包氣帶中黏土透鏡體并非都是LNAPL運(yùn)移的阻礙，當(dāng)黏土透鏡體中含水量低時(shí)，LNAPL可以穿透黏土透鏡體；而黏土透鏡體含水量較高時(shí)，將阻礙LNAPL的入滲。

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

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

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

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

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

3 結(jié)論

(1)當(dāng)泄漏量較大時(shí)，LNAPL可運(yùn)移至潛水面。

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

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

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