馬騰飛，劉漢樂(lè)，張 閃

(1.桂林理工大學(xué)環(huán)境科學(xué)與工程學(xué)院，廣西 桂林541004;2.廣西環(huán)境工程與保護(hù)評(píng)價(jià)重點(diǎn)實(shí)驗(yàn)室，廣西 桂林541004)

層狀非均質(zhì)多孔介質(zhì)中LNAPLs污染物運(yùn)移與分布的鋒面擴(kuò)展模型

馬騰飛1，2，劉漢樂(lè)1，2，張 閃2

(1.桂林理工大學(xué)環(huán)境科學(xué)與工程學(xué)院，廣西 桂林541004;2.廣西環(huán)境工程與保護(hù)評(píng)價(jià)重點(diǎn)實(shí)驗(yàn)室，廣西 桂林541004)

通過(guò)建立二維物理模型，研究非飽和層狀非均質(zhì)多孔介質(zhì)中輕非水相液體LNAPLs(Light Non－aqueous Phase Liquids)的入滲機(jī)制與變化特征，建立LNAPLs運(yùn)移與分布的鋒面擴(kuò)展模型，探討入滲階段油流鋒面擴(kuò)展速率的變化規(guī)律。結(jié)果發(fā)現(xiàn):模型中砂土的滲透系數(shù)、孔隙度及油的相對(duì)滲透系數(shù)對(duì)LNAPLs運(yùn)移的鋒面擴(kuò)展速率影響較大，而油的飽和度與壓頭變化對(duì)其影響較小。當(dāng)LNAPLs由中砂進(jìn)入細(xì)砂與細(xì)砂傾斜透鏡體時(shí)，砂的滲透系數(shù)、孔隙度及油的相對(duì)滲透系數(shù)減小，這種介質(zhì)結(jié)構(gòu)面的突變改變了 LNAPLs鋒面擴(kuò)展速率，其在中砂中鋒面擴(kuò)展速率快而在細(xì)砂與細(xì)砂傾斜透鏡體中鋒面擴(kuò)展速率慢。此外，當(dāng)LNAPLs進(jìn)入干濕界面后，由于毛細(xì)作用的增強(qiáng)與含水量的加大，其橫向比垂向運(yùn)移的平均鋒面擴(kuò)展速率大。

輕非水相液體;鋒面擴(kuò)展模型;多孔介質(zhì);運(yùn)移與分布;毛細(xì)作用

近些年，大量的加油站、地下儲(chǔ)油設(shè)施使用過(guò)程中石油產(chǎn)品的泄露給土壤和地下環(huán)境造成了嚴(yán)重的污染。石油產(chǎn)品中的非水相液體NAPLs(Non－aqueous phase liquids)在水中的溶解度低，其中密度比水小的稱為輕非水相液體LNAPLs(Light NAPLs)，多為揮發(fā)性有機(jī)物，常見(jiàn)的有汽油、煤油等［1，2］。這些污染物與水不混溶，進(jìn)入地下環(huán)境后會(huì)進(jìn)一步擴(kuò)大污染范圍，使危害進(jìn)一步加劇。因此，對(duì)其污染與分布、控制與治理等方面的研究已成為當(dāng)前地下水環(huán)境領(lǐng)域的前沿課題。

國(guó)外在該方面的研究已開(kāi)展多年，如Pantazidou和Sitar［3］(1993)用二維物理試驗(yàn)?zāi)M了 LNAPLs在非飽和帶的運(yùn)移，并建立相應(yīng)的數(shù)學(xué)模型預(yù)測(cè)LNAPLs在入滲期間的傳播速率和油透鏡體達(dá)到的最終厚度;Schroth［4］(1995)用二維砂槽的方法預(yù)測(cè)了多孔介質(zhì)中LNAPLs透鏡體的形狀和厚度;Oostrom［5］等(1997)研究了飽和砂土中LNAPLs運(yùn)移的特性，從汽油溢漏區(qū)域獲得數(shù)據(jù)，應(yīng)用基本理論建立模型;E.L.Wipfler［6］(2004)通過(guò)兩個(gè)對(duì)比試驗(yàn)探索了 LNAPLs運(yùn)移和分布的規(guī)律，結(jié)果發(fā)現(xiàn)毛細(xì)作用對(duì)LNAPLs的運(yùn)移起決定因素。相對(duì)而言，國(guó)內(nèi)在該方面的研究起步較晚，武曉鋒［7］(2000)在二維砂槽模型中模擬了LNAPLs在均勻多孔介質(zhì)地下水非飽和區(qū)和飽和區(qū)中的運(yùn)移過(guò)程，并利用多孔介質(zhì)的毛細(xì)管模型估算LNAPLs透鏡體的厚度，但未考慮介質(zhì)層狀非均質(zhì)性對(duì)LNAPLs運(yùn)移的影響;束善治［8］(2002)采用數(shù)值分析的方法，對(duì)LNAPLs在包氣帶非均質(zhì)介質(zhì)中的遷移和分布特性進(jìn)行了分析研究，但并未與真實(shí)的試驗(yàn)結(jié)果相對(duì)比，無(wú)從驗(yàn)證數(shù)值模型的有效性;劉漢樂(lè)［9］(2006)通過(guò)室內(nèi)試驗(yàn)建立物理模型，研究了非飽和帶中層狀非均質(zhì)條件下LNAPLs的運(yùn)移機(jī)制和分布特性，但并未建立相應(yīng)的數(shù)學(xué)模型。

Pantazidou和Sitar［3］(1993)建立一個(gè)單相流的數(shù)學(xué)模型，描述了LNAPLs在均質(zhì)條件下的入滲鋒面擴(kuò)展情形。鄭德鳳［10］(2002)用二維砂箱模擬了LNAPLs在濕潤(rùn)多孔介質(zhì)中的泄漏，得出LNAPLs在均質(zhì)條件下的入滲、遷移及分布規(guī)律，在實(shí)驗(yàn)基礎(chǔ)上利用此模型，預(yù)測(cè)滲流帶中污染鋒面的擴(kuò)展速度及透鏡體的最終形狀，而將該模型建立在非均質(zhì)條件下模擬LNAPLs運(yùn)移與分布的研究還未見(jiàn)報(bào)道。因此本文通過(guò)室內(nèi)試驗(yàn)?zāi)MLNAPLs入滲鋒面的變化特征，在層狀非均質(zhì)條件下運(yùn)用單相流模型的簡(jiǎn)化方程，建立LNAPLs運(yùn)移與分布的鋒面擴(kuò)展模型，探索油流鋒面擴(kuò)展速率的變化規(guī)律，對(duì)理解LNAPLs在層狀非均質(zhì)條件下的運(yùn)移機(jī)制具有重要意義。

1 試驗(yàn)材料及方法

試驗(yàn)裝置為無(wú)色有機(jī)玻璃制成的二維砂槽，尺寸為L(zhǎng)×H×W=112×90×4 cm，兩邊是水室，用于控制地下水位。水室與砂槽之間的隔板上打有小孔并貼有金屬絲網(wǎng)。砂槽中裝入三種粒徑不同的石英砂，用于模擬層狀非均質(zhì)多孔介質(zhì)。細(xì)砂粒徑為 0.1 ～0.2 mm，中砂粒徑為 0.2 ～0.4 mm，粗砂粒徑為0.4～0.5 mm。砂槽主要部分為中砂，其中設(shè)置有細(xì)砂、粗砂水平透鏡體和細(xì)砂傾斜透鏡體，如圖1所示。

圖1 LNAPLs入滲鋒面變化圖

試驗(yàn)采用汽油作為L(zhǎng)NAPL的代表，用溶于有機(jī)溶劑而不溶于水的蘇丹紅將LNAPLs染色，以提高LNAPLs與背景之間的對(duì)比度。試驗(yàn)過(guò)程中，以一定的時(shí)間間隔用數(shù)碼相機(jī)進(jìn)行拍照，把LNAPLs的擴(kuò)散過(guò)程以圖像的形式記錄下來(lái)。注油試驗(yàn)從開(kāi)始至結(jié)束歷經(jīng)132 min，試驗(yàn)過(guò)程中室內(nèi)溫度為19℃ ～20℃。試驗(yàn)結(jié)束后，對(duì)圖片進(jìn)行分析，并用AutoCAD軟件進(jìn)行處理，然后將不同時(shí)刻LNAPLs的鋒面重疊到一幅圖上，合成不同時(shí)刻LNAPLs鋒面的空間分布變化圖，見(jiàn)圖1。

2 預(yù)測(cè)模型

2.1 模型的建立

根據(jù)Pantazidou和Sitar［6］(1993)描述入滲期間的鋒面擴(kuò)展情形，可把LNAPLs在非均質(zhì)多孔介質(zhì)中入滲期間的鋒面曲線近似為橢圓。當(dāng)污染質(zhì)從一污染源向下滲漏并在淺層滲流帶中擴(kuò)散，對(duì)三維情況，用直徑為的球體來(lái)近似油流體積，在時(shí)間內(nèi)油流向前移動(dòng)的距離。其縱剖面如圖2所示。

圖2 油入滲鋒面理想模式

污染源用一面積為的點(diǎn)源近似，油經(jīng)過(guò)剖面為的限定面積向砂樣滲透，用達(dá)西定律描述為:

式中:dq/dt為油的入滲速率(cm3/min)，k為砂的滲透系數(shù)，kr0為油的相對(duì)滲透系數(shù)，ρ0油的密度，g為重力加速度，u0為油的粘度，A為污染源的剖面積，L為油進(jìn)入砂樣前的垂向運(yùn)移距離，ΔP為經(jīng)過(guò)距離的壓頭變化。于是在dt時(shí)間內(nèi)，進(jìn)入到介質(zhì)中的油體積為:

假設(shè)油流污染范圍內(nèi)油的平均飽和度為S0，則進(jìn)入到介質(zhì)中的油體積也等于dt時(shí)間內(nèi)鋒面擴(kuò)展的污染面積內(nèi)孔隙中的油體積

式中:n為介質(zhì)孔隙度，S0為污染范圍內(nèi)油的平均飽和度，x、x+dx分別代表t、t+Δt時(shí)刻油流的垂向高度。將(2)式與(3)式合并，忽略微分高次冪項(xiàng)并積分可得

污染源若用一個(gè)面積A為的點(diǎn)源來(lái)近似，則w與此滲漏面積的直徑相對(duì)應(yīng);若用一線源近似，則w為線源長(zhǎng)度。對(duì)二維情形來(lái)說(shuō)，用一個(gè)直徑為x，長(zhǎng)度為w水平放置的圓柱體來(lái)近似油流體積?？捎门c上面類似的方法得出

式(4)和式(5)描繪了油流鋒面隨時(shí)間的擴(kuò)展情形，根據(jù)(5)式把作為油流鋒面運(yùn)移速率的特征比值。

2.2 模型參數(shù)的確定

根據(jù)介質(zhì)的特性［11］，模型與試驗(yàn)中相關(guān)參數(shù)的選擇如表1所示。

表1 模型輸入?yún)?shù)

表2 砂槽中相關(guān)參數(shù)的變化

3 結(jié)果分析與討論

3.1 模型預(yù)測(cè)結(jié)果

在LNAPLs由中砂進(jìn)入傾斜細(xì)砂透鏡體時(shí)，油的相對(duì)滲透系數(shù)、砂樣的滲透系數(shù)及孔隙度各為中砂與細(xì)砂層中的1/2，下表列出了計(jì)算式(5)中所包含的輸入?yún)?shù)值，模型預(yù)測(cè)的結(jié)果和相應(yīng)的試驗(yàn)測(cè)量結(jié)果以及相對(duì)誤差，如表3、表4所示。

表3 模型垂向鋒面擴(kuò)展速率的計(jì)算值與實(shí)測(cè)值

表4 模型橫向鋒面擴(kuò)展速率的計(jì)算值與實(shí)測(cè)值

表5 垂向與橫向鋒面擴(kuò)展的平均速率

3.2 砂體中 LNAPLs的運(yùn)移分析

入滲初始，LNAPLs主要在重力作用下垂直向下遷移，同時(shí)在毛細(xì)力作用下產(chǎn)生橫向遷移，但橫向遷移不顯著。LNAPLs發(fā)生遷移時(shí)鋒面線形狀不規(guī)則，呈現(xiàn)出指狀擴(kuò)散，這是由于介質(zhì)局部微觀存在非均質(zhì)性，流體在多孔介質(zhì)中某一部分運(yùn)移的速率比平均鋒面的速率大(如第10～13min大1.19、第54～61min大 0.45、第61～74 min大0.07與第102～112 min大0.77)，即發(fā)生突進(jìn)部分，導(dǎo)致界面移動(dòng)的不規(guī)則。第7～10 min與第10～13 min LNAPLs在中砂體中運(yùn)移，鋒面擴(kuò)展速率由慢到快，在第13 min LNAPLs到達(dá)第一個(gè)細(xì)砂層上界面時(shí)，鋒面擴(kuò)展速率開(kāi)始變慢。這是由于細(xì)砂層中砂的滲透系數(shù)、油的相對(duì)滲透系數(shù)及孔隙度比中砂中小(表2)，其對(duì)LNAPLs運(yùn)移的鋒面擴(kuò)展速率影響較大，而油的飽和度與壓頭變化(表3)對(duì) LNAPLs的運(yùn)移影響較小，砂槽中這種介質(zhì)結(jié)構(gòu)面的突變使 LNAPLs沿著細(xì)砂透鏡體產(chǎn)生橫向運(yùn)移，縱向運(yùn)移受到阻礙。

第21 min LNAPLs到達(dá)第一個(gè)細(xì)砂層下界面，隨著時(shí)間的推移，LNAPLs運(yùn)移的鋒面擴(kuò)展速率再次加快。由于砂土裝填不均勻，LNAPLs沒(méi)有同時(shí)進(jìn)入左右兩側(cè)的細(xì)砂與粗砂透鏡體。第61～74 min LNAPLs逐漸由中砂層進(jìn)入傾斜細(xì)砂透鏡體，此時(shí)鋒面擴(kuò)展速率又開(kāi)始變慢。這是由于傾斜細(xì)砂透鏡體中介質(zhì)特性，如砂的滲透系數(shù)、孔隙度、油的相對(duì)滲透系數(shù)等與中砂中相比較小(表2)，并且毛細(xì)作用作為 LNAPLs橫向運(yùn)移的驅(qū)動(dòng)力，使其縱向運(yùn)移的速率變慢，開(kāi)始產(chǎn)生橫向運(yùn)移。

第86～94 min LNAPLs到達(dá)干濕界面，其鋒面擴(kuò)展速率變慢，持續(xù)到第102 min。這是由于隨著含水量的加大，壓頭變化相對(duì)較小，由16 cm逐漸減小到8 cm(表3)，此外毛細(xì)作用力也阻礙了 LNAPLs的縱向運(yùn)移。第102～112 min LNAPLs隨著壓頭變化的增大，由8 cm增大到12 cm，其鋒面擴(kuò)展速率加快。

第112～185 min，LNAPLs以橫向運(yùn)移的形式在地下水位上部聚積，其橫向運(yùn)移的平均速率比垂向運(yùn)移的平均速率大(表5)，說(shuō)明LNAPLs進(jìn)入干濕界面后，其與水之間的壓力差受到顆?？紫吨兴拿?xì)束縛，LNAPLs不能驅(qū)替孔隙水垂直下滲，于是產(chǎn)生側(cè)向傳播。當(dāng) LNAPLs聚積到一定厚度時(shí)，足以克服孔隙水的毛細(xì)束縛到達(dá)地下水位，在地下水位上產(chǎn)生 LNAPLs池。

3.3 誤差分析

模型中的誤差主要由以下幾方面引起:(1)第41～54 min的相對(duì)誤差比第7～10 min、第10～13 min與第13～21 min的相對(duì)誤差要大(表3中的1、2、3、4)，這是由于在第41～54 min鋒面線的兩側(cè)分別設(shè)置了細(xì)砂與粗砂透鏡體，砂土裝填時(shí)會(huì)造成介質(zhì)的不均勻。當(dāng)LNAPLs在砂槽中運(yùn)移時(shí)會(huì)受到大顆粒的阻滯作用，這種人為因素產(chǎn)生的誤差是不可避免的。(2)在 LNAPLs到達(dá)干濕界面后(表3中的7、8、9)，由于含水量的加大模型中油的飽和度及壓頭變化減小，同時(shí)受到孔隙水毛細(xì)作用力的束縛，第86～94 min與第94～102 min的相對(duì)誤差較大。(3)介質(zhì)結(jié)構(gòu)面的突變會(huì)產(chǎn)生一定的誤差(表3中的3、6)，在 LNAPLs由中砂進(jìn)入細(xì)砂與細(xì)砂透鏡體時(shí)，砂的滲透系數(shù)、孔隙度及油的相對(duì)滲透系數(shù)、飽和度等受到介質(zhì)特性的影響發(fā)生變化(表2)，使得計(jì)算值與實(shí)測(cè)值之間存在誤差。(4)模型中的參數(shù)如砂的滲透系數(shù)、孔隙度，油的飽和度及相對(duì)滲透系數(shù)等，與真實(shí)介質(zhì)的特性間存在一定差異，因此 LNAPLs鋒面擴(kuò)展速率的計(jì)算值與實(shí)測(cè)值間的誤差是不可避免的。

由預(yù)測(cè)模型得到的不同時(shí)間 LNAPLs鋒面擴(kuò)展速率計(jì)算值與實(shí)測(cè)值之間的誤差可看出，除了少數(shù)幾個(gè)相對(duì)誤差偏小外，多數(shù)的相對(duì)誤差都偏大。說(shuō)明計(jì)算值與實(shí)測(cè)值之間存在較大的差異，表明利用該模型模擬 LNAPLs在非飽和非均質(zhì)條件下的運(yùn)移仍需進(jìn)一步研究。

4 結(jié)語(yǔ)

通過(guò)試驗(yàn)結(jié)果的分析與討論，可得出以下結(jié)論:

1)本文通過(guò)試驗(yàn)?zāi)M了LNAPLs在非飽和層狀非均質(zhì)條件下的運(yùn)移，得到 LNAPLs鋒面的空間分布變化圖，圖像清晰地反映了 LNAPLs的入滲過(guò)程。在試驗(yàn)的基礎(chǔ)上，利用均勻介質(zhì)條件下單向流模型的簡(jiǎn)化方程，建立 LNAPLs在非飽和層狀非均質(zhì)條件下的鋒面擴(kuò)展模型，獲得其運(yùn)移的鋒面擴(kuò)展速率。

2)通過(guò)誤差分析看出，由于人為因素、介質(zhì)自身特性間的差異，如模型中砂的滲透系數(shù)、孔隙度及油的相對(duì)滲透系數(shù)、飽和度等對(duì)LNAPLs運(yùn)移的鋒面擴(kuò)展速率影響較大，這種不均勻性及介質(zhì)結(jié)構(gòu)面的突變都會(huì)產(chǎn)生一定的誤差。

3)通過(guò)模型獲得LNAPLs運(yùn)移的鋒面擴(kuò)展速率預(yù)測(cè)值與實(shí)測(cè)值的結(jié)果顯示，在非飽和層狀非均質(zhì)下該模型還難以應(yīng)用，有待建立更可靠的預(yù)測(cè)模型。此外，文中僅模擬了 LNAPLs在室內(nèi)小尺度上的運(yùn)移，而在野外大尺度條件下獲取LNAPLs的入滲鋒面變化圖，并建立相應(yīng)的數(shù)學(xué)模型對(duì)其污染范圍進(jìn)行預(yù)測(cè)、圈定仍有待進(jìn)一步研究。

［1］Borden R C，Gomez C A，Becker M T.Geochemical indicators of intrinsic bioremediation［J］.Ground Water，1995，33(2):180 －189.

［2］李永濤，王文科，王麗等.輕非水相液體在地下多相流的K－SP關(guān)系研究［J］.灌溉排水學(xué)報(bào).2008，(06):1672－3317.

［3］Pantazidou M，Sitar N Emplacement of nonaqueous liquids in the vadose zone［J］.Water Resources，1993，(29):705 － 721.

［4］Schroth M H ，Istok J D ，Ahearn S J，et al.Geometry and position of light nonaqueous－phase liquid lenses in water－wetted porous media［J］.Journal of Contaminant Hydrology ，1995，(19):269 －287.

［5］Oostrom，M.，C.Hofstee，and J.H.Dane，Light nonaqueousphase liquid movement in a variably saturated sand，Soil Sci.Soc.Am.J.，1997，61:1547 －1554.

［6］E.L.Wipfler.Infiltration and redistribution of LNAPL into unsatured layered porous media［J］.Journal of Contaminant Hydrlogy，2004，(7):47－19.

［7］武曉鋒，唐杰，藤間幸久.地下水中輕質(zhì)有機(jī)污染物(LNAPL)透鏡體研究［J］.環(huán)境污染與防治.2000，(3):17－5.

［8］束善治，梁宏偉，袁勇.輕非水相液體在非均質(zhì)地層包氣帶中運(yùn)移和分布特征數(shù)值分析［J］.水利學(xué)報(bào).2002，(11):31－7.

［9］劉漢樂(lè)，周啟友，徐速.非飽和帶中非均質(zhì)條件下 LNAPL運(yùn)移與分布特性實(shí)驗(yàn)研究［J］.水文地質(zhì)工程地質(zhì).2006，(5):52－56.

［10］鄭德鳳，趙勇勝，王本德.輕非水相液體在地下環(huán)境中的運(yùn)移特征與模擬預(yù)測(cè)研究［J］.水科學(xué)進(jìn)展.2002，(3):321 －5.

［11］王大純等.水文地質(zhì)學(xué)基礎(chǔ)［M］.北京出版社.1995，6.

［12］White M D ，Ostrom M，Lenhard R J.Modeling fluid flow and transport in variably saturated porous media with the STOMP simulator.1.Nonvolatile three － phase model description［J］.Advances in Water Resources，1995，18(6):353 － 364.

［13］Lenhard R J，Ostrom M，White M D.Modeling fluid flow and transport in variably saturated porous media with the STOMP simulator.2.Nonvolatile three － phase model description［J］.Advances in Water Resources，1995 ，18(6):365 – 373.

［14］M.Y.Corapcioglu etal.Pore－ Scale Analysis of NAPL Blob Dissolution and Mobilization in Porous Media［J］.Transp Porous Med(2009)79:419－442.

Infiltration＆Distribution front Spreading Model of LNAPLs in Layered Heterogenious Porous Media

MA Teng － fei1，2，LIU Han － le1，2，ZHANG Shan2
(1.School of Environment Sciences and Engineering，Guilin University of Technology，Guilin 541004，Guangxi;2.The Guangxi Key Laboratory of Environmental Engineering，Protection and Assessment，Guilin 541004，Guangxi)

The paper is based on a 2－D physical model，the purpose is to explore characteristics of LNAPLs infiltration in unsaturated layered heterogeneous porous media，and establish the inflitration and distribution of LNAPLs front spreading model，forecast the characteristic ratio of oil spreading rate in infiltration stage.The result show that capillary，porosity of sand and relative capillary of oil in model affect LNAPLs front spreading rate，but saturation of oil and the pressure head changed influence is small.When LNAPLs transfer from medium sand to fine，the capillary、porosity of sand and relative capillary of oil，this medium structure of mutation change LNAPLs frontal migration of rate，in the medium sand layered front spreading velocity is quickly and slow in the fine sand layered.In addition，when LNAPLs get to the dry and wet interface，for the capillary of pore water is stronger.

Non－aqueous phase liquids;front spreading model;porous media;infiltration and distribution

X132

A

1004－1184(2012)03－0106－04

2012－03－09

廣西自然科學(xué)基金(No.2011GXNSFB018004);廣西環(huán)境工程與保護(hù)評(píng)價(jià)重點(diǎn)實(shí)驗(yàn)室研究基金(No.桂科能0801Z007);桂林理工大學(xué)博士基金(No.002401003254)

馬騰飛(1988－)，女，河南南陽(yáng)人，碩士研究生，研究方向:水污染與水環(huán)境。